Yesterday, a colleague of mine pointed me to the article “How to fit nonlinear plant growth models and calculate growth rates: an update for ecologists” (Paine et al., 2012). It addresses a relevant topic: many plant scientists are involved in growth analyses and need to determine Absolute Growth Rates (AGRs) and Relative Growth Rates (RGRs).

The main point made by Paine et al. is that we can use the observed data to fit a growth model via nonlinear regression and then calculate model-derived growth rates together with their standard errors. In principle, the process is straightforward: we select a suitable growth model to predict biomass at any given time \(t\); the AGR at time \(t\) is the derivative of the selected growth function with respect to time, and the RGR is simply the AGR at time \(t\) divided by the biomass at that same time point.

In practice, however, implementing these calculations is not immediately clear from the manuscript—especially for us biologists. The main difficulty is that the calculations and the R code provided in the supplemental materials are algebra-based and specific to each growth model, making it hard to extract a general procedure that will always work with any dataset and any growth model.

I was intrigued by the idea that I could develop such a general procedure, and so I decided to sit down and write this post.

Before starting, I would like to offer the following suggestion: once we have a suitable growth model in mind (linear, exponential, logistic, or otherwise), we should select a parameterisation with good fitting properties (see Ratkowsky, 1990; Ritz, 2010) and that is already implemented for nonlinear regression in R, ideally with self-starting routines. The choice of parameterisation will not affect our estimates of AGRs and RGRs.

Now, let’s move on to an example, using the growth data provided in the supplemental R code of Paine et al. (2012). These data are reported and loaded in the box below, together with the ‘statforbiology’ package (which also loads the ‘drc’ package).

library(statforbiology)
dat_asymp <-
structure(list(X = c(14L, 14L, 14L, 35L, 35L, 35L, 70L, 133L, 
70L, 161L, 98L, 161L, 133L, 98L, 189L, 189L, 189L, 98L, 133L, 
161L), Y = c(0.031, 0.087, 0.09, 0.261, 0.291, 0.437, 2.104, 
2.736, 2.814, 5.832, 2.979, 9, 6.309, 3.483, 6.103, 5.655, 5.889, 
4.462, 8.392, 7.043)), .Names = c("X", "Y"), class = "data.frame", row.names = c(NA, 
-20L))
head(dat_asymp)
##    X     Y
## 1 14 0.031
## 2 14 0.087
## 3 14 0.090
## 4 35 0.261
## 5 35 0.291
## 6 35 0.437

Now, let’s imagine that we want to fit a logistic growth model. Among the available parameterisations—including the one provided in Paine et al. (2012)—we will select the following form, which has good fitting properties and is already implemented as the function L.3() in the ‘drc’ package:

\[M = \frac{d}{1 + \exp\left[ b \left( X - e\right) \right]}\]

In the above function, \(d\) is the upper asymptote, \(b\) is the slope at inflection point and \(e\) is the abscissa of the inflection point. Once we have selected this function, all we need to do is write it in R as a character string (being careful to code it correctly), compute its derivative using the D() function (note the use of parse(), which transforms the text string into an expression for D(), and deparse(), which converts the resulting derivative expression back into a text string), and then construct the ratio between the derivative and the growth function—again as a text string (using the paste() function).

fun <- "d/(1 +exp(b*(X - e)))"
derFun <- deparse(D(parse(text = fun), "X"))
rgrFun <- paste("(", derFun, ")/(", fun, ")")
derFun
## [1] "-(d * (exp(b * (X - e)) * b)/(1 + exp(b * (X - e)))^2)"
rgrFun
## [1] "( -(d * (exp(b * (X - e)) * b)/(1 + exp(b * (X - e)))^2) )/( d/(1 +exp(b*(X - e))) )"

We are ready to fit the model by using the drm() function, as shown below.

fit.logis.drm <- drm(Y ~ X, fct = L.3(),
                     data = dat_asymp)
summary(fit.logis.drm)
## 
## Model fitted: Logistic (ED50 as parameter) with lower limit fixed at 0 (3 parms)
## 
## Parameter estimates:
## 
##                Estimate Std. Error t-value   p-value    
## b:(Intercept) -0.046214   0.014877 -3.1065  0.006415 ** 
## d:(Intercept)  6.652995   0.598588 11.1145 3.221e-09 ***
## e:(Intercept) 89.621517   9.502619  9.4312 3.627e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error:
## 
##  1.279088 (17 degrees of freedom)

At this point, we choose the time points at which we want to calculate the AGRs and RGRs, and we use the gnlht()function from the ‘statforbiology’ package (Onofri, 2026). This function takes a nonlinear combination of the model parameters (specified as a text string), evaluates it, and provides standard errors obtained via the ‘delta’ method (Fox, 2019).

The calculations can be somewhat slow—gnlht()was not designed for such large numbers of evaluations—but if we restrict ourselves to a few hundred time points, the computation time remains acceptable. Please note the parameterNames argument, which matches the order of the parameters in the fitted model object with the names used in the nonlinear function func.

# Time points to estimate
nPoints <- 100
timePoints <- data.frame(X = seq(min(dat_asymp$X), max(dat_asymp$X),
                                       length = nPoints))
# Get AGR
AGRs <- gnlht(fit.logis.drm,
              func = list(derFun),
              const = timePoints,
              parameterNames = c("b", "d", "e"))
head(AGRs[,-1])
##          X    Estimate          SE  t-value   p-value
## 1 14.00000 0.008791161 0.006833979 1.286390 0.2155509
## 2 15.76768 0.009491815 0.007122149 1.332718 0.2002126
## 3 17.53535 0.010244087 0.007410021 1.382464 0.1847231
## 4 19.30303 0.011051063 0.007695797 1.435987 0.1691525
## 5 21.07071 0.011915887 0.007977492 1.493688 0.1535855
## 6 22.83838 0.012841743 0.008252938 1.556021 0.1381226
#
# Get RGRs
RGRs <- gnlht(fit.logis.drm,
              func = list(rgrFun),
              const = timePoints,
              parameterNames = c("b", "d", "e"))
head(RGRs[,-1])
##          X   Estimate         SE  t-value    p-value
## 1 14.00000 0.04485295 0.01580079 2.838652 0.01134277
## 2 15.76768 0.04474076 0.01583658 2.825152 0.01167042
## 3 17.53535 0.04461964 0.01587175 2.811261 0.01201709
## 4 19.30303 0.04448896 0.01590596 2.797000 0.01238334
## 5 21.07071 0.04434801 0.01593880 2.782393 0.01276959
## 6 22.83838 0.04419608 0.01596984 2.767471 0.01317616

Here we go: the results match those obtained in the supplemental code of Paine et al. (2012). A few remarks before closing:

1. This method is feasible with any growth model—we only need to change the name of the R fitting function. For example, we can replace the call to L.3() with any other function provided in drc or statforbiology. To obtain the corresponding equation to write as a text string, you can refer to my post on useful nonlinear regression functions at this link. Or consult my new book, available on the Springer website.
2. As noted by Paine et al. (2012), do not forget to check the basic assumptions of nonlinear regression. For example, in this dataset there are clear signs of heteroscedasticity, which can be addressed using a Transform-Both-Sides approach (see Ritz et al., 2019) or by modelling the variance structure with a Generalised Nonlinear Least Squares method (Pinheiro & Bates, 2000).

Thanks for reading—and don’t forget to check out my new book below!

Andrea Onofri
Department of Agricultural, Food and Environmental Sciences
University of Perugia (Italy)
Send comments to: [email protected]