A functional growth model with intraspecific competition applied to a sea urchin, Paracentrotus lividus

Color profile: Disabled Composite Default screen

237

A functional growth model with intraspecific competition applied to a sea urchin, Paracentrotus lividus Ph. Grosjean, Ch. Spirlet, and M. Jangoux

Abstract: A new growth model is fitted on data from reared sea urchins, Paracentrotus lividus. Quantile regressions are used instead of least-square regressions, because they are insensitive to the dimension of the measurement and accommodate more than just symmetrical distributions. Quantile regressions allow comparison of fittings on various parts of the size distributions, including large competitors versus small, inhibited animals, in the presence of a sizebased intraspecific competition. The model has functionally interpretable parameters and allows quantifying of the intensity of growth inhibition. An extension of this model, called “envelope model”, fits the whole data set at once, including size distributions. Its parameters are constrained using information about underlying biological processes involved, namely asymptotic growth with inhibition in early ages as the result of intraspecific competition, the intensity of which depends on the relative size of the individual in the cohort. The new model appears most adequate to describe growth of P. lividus and probably many other sea urchins species as well as other animals or plants. Résumé : Un nouveau modèle est ajusté sur des données de croissance de l’oursin Paracentrotus lividus en élevage. Des régressions quantiles sont utilisées à la place de la régression par les moindres carrés, parce qu’elles sont insensibles à la dimension de la mesure et restent parfaitement utilisables en présence de distributions asymétriques. Les régressions quantiles permettent une comparaison d’ajustements réalisés sur différentes parties de la distribution des tailles, y compris au niveau des plus gros et des plus petits individus, c’est-à-dire respectivement les inhibiteurs et les inhibés en présence de compétition intraspécifique basée sur la taille. Les paramètres du modèle sont fonctionnellement interprétables et permettent de quantifier l’intensité de l’inhibition. Une extension du modèle, appelée « modèle enveloppe », s’ajuste sur l’ensemble des données et prend en compte la variabilité individuelle. Ses paramètres sont contraints en utilisant les informations disponibles sur les processus biologiques sous-jacents, c’est-à-dire, une croissance asymptotique avec une inhibition dans le jeune âge due à la compétition intraspecifique et dont l’intensité dépend de la taille relative de l’individu dans la cohorte. Le nouveau modèle s’avère tout à fait adéquat pour décrire la croissance de P. lividus et probablement d’autres espèces d’oursins aussi bien que d’autres animaux ou de plantes. Grosjean et al. 246

Introduction For the last two centuries, after Malthus (1798) demonstrated the exponential nature of growth, the diversity of growth models has steadily increased (e.g., von Bertalanffy 1938; Richards 1959; Schnute 1981). However, these models compete rather than complement each other. Choosing a growth model often remains arbitrary (Fletcher 1974). Sea urchin individual growth is an example of such a problem (e.g., Ebert and Russell 1993; Gage and Tyler 1985; Lamare and Mladenov 2000). Papers on growth models applied to sea urchins compare the fit of different curves and are lim-

ited to a discussion of how they represent the growth of a mean individual. They all conclude that none of them is fully satisfactory. The difficulty with these models is that parameters are not all comparable and lack biological meaning (or lose it when applied to real data). Moreover, elaborated growth models have three or more parameters that are not independent from each other (intercorrelations) (Gallucci and Quinn 1979). Thus it is not possible to extract the estimator of one parameter from one fitting and compare it with the value obtained for the same parameter with another data set. Indeed, its value depends on the estimation of all other parameters in

Received 29 January 2002. Accepted 20 February 2003. Published on the NRC Research Press Web site at http://cjfas.nrc.ca on 4 April 2003. J16738 Ph. Grosjean.1,2 Marine Biology Laboratory (CP 160/15), Université Libre de Bruxelles, 1050 Bruxelles, Belgium. Ch. Spirlet. Marine Biology Laboratory, Université de Mons-Hainaut, 7000 Mons, Belgium. M. Jangoux. Marine Biology Laboratory (CP 160/15), Université Libre de Bruxelles, 1050 Bruxelles, Belgium, and Marine Biology Laboratory, Université de Mons-Hainaut, 7000 Mons, Belgium. 1 2

Corresponding author (e-mail: [email protected]). Present address: LOV, UMR 7093, Observatoire Océanologique, BP 28, 06234 Villefranche sur mer CEDEX, France.

Can. J. Fish. Aquat. Sci. 60: 237–246 (2003)

J:\cjfas\cjfas60\cjfas6003\F03-017.vp Tuesday, April 01, 2003 1:58:45 PM

doi: 10.1139/F03-017

© 2003 NRC Canada

Color profile: Disabled Composite Default screen

238

each respective fitting. This contrasts with linear models in which slopes can be compared and usually convey a biological or physical meaning about the relationship between the considered variables (proportionality). In an attempt to escape intercorrelation problems, some authors have tried to simplify the growth model to a single parameter (Gallucci and Quinn 1979). However, because growth is a complex emergent property of various physiological processes like feeding, digestion, respiration, etc., oversimplified models are probably of very little utility. At the other extreme, complexification, a few authors (e.g., Richards 1959; Fletcher 1974) have built up flexible and general growth models that include some other existing models as special cases. Schnute (1981) developed such a general model in which parameters carry, in theory, higher biological meanings. None of these models has proved to be fully efficient when fitting real data, partly because of the problem of intercorrelation between parameters. It is thus only possible either to carry on a global comparison of various models for the same data set or to use a single model applied to various data sets. Usually, either the R2 value or the residual sum of squares of the nonlinear leastsquare regression is used to evaluate how well a model fits the data (Cellario and Fenaux 1990; Lamare and Mladenov 2000). A visual comparison of graphs is also commonly used (Gage and Tyler 1985; Ebert and Russell 1993). These techniques are not considered rigorous by statisticians. Among all papers cited, only Lamare and Mladenov (2000) performed a complete residual analysis. In this context, a growth model, rather arbitrarily chosen, only summarizes the data cloud into a “best-fit” curve supposedly representing the growth of a mean individual. It is very difficult to use such a growth curve as a tool to investigate underlying processes (e.g., to understand the impact of a considered factor on the curve’s shape) in relation with individual variations. For instance, current models hardly cope with a superimposed effect of intraspecific competition on growth, though the latter was evidenced in sea urchins (Himmelman 1986; Grosjean et al. 1996) as in many other species (e.g., Timmons and Shelton (1980) for the largemouth bass, Micropterus salmoides, and Kautsky (1982) for the mussel Mytilus edulis). The aim of this paper is to propose a deterministic growth model that goes beyond the pure fitting of a mean curve on experimental data and provides functionally interpretable parameters, taking into account both individual variations and interactions. To reach this goal, we will cut with the traditional approach of designing growth curves from their differential equations (dynamic modelling) and fitting them using least-square regression. We will introduce fuzzy sets and quantile regression as alternative tools.

Materials and methods Experimental data: growth of a cohort of reared sea urchins A single cohort of sea urchins, Paracentrotus lividus, issued from a single artificial fertilization was reared over a period of 7 years in a controlled environment (see Grosjean et al. (1998) for a detailed description of the rearing protocol). All individuals had unlimited access to food. Echinoids

Can. J. Fish. Aquat. Sci. Vol. 60, 2003

were neither size-sorted nor individually tagged. All sea urchins in the cohort were measured every 6 months starting at 6 months old (younger echinoids are too fragile to be measured alive) until 7 years old. Size is expressed by the ambital test diameter D, which corresponds to the external diameter of the test at its largest region (the ambitus) excluding spines. D is measured with an electronic sliding caliper at the nearest 0.1 mm (Grosjean et al. 1999) and recorded into 1-mmwide size classes. Time t is measured from the metamorphosis that was artificially triggered at the same moment for all individuals using living Coralina elongata as a stimulating factor (see Grosjean et al. 1998). t is thus the age of the sea urchins after the metamorphosis, expressed in years. To estimate size and size distribution of the sea urchins just after metamorphosis, another batch of siblings is sacrificed and individual test diameter is measured by macrophotography with an image analysis system (Grosjean et al. 1996). Quantile regression Quantile regression, as defined by Koenker and Bassett (1978), is an extension to quantiles of the least-absolute deviation regression that fits a function for a median individual. In quantile regression, the objective function (called here deviance δ1) to be minimized is n

(1)

δ1 =

∑ ρτ (Di – ξ1) i=1

where ξ1 is the solution returned by the model and Di is each actual observation (diameter of a sea urchin) with i = 1,…,n observations. ρτ (u) is a piecewise linear function defined as (2)

ρτ (u) = u(τ – I(u < 0))

In eq. 2, τ is the quantile and I(u < 0) equals 1 if (u < 0) is true and 0 if it is false. Quantile τ defines the fraction of all observations that lie beneath the curve. Equation 2 is thus merely a way to apply a different weighting to positive (u > 0; weighting equals τ) and negative (u < 0; weighting equals τ – 1) residuals. If τ = 0.5, half of the observations is beneath and the other half is above the curve, both weighted the same in absolute value (*τ* = *τ – 1*), and the objective function (eq. 1) simplifies to half of the least-absolute devia1 tion δ1 = ∑*Di – ξ1*. With a different τ value than 0.5 (0 < 2 τ < 1), the fitted curve represents another fraction of the population. For instance, τ = 0.975 fits a curve that is closest to the 2.5% largest individuals because residuals above the curve are more influential (larger weighting) than those below it. They thus have to be reduced to minimize the objective function (eq. 1), and as a consequence, the fitted curve will be located nearby these large individuals. Similarly, a curve fitted using τ = 0.025 models the 2.5% smallest individuals. The surface between these two curves represents 95% of the whole batch. A third curve fitted on median individuals, using τ = 0.5, indicates asymmetry in the distribution. If this last curve is located in the middle of the two previous ones, the distribution of sizes is symmetrical. If it is closer to the lowest curve, the distribution is skewed toward small individuals. If it is closer to the highest curve, it is skewed toward large individuals. These three curves plotted together © 2003 NRC Canada

J:\cjfas\cjfas60\cjfas6003\F03-017.vp Tuesday, April 01, 2003 1:58:45 PM

Color profile: Disabled Composite Default screen

Grosjean et al.

on the same graph thus give many insights both on the relationship between the dependent (size) and the independent variable (time) and on the way in which individual variation spreads. As far as we know, only one method is currently described to fit a nonlinear model using quantile regression (Koenker and Park 1996). It uses a particular interior point algorithm (Nocedal and Wright 1999). The differential weighting of residuals above and below the curve, among other constraints, does not easily translate in most classical nonlinear regression modules (e.g., Statistica, Systat, SPSS, SAS). Special dedicated modules are required, like the “nlrq” package for R (http://cran.r-project.org). Analyses in this paper are performed using it and also some additional R code (http://www.sciviews.org/_phgrosjean/growth/index.htm) developed by the authors.

Model development Designing the model Fuzzy logic can deal efficiently with complex nonlinear problems (Zimmermann 1996; Cox 1999) and is thus another possible approach for creating growth models other than dynamic modelling (differential equations) commonly used in this field, though, it has not been employed yet. Fuzzy systems can often be formulated rather intuitively in a linguistic way (Zimmermann 1996) and we will start this way, introducing mathematics subsequently. For individual growth without competition, a trivial semantic description could be “a young, small individual gradually becomes larger with age”. In terms of fuzzy sets, this translates into a temporal transition between two sets: small (S) and large (L). Young animals belong to the S set, old ones belong to the L set, and middle-aged individuals belong partly to each of them. The degree of membership to each set depends on the amount of growth achieved and thus gradually shifts with time from set S to set L. This change is characterized by a membership function to each set, thus MS and ML, respectively. The sum of all membership values at any given time is one, because we are dealing with a single individual in its integrity. The next step is to incorporate the concept of growth inhibition to represent the effect of intraspecific competition. Grosjean et al. (1996) showed that it is a size-based competition in the case of reared sea urchins, P. lividus. Ten percent to 15% of the largest individuals (the inhibitors) in the populations grow at their maximal speed (that is, the growth speed that they would have if they were alone in the same food and environmental conditions). The growth of others depends on their relative size in the population (the smaller they are, the slower they grow). Inhibition progressively fades out when larger individuals reach their asymptotic size and are caught up with smaller ones that are still growing. According to these observations, a semantic formulation of the problem becomes “a young, small individual is potentially inhibited in its growth but still gradually reaches its maximum growth speed with age”. It can also be represented by two sets and one transition, but now set S is the minimal size with time (with maximum inhibition) and set L is the size at the age when the growth speed is maximal (with no inhibition at all). The transition is now the expression of a

239

progressive release of the inhibition, instead of a representation of the entire growth process. The difficulty resides in the proper characterization of sets and membership functions with age. Set S corresponds to the minimum possible growth, which is simply no growth at all. Thus, in set S, size remains constant at its minimum initial value just after metamorphosis (D0, see Fig. 1). In this model, no negative growth is allowed (it does not occur with animals fed ad libitum). Equation for set S is (3)

D(t) = D0

Set L describes the largest size reached by sea urchins at maximum growth speed with time. Because P. lividus has a determinate or asymptotic growth, final increase of size ∆D∞ = D∞ – D0 is finite for t → ∞. However, maximum size cannot be reached instantaneously. If the largest individuals in the actual data set are not inhibited at all, they can be used as a reference for the whole cohort to define this maximum growth curve. Supposing that a von Bertalanffy 1 curve adequately fits the growth of largest echinoids of the cohort (see next section), it is an appropriate model for set L. To make sure this curve starts from D(t = 0) = D0, we use the following parameterization: (4)

D(t) = D0 + ∆D∞(1 – e–kt)

where k quantifies the speed of maximum growth. Suppose that a good model for the membership to set L with time, ML(t), is a logistic function (as commonly done to describe transition in fuzzy sets (Cox 1999); Fig. 1): (5)

ML(t) =

1 1 + li e− ki t

where li and ki are the two parameters characterizing the way in which a possible inhibition is released with time t (i denotes parameters characterizing the inhibition). Membership to set S with time, MS(t), is complementary so that MS and ML add up to one: (6)

MS(t) = 1 – ML(t) = 1 –

1 1 + li e− ki t

The effect of intraspecific competition (or any other inhibition mechanism having a similar effect on growth) is integrated in the model as a delayed transition from set S to set L, as explicitly quantified by parameter li (the lag or position of the inflexion point in the membership curves, see Fig. 1), and ki represents the speed at which the inhibition is released with time. If li = 0, there is no inflexion point and ML(t) = 1; growth occurs at maximum speed. Thus, full-growing individuals belong to set L from the beginning and exhibit a perfect von Bertalanffy 1 growth curve. The stronger the growth inhibition, the longer other individuals remain in set S before gradually shifting to set L, and their growth curve is a more marked sigmoid. With the chosen equations, all parameters in this model carry a clear biological meaning, considering the hypotheses that were formulated to build it. Usually, fuzzy sets are manipulated using fuzzy arithmetic. The output is then “defuzzified” by one of several methods (Cox 1999) to provide a crisp number (the most probable size of an individual at a determined age). Being simple enough, © 2003 NRC Canada

J:\cjfas\cjfas60\cjfas6003\F03-017.vp Tuesday, April 01, 2003 1:58:46 PM

Color profile: Disabled Composite Default screen

240

Can. J. Fish. Aquat. Sci. Vol. 60, 2003

Fig. 1. Construction of the fuzzy growth model (diameter D versus age t). Two sets, called “large” (L, dashed curve in panel a) and “small” (S, dash-dotted curve in panel a) are used. Actual size is always located between these two curves (“fuzzy”, solid curve). Growth starts at t = 0 with the initial size D0 and is asymptotic to the maximal size D∞. Maximum size increase is thus finite (∆D∞). The membership functions (ML and MS, dashed and dash-dotted curves in panel b, respectively) model the belonging to each set with two logistic functions. If a membership is close to 1, like MS for young juveniles at t = 0 to 1 or ML for large adults at t = 8 to 10, the resulting fuzzy growth curve is close to the corresponding set. When memberships are 0.5 for each set, here at t = 4, the value returned by the fuzzy model is in the middle of values returned by both sets. Depending on the position of the inflexion point (parameter li in eq. 8) in the membership functions, the fuzzy growth curve will shift more or less rapidly from set S to set L.

the current model can also be transformed into a classical analytic equation: (7)

D(t) = MS(t)S(t) + ML(t)L(t)

which gives, after combination of eqs. 3–7 and simplification: (8)

1 − e− kt D(t) = D0 + ∆D∞ 1 + li e− ki t

In this way, the model can be treated with classical (crisp) arithmetic that offers a larger panel of statistical tools than fuzzy arithmetic. Fitting the data set Because echinoids are not tagged individually, it is not possible to track animals across measurement sets. Consequently, one will consider virtual individuals according to their relative position in the entire size distribution at each

sampled time, that is, virtual individuals corresponding to fixed quantiles in each size distribution. Virtual individuals exactly match real ones only if (i) growth is smooth and animals do not jump from one quantile to the other and (ii) mortality is randomly distributed among the cohort. This is unverifiable in this data set. However, the total number of individuals in the cohort dropped from 725 at 6 months old to 221 at 4.5 years old and to 67 at 7 years old; thus, mortality is not negligible. If virtual individuals do not match real ones, the model remains valid to describe the cohort as a whole, but one should be careful when interpreting it down to the individual level. The whole data set amounts to 4016 observations, all considered as independent by both the traditional least-square regression and quantile regression, although most individuals where measured up to 14 times (for those living to 7 years old) at 6-month intervals. There is thus probably some autocorrelation, but we will not consider it (as it is often done implicitly in similar studies). Yet, we will not attempt to diagnose the regressions using tools that rely on the independency of the error term, e.g., analysis of variance or tests for the random distribution of the residuals. Keeping these restrictions in mind, a few quantile regressions (τ = 0.95, 0.5, and 0.05) are fitted on the data set with some usual growth models (Table 1) and with a simplified four-parameter version of our new model (with D0 being fixed to zero; Table 1 and Fig. 2). It is important to note that the regressions for the three different quantiles are independent from each other. There is currently no other means to relate the different regressions than to plot them on the same graph (Fig. 2): the various equations remain independent from each other. Growth of a median individual, as well as change in the distribution of sizes with time (more particularly, the dissymmetry that appears when large individuals detach from the rest around 1–2 years old and that disappears later on), is correctly represented by all models except the logistic and von Bertalanffy 1 functions. The von Bertalanffy 1 model is adequate to describe growth of the largest animals, with τ = 0.95, but not for other quantiles. Criteria to decide which model best fits the data (deviance δ1 and visual impression on a graph) are neither rigorous nor very discriminant. Many models among the seven tested seem adequate. Our new model is adequate too, but it does not perform better that the average of the other models tested. Many models fail to extrapolate the size at metamorphosis, D(t = 0), and give sometimes totally aberrant estimates (negative values). Our new model has a parameter that explicitly expresses this size: D0. It was fixed at zero here but can be fixed to any reasonable estimate of the size at metamorphosis. The other models have to be reparameterized to quantify size around t = 0 explicitly, and even if D(t = 0) is imposed, they will not all fit the whole data set equally well. This is one advantage of a functional model over a purely empirical one: the possibility to constrain its parameters around reasonable estimates, giving knowledge and making additional data available. However, unconstrained regressions (i.e., those used to generate data in Table 1) do not trigger this advantage for our model here, except for D0, which was purposely fixed. © 2003 NRC Canada

J:\cjfas\cjfas60\cjfas6003\F03-017.vp Tuesday, April 01, 2003 1:58:46 PM

Fig. 1. C sus age t and “sma

Color profile: Disabled Composite Default screen

Grosjean et al.

241

Table 1. Results of quantile regressions with different growth models for three quantiles τ. Model τ = 0.95 New model Weibull Four-parameter logistic Logistic von Bertalanffy 1 von Bertalanffy 2 Gompertz τ = 0.50 New model Weibull Four-parameter logistic Logistic von Bertalanffy 1 von Bertalanffy 2 Gompertz τ = 0.05 New model Weibull Four-parameter logistic Logistic von Bertalanffy 1 von Bertalanffy 2 Gompertz

a

b

c

d

Deviance δ1

Fitting

D(t = 0)

65.3 65.3 63.9 60.7 66.4 62.9 61.8

0.531 0.509 0.727 1.33 0.535 0.809 6.85 × 10–2

0.92 1.07 –0.549 1.44 0.168 –0.531 0.387

1.51 69.5 –103 . . . .

2253 2232 2230 2356 2236 2243 2262

+ + + – + + +

+ + +

0.00 –4.17 –3.18 7.71 –6.26 2.67 4.24

57.2 55.9 56.4 54.6 66.3 58.8 56.7

0.612 0.199 1.02 1.39 0.360 0.688 1.36 × 10–2

5.67 1.78 1.67 2.06 0.383 –0.230 0.407

1.25 55.5 –12.5 . . . .

8135 8085 8079 8434 8668 8150 8131

+ + ++ –– –– + +

0.00 0.479 –1.81 2.92 –9.82 0.183 0.769

49.0 47.9 48.4 46.5 61.9 51.7 49.4

0.641 0.108 1.18 1.80 0.300 0.692 2.82 × 10–3

14.5 2.20 2.17 2.38 0.658 0.154 0.403

1.37 47.9 –6.08 . . . .

1894 1847 1835 1995 2067 1834 1856

+ + + – – + +

0.00 –2.96 × 10–2 –2.19 0.634 –13.5 0.139 3.62 × 10–2

+ +

+ + + – – +

Note: Weibull model, D = a – d e −bt ; four-parameter logistic, D = (a – d)/(1 + e–b(t–c)) + d; logistic, D = a/(1 + e–b(t–c)); von Bertalanffy 1, D = a(1 – t –b(t–c) e ); von Bertalanffy 2, D = a(1 – e–b(t–c))3; Gompertz, D = ab c (parameters are not all comparable between models). For comparison with the others, –bt –dt the new model (eq. 8) is D = a(1 – e )/(1 + ce ), after substituting a = 䉭D∞, b = k, c = li, and d = ki and fixing D0 = 0. “Fitting” is a visual impression of adequacy of the model on a graph (see Fig. 2 for an example). c

Constraining parameters of the model The ambital test diameter of reared sea urchins was measured on 296 individuals originating from another batch 7 days after metamorphosis (see Grosjean et al. 1996). At this age, the size is normally distributed, with a mean of 0.497 mm and a standard deviation of 0.056 mm. Because this spreading in initial sizes is negligible compared with the size scale during the whole growth process (compare 0.056 mm with the range 0.50 to 50–65 mm), one could consider an equal initial size for all echinoids just after metamorphosis (D0 = 0.50 mm). Intuitively, there should be a relationship between the curves fitted on the same data set with different quantiles, because they originate from the same conditional distribution (size distribution given the age) and also because they represent growth of virtual individuals related to the presence of other ones (inhibitors–inhibited relations). With the current model (eq. 8) and quantile regression method (eqs. 1 and 2), it is only possible to fit one curve at a time. An adaptation of both the model and the regression method, with further constraints on the parameters, are required to link these curves together. We expect that li = 0 for τ = 1, according to the hypothesis that larger animals in the batch are not inhibited at all (see above, Designing the model). We would also expect a monotonous increase of li with a decrease of τ, because fractions of smaller individuals should be more inhibited than fractions of larger ones (recall that this is a size-based com-

petition mechanism). Indeed, fitting quantile regressions for τ = 0.05 to 0.95 for every 0.05 step and plotting the corresponding estimates of parameter li in the function (1 – τ), evidence a quasilinear relationship (Fig. 3a). An acceptable expression of li across these various quantile regressions is thus (9)

li (τ) = si(1 – τ)

where si is the slope of that linear relation (the change in the degree of inhibition along the quantiles). k and ki appear to be negatively correlated, and consequently, their mean value is quite constant across τ values (Fig. 3b). A higher estimate for k is compensated by a proportionally lower estimate for ki in the same regression. Values of k and ki seem to change with τ, but this could be an artifact of the intercorrelation. Because of these considerations, we will keep relationships between k and ki and τ as simple as possible: (10)

k(τ) = constant = k ki(τ) = constant = ki

where k is probably slightly different (and lower) than ki. Finally, ∆D∞ should follow a normal distribution, as size distribution when approaching asymptotic maximum size recovers a Gaussian or near-Gaussian shape (see Fig. 2, after 4 years and also Grosjean et al. (1996)): (11)

∆D∞(τ) ~ N(µ ∆ D∞ , σ∆ D∞ ) © 2003 NRC Canada

J:\cjfas\cjfas60\cjfas6003\F03-017.vp Tuesday, April 01, 2003 1:58:47 PM

Color profile: Disabled Composite Default screen

242

Can. J. Fish. Aquat. Sci. Vol. 60, 2003

Fig. 2. Quantile regressions (D versus t) for quantiles τ = 0.05, 0.50, and 0.95 for the lower, middle, and top curves, respectively, in each graph and size distribution (a) boxplots or (b) histograms of the data set at each sampled time. The new growth model constrained to the origin (eq. 8, with D0 = 0) is used (parameter estimates and deviance δ1 are reported in Table 1). Panel b presents the fitted growth curves, whereas panel a is a diagnostic of residuals after subtraction of the median curve (∆D = D – Dmedian versus t). Extreme hinges of boxplots are quantiles τ = 0.05 and 0.95 of the observed distributions and can be compared with the curves fitted on the same quantiles. The box ranges from 1st to 3rd quarters and the division inside it is the median τ = 0.50, which can also be compared with the median curve. Changes of size spreading and asymmetries in time are emphasized in this representation.

where µ ∆D∞ is the mean and σ∆D∞ is the standard deviation of the normal distribution of ∆D∞(τ). Using eq. 11 to describe the distribution of ∆D∞(τ), the distribution of D at any time t is conditioned by the standard deviation of maximum size σ∆D∞ . As t becomes larger, the distribution of D increasingly converges to that of D∞, the value of which depends thus on τ. In such a respect, σ∆D∞ quantifies variation of individual growth in the batch, without considering the supplementary inhibition effect. Importing eqs. 9–11 into eq. 8, we obtain (12)

1 + e− kt D(t, τ) = D0 + ∆D∞(τ) 1 − si(1 − τ)e− ki t

which links curves for all quantiles 0 < τ < 1 and has six parameters: D0, µ∆ D∞ , σ∆D∞ , k, ki, and si. It incorporates individual variations into the model and is a kind of threedimensional (3D) surface that envelops data (Fig. 4). For this reason, it will be called the “envelope model”. The quantile regression method is modified as follows. Considering that every individual present in the batch is measured at each sampling time, unconditional quantiles of

Fig. 3. (a) Variation of li as a function of (1 – τ) for several quantile regressions performed separately with the new growth model constrained to the origin (eq. 8, with D0 = 0). A simple linear relationship appears suitable in a first approximation (eq. 9, with indicative values si = 12.3 and R2 = 0.980). (b) Variation of k (open triangles) and ki (solid triangles) as functions of (1 – τ). The line connecting the mean values (k + ki)/2 for each regression is almost perfectly horizontal. This illustrates how k and ki are intercorrelated: a higher value of k is associated with a proportionally lower value of ki, and vice versa.

each size distribution at a given time ti can be regarded as estimators of conditional quantiles τ at the corresponding time ti in eq. 12. Note that this is fundamentally different than the previous quantile regression method in eqs. 1 and 2 in which unconditional quantiles were not used at all. Now we consider that all observations are not independent, but we emphasize independence within each size distribution and dependence between them, which is exactly the way in which they are interdependent. Estimators of conditional quantiles τ, using unconditional quantiles noted τ$, are calculated as n( ti )

(13)

τ$ i =

∑ I(Dj (ti) < Di) j =1

n(t i)

where ti is t corresponding to the ith observation, n(ti) is the total number of individuals measured at time ti, Dj(ti) is the jth observation among all measures made at time ti, and I(u < v) returns 1 if true and 0 if false, as in eq. 2. Parameters of the envelope model to fit, ξ2 , are estimated by minimizing the following objective function δ2 : n

(14)

δ2 =

∑ Di − ξ2 (ti, τ$ i) i =1

n © 2003 NRC Canada

J:\cjfas\cjfas60\cjfas6003\F03-017.vp Tuesday, April 01, 2003 1:58:47 PM

Color profile: Disabled Composite Default screen

Grosjean et al.

Deviance δ2 is the mean absolute deviation between observed and predicted sizes for all observations with matching τ$. A robust simplex minimization algorithm is used to converge to the solution (Nelder and Mead 1965; Nocedal and Wright 1999). Fitting of the envelope model (eq. 12) by minimizing δ2 (eq. 14) emphasizes how individual variations are now included in the model itself (Fig. 4). Gain is obvious by comparing it with usual quantile regressions (Fig. 2) that led to a series of unlinked two-dimensional (2D) curves. The envelope model adequately represents size-at-age, including size spreadings and asymmetries. It is, however, a little bit less accurate for smaller sizes at 0.5 to 1.5 years old. The mean error in size prediction (δ2 ) is less than 1 mm for all quantiles and all ages.

243 Fig. 4. Envelope model (eq. 12) fitted to the whole data set. D0 = 0.500 (imposed) and µ ∆D ∞ = 57.4, σ∆D ∞ = 3.91, k = 0.555, ki = 1.33, and si = 12.3 (calculated); deviance δ2 = 0.965. All quantiles are fitted at once and the model is a three-dimensional surface around the size distributions. The solid curves are quantiles τ = 0.05, 0.50, and 0.95 (from bottom to top in each graph) extracted from the model for comparison with Fig. 2. Predicted size distributions at each sampled time (shaded curves) are mapped over the observed size distribution histograms and corresponding boxplots (shaded boxplots) are plotted behind and slightly lagged from those of the observed data (open boxplots). Globally, the model respects size distributions, that is, both spreading and asymmetries at all times but for the smallest sizes at 0.5– 1.5 years old.

Discussion Which regression method for growth? By far the most widespread method to fit a curve is the least-square regression with one of the many minimization algorithms available (simplex, (quasi-)Newton, etc.; Draper and Smith 1998; Nocedal and Wright 1999). The algorithm finds the combination of values for the various parameters in the model (the solution) that leads to a minimal value for the objective function, which is here the sum of the square of the residuals (that is, the sum of squared distances between observed values for the dependent variable and values predicted by the model at the same levels of the independent variables). Least-square regression has many advantages over other methods. In particular, when partial first (gradient matrix) and second (Hessian matrix) derivatives of the function are calculable for each parameter, convergence through a solution is accelerated and can be verified (at least for a local solution; Nocedal and Wright 1999). In the counterpart, that regression supposes that the fluctuations around the model (called the error term) are additive, independent, normally distributed, and with a constant standard deviation (heteroscedasticity). It is also very sensitive to outliers because it uses the squared residuals. Those constraints, even if not strictly met every time, particularly in many nonlinear phenomena like growth, appear to be of minor importance for many authors. Indeed, outliers are eliminated, or weighing methods are applied to limit their impact. Yet, there are two additional arguments against the use of least-square regression with growth models. First, it can only represent a mean effect, by definition, and this mean effect implies losing any information on the individual variations of growth. Second, least-square regression is sensitive to the dimension of the size measurement: models fitted on length data (dimension 1) never represent the same mean individual as those fitted on weight or volume data (dimension 3). There is more than just a subtle difference. So, what is the exact meaning of a “mean-sized” individual if it depends on the way in which size was measured? In the contrary, in the absence of measurement errors or individual deviations from allometric relationships, median (as any other quantile) always represents the same individual, no matter the dimension of the measurement used to quantify size. Because it can also fit different parts of the size distribution, there are

many reasons to use quantile regression instead of leastsquare regression for fitting growth models. We showed that nonlinear quantile regression using interior point algorithm proposed by Koenker and Park (1996) is usable with many growth models. Modelling variations of growth When a size-based intraspecific competition occurs, the largest animals are inhibitors, whereas the smallest ones represent the most inhibited fraction. Growth curves fitted on large and small individuals thus contain information about the impact of the intraspecific competition (by comparison). When individuals are not tracked from one measurement to the other in the population, this information is only available using two different quantile regressions, based on large and small quantiles. For more descriptive fittings, we introduced the triple τ = 0.95, 0.50, 0.05 quantile regression representation as an informative summary of the data together with a representation of residuals and size distributions using boxplots around the median model. Ten percent is a commonly used critical level in statistics. The curves for τ = 0.95 and t = 0.05 create © 2003 NRC Canada

J:\cjfas\cjfas60\cjfas6003\F03-017.vp Tuesday, April 01, 2003 1:58:48 PM

Color profile: Disabled Composite Default screen

244

a kind of two-tailed 10% nonparametric conditional confidence interval of the data set: 90% of data are inside this interval. Additionally, a τ = 0.50 median curve visualizes possible asymmetry and its change with time. When a representative sample of the whole size distribution at each time is available, typically with n > 50, histograms of size distributions can be superimposed on the graph to show how well quantile regressions match them. This is even better evidenced by boxplots around the median curve (∆D in the graphs). As useful as they are, usual quantile regressions always give curves independently fitted for various quantiles. Both the model and the fitting method must be adapted to link these curves into a single model. For the cases in which representative samples are measured at each time interval, we have proposed a modified method to use with envelope models, that is, models of the form Y = f(t, τ) as eq. 12. These models calculate a 3D surface enveloping data. Curves for all quantiles 0 < τ < 1 are calculated at once. They contain most of the information in the initial data set, including individual variability. Functional analysis of constrained parameters in the envelope model Flexible, unconstrained envelope models can be incredibly complex, with dozens of parameters. They are impossible to fit in practice. Constraining parameters as we did in eqs. 9– 12 leads to a double benefit. First, it keeps the model simple. In the present case, we started with a five-parameter, unconstrained classical model (eq. 8) and we ended with a sixparameter, constrained envelope model (eq. 12). Yet, the latter contains much more information than the former. Second, if constraints are formulated according to some knowledge about the underlying phenomenon (value of the intercept corresponding to actual initial size) or to some reasonable hypotheses (relation between li and τ from information in the literature), parameters remain meaningful in the fitted model. Constraining the model to the origin is very easy, in theory. Most models in our study and also eq. 8 have a free intercept. It could be formally expressed (D0 in eq. 8) or hidden in the parameterization (a(1 – ebc) in von Bertalanffy 1 or a – d in Weibull models in this study). Unconstrained intercept means that the parameter representing size when the growth process initiates is estimated at the same time as all others and is thus influenced by their values (intercorrelation). In real life, the initial size can influence following growth (for some experimental studies on P. lividus, see Vaïtilingon et al. (2001)). We believe that a meaningful model should follow the same logic: initial size is fixed first and parameters that characterize growth are estimated afterwards. This is achieved with the new model by fixing D0. Of course, initial size just after metamorphosis is the result of another growth process during larval life but the model describes postmetamorphic growth, not larval growth. If a model does not fit correctly after fixing the origin, it means that it is not adapted to describe growth in this case. The only good reason to avoid constraint is when time or size or both are unknown at the origin of the growth process. It is unfortunately common with data collected in the field, when it is not possible to estimate age accurately (Ebert 1998; Russell and Meredith 2000). In this case, only relative

Can. J. Fish. Aquat. Sci. Vol. 60, 2003

growth can be studied and the problem of origin is thus eliminated de facto. However, the model must be reworked to fit relative growth data. At the other extreme of the growth process, a single parameter characterizes its completion when growth is asymptotic in all models (parameter a or ∆D∞). Several authors questioned whether asymptotic growth is a biological reality or just a mathematical artifact. Ricker (1979) wrote a section entitled “Asymptotic growth: is it real?” in a chapter of a book; Knight (1968) devoted a whole article to demonstrating that it is biological nonsense. Some models with infinite growth appeared (for instance, Tanaka 1982). They were also tested on sea urchins (Ebert and Russell 1993; Ebert 1998). Some P. lividus were reared in our installations for 15 years. They reached their maximum size at 4–5 years old. They thus displayed exactly the same size for almost 10 years, indicating asymptotic growth for this species. For other species, where no plateau is observed, a lifetime could be simply too short to reach it. Yet, it is then impossible to tell if growth is determinate or indeterminate. Anyway, if maximum size is not actually reached, it is very difficult to estimate the corresponding parameter in the model. We constrained ∆D∞ to be normally distributed in the envelope model (eq. 12). It is in agreement with the analysis of size distributions for full-grown animals (Grosjean et al. 1996; current data set). It is also a consequence of the genetic homogeneity of the batch as all individuals are siblings in this data set. As a consequence of fixing k (eq. 10), ∆D∞ is the only parameter to contain information on relative growth potential of the individuals among the cohort in eq. 12. The kinetic parameter k could be viewed as environment dependent (temperature, food, water quality, etc.). Because these are the same for all animals as they are in the same aquarium, they are fed ad libitum and have access to the food in the same way, it appears logical to fix k. Fixing ki is motivated by a similar reason: we want it to express one global aspect of the inhibition. When homogeneous batches of animals of the same age and same genetic origin are reared together, speed at which inhibition is released is supposed to be about the same for all individuals. In this way, only li quantifies changes between virtual individuals (inhibitors versus inhibited). Of course, many other variants are possible, but at the cost of an increasing complexity of the model. Even in the case of a single cohort reared in constant artificial conditions, eq. 9 seems to be a very simplified relationship between li and τ. There are some visually detectable departures for linearity with extreme quantiles in Fig. 3a. This could probably explain why the model is less accurate for smallest sizes at 0.5 to 1.5 years old. Again, the simplest possible model was presented, but many other variations can be conceived for more accurate results. Relations between li and τ could be even more complex in other circumstances: different rearing methods, large and small animals of different ages and (or) genetic origins maintained together, periodic size sorting in culture, etc. Profiles of li in the function τ must be studied in each particular case. It is also probably very different in the field, or for other species. The model still needs to be reformulated and tested before being used with field-collected data (mainly cohort separation or mark– recapture) to confirm it. © 2003 NRC Canada

J:\cjfas\cjfas60\cjfas6003\F03-017.vp Tuesday, April 01, 2003 1:58:48 PM

Color profile: Disabled Composite Default screen

Grosjean et al.

An interesting potential of this model, because of variations in the relations between li and τ, is the possibility of predicting growth of the remaining fraction after elimination of the largest animals (fisheries or harvesting of largest fraction in aquaculture). Virtual individuals just below minimum harvesting size suddenly become the largest fraction and will exhibit a very rapid catch-up growth to reach the maximum growth speed curve (as evidenced by Grosjean et al. 1996). However, this goes far beyond the scope of this paper. Not considering individual variation and asymmetries in the size distributions could lead to the rejection of the von Bertalanffy model (Sainsbury 1980). In the case of P. lividus, Cellario and Fenaux (1990) for reared and Turon et al. (1995) for wild populations both rejected the von Bertalanffy model in favor of the Gompertz curve. We would reach the same conclusion if we considered only median quantile regression with τ = 0.5 in the present study. However, taking intraspecific competition into account using the new growth model leads to a different conclusion when inhibition is eliminated (the new model appears very suitable and the von Bertalanffy 1 model fits reasonably quantile τ = 0.95). Thus, the present study rehabilitates von Bertalanffy’s theory (von Bertalanffy 1957) that was too often rejected after observation of a sigmoidal growth (and now we know it could result from an inhibition). In conclusion, the new growth model with intraspecific competition is a very flexible one. It can accommodate different situations and has meaningful parameters that allow exploring and quantifying various aspects of growth. Using a quantile regression method, modified for envelope modelling, and constraining parameters ensures the meaning of the latter is saved into the fitted model. It takes also individual variability into account. This is particularly useful to model growth of P. lividus and probably of many other sea urchin species and other animals or plants. It allows for quantifying the degree of inhibition in the case of size-based intraspecific competition. Other aspects of this model are original. It was designed by defuzzifying a fuzzy model, whereas most of the other growth models were built from their differential equations. Fuzzy sets and transitions (membership functions) can be combined in countless ways to create many other similar models. We propose to call the family of models designed this way “fuzzy-remnant functions”. From their fuzzy origin, they keep nothing in appearance, but the biological meaning of their parameters is still there. Recalling the fuzzy model from which they originated, one has a much clearer idea of how various components (sets and membership functions) interact to produce the final result. Fuzzy logic is closer to the way in which the human brain conceptualizes complex objects. Statistical tools handle defuzzified analytic functions more conveniently. By their bivalence, fuzzy-remnant functions promise to be powerful tools for modelling complex nonlinear phenomena, e.g., growth, in a functional way.

Acknowledgments We thank the Centre Régional d’Etudes Côtières and the University of Caen for their contribution in building a specific sea urchin rearing facility. We are grateful to Didier Bucaille, who performed the laboratory measurements. We

245

thank also Michael Russell for constructive criticisms and Jean-Pierre Van Noppen for proofreading the manuscript. This study was conducted in the framework of the European Contracts AQ2.530 “Sea urchins cultivation” and FAIRCT96-1623 “Biology of sea urchins under intensive cultivation (closed cycle echiniculture)”. This is a contribution of the Centre Interuniversitaire de Biologie Marine.

References Cellario, Ch., and Fenaux, L. 1990. Paracentrotus lividus (Lamarck) in culture (larval and benthic phases): parameters of growth observed during two years following metamorphosis. Aquaculture, 84: 173–188. Cox, E. 1999. The fuzzy systems handbook. 2nd ed. Academic Press, San Diego, Calif. Draper, N.R., and Smith, H. 1998. Applied regression analysis. 3rd ed. Wiley & Sons, New York. Ebert, T.A. 1998. An analysis of the importance of Allee effects in management of the red sea urchin Strongylocentrotus franciscanus. In Echinoderms. Edited by R. Mooi and M. Telford. Balkema Publishers, Rotterdam. pp. 619–627. Ebert, T.A., and Russell, M.P. 1993. Growth and mortality of subtidal red sea urchins (Strongylocentrotus franciscanus) at San Nicolas, California, U.S.A.: problems with models. Mar. Biol. 117(1): 79–89. Fletcher, R.I. 1974. The quadratic law of damped exponential growth. Biometrics, 30: 111–124. Gage, J.D., and Tyler, P.A. 1985. Growth and recruitment of the deep-sea urchin Echinus affinis. Mar. Biol. 90: 41–53. Gallucci, V.F., and Quinn, T.J. 1979. Reparameterizing, fitting, and testing a simple growth model. Trans. Am. Fish. Soc. 108: 14–25. Grosjean, Ph., Spirlet, Ch., and Jangoux, M. 1996. Experimental study of growth in the echinoid Paracentrotus lividus (Lamarck, 1816) (Echinodermata). J. Exp. Mar. Biol. Ecol. 201: 173–184. Grosjean, Ph., Spirlet, Ch., Gosselin, P., Vaïtilingon, D., and Jangoux, M. 1998. Land-based closed cycle echiniculture of Paracentrotus lividus (Lamarck) (Echinoidea: Echinodermata): a long-term experiment at a pilot scale. J. Shellfish Res. 17: 1523–1531. Grosjean, Ph., Spirlet, Ch., and Jangoux, M. 1999. Comparison of three body-size measurements for echinoids. In Echinoderm research 1998. Edited by M.D. Candia Carnevali and F. Bonasoro. Balkema Publishers, Rotterdam. pp. 31–35. Himmelman, J.H. 1986. Population biology of green sea urchins on rocky barrens. Mar. Ecol. Prog. Ser. 33: 295–306. Kautsky, N. 1982. Growth and size structure in a Baltic Mytilus edulis population. Mar. Biol. 68: 117–133. Knight, W. 1968. Asymptotic growth: an example of non-sense disguised as mathematics. J. Fish. Res. Board Can. 25: 1303–1307. Koenker, R., and Bassett, G. 1978. Regression quantiles. Econometrica, 46: 33–50. Koenker, R., and Park, B.J. 1996. An interior point algorithm for nonlinear quantile regression. J. Econometrics, 71(1–2): 265–283. Lamare, M.D., and Mladenov, P.V. 2000. Modelling somatic growth in the sea urchin Evechinus chloroticus (Echinoidea: Echinometridae). J. Exp. Mar. Biol. Ecol. 243: 17–43. Malthus, T.R. 1798. An essay on the principal of population. Reedition (1970). Penguin Books, New York. Nelder, J.A., and Mead, R. 1965. A simplex algorithm for function minimization. Comput. J. 7: 308–313. Nocedal, J., and Wright, S.J. 1999. Numerical optimization. SpringerVerlag, New York. © 2003 NRC Canada

J:\cjfas\cjfas60\cjfas6003\F03-017.vp Tuesday, April 01, 2003 1:58:48 PM

Color profile: Disabled Composite Default screen

246 Richards, F.J. 1959. A flexible growth function for empirical use. J. Exp. Bot. 10(29): 290–300. Ricker, E. 1979. Growth rates and models. In Fish physiology. Vol. 8. Academic Press, New York. Russell, M.P., and Meredith, R.W. 2000. Natural growth lines in echinoid ossicles are not reliable indicators of age: a test using Strongylocentrotus droebachiensis. Invert. Biol. 119: 410–420. Sainsbury, K.J. 1980. Effect of individual variability on the von Bertalanffy growth equation. Can. J. Fish. Aquat. Sci. 37: 241–247. Schnute, J. 1981. A versatile growth model with statistically stable parameters. Can. J. Fish. Aquat. Sci. 38: 1128–1140. Tanaka, M. 1982. A new growth curve which expresses infinite increase. Publ. Amakusa Mar. Biol. Lab. 6: 167–177. Timmons, T.J., and Shelton, W.L. 1980. Differential growth of largemouth bass in West Point reservoir, Alabama–Georgia. Trans. Am. Fish. Soc. 109: 176–186.

Can. J. Fish. Aquat. Sci. Vol. 60, 2003 Turon, X., Giribet, G., Lopez, S., and Palacin, C. 1995. Growth and population structure of Paracentrotus lividus (Echinodermata: Echinoidea) in two contrasting habitats. Mar. Ecol. Prog. Ser. 122: 193–204. Vaïtilingon, D., Morgan, R., Grosjean, Ph., Gosselin, P., and Jangoux, M. 2001. Influence of delayed metamorphosis and food intake on the perimetamorphic period of the echinoid Paracentrotus lividus. J. Exp. Mar. Biol. Ecol. 262(1): 41–60. von Bertalanffy, L. 1938. A quantitative theory of organic growth (inquiries of growth laws II). Hum. Biol. 10(2): 181–213. von Bertalanffy, L. 1957. Quantitative laws in metabolism and growth. Quart. Rev. Biol. 32(3): 217–231. Zimmermann, H.J. 1996. Fuzzy set theory and its applications. 3rd ed. Kluwer Academic, Boston.

© 2003 NRC Canada

J:\cjfas\cjfas60\cjfas6003\F03-017.vp Tuesday, April 01, 2003 1:58:48 PM