5  Nonlinear regressions

NoteIn this chapter

In this chapter, you will learn

Workflow steps (Chapter 12): step 3 (a-priori model), step 5 (residual checks), step 6 (inference).

Prerequisites: the linear model (Chapter 2).

5.1 What is a nonlinear regression?

The most important thing first: the main distinction between a linear and a nonlinear regression is NOT whether you fit a straight line. The defining criterion is whether the model is linear in the parameters. Quadratic, cubic and other polynomial functional forms (although sometimes also referred to as “polynomial regressions”) are all still linear regressions, and so are transformations such as \(\log(x)\) or spline bases.

The reason is that as long as the function of the mean is linear in the parameters - for example a polynomial of the form

\[ y = a_0 + a_1 \cdot x + a_2 \cdot x^2 + ...\] you are considering linear effects (\(a_0, a_1, a_2\)), where the predictor just happens to be \(x^2\), \(\sqrt{x}\) or something like that. You can see this by noting that instead of writing the \(x^2\) in the regression formula, you could pre-calculate it in the data as \(x_2 = x*x\) and it would just look like a multiple regression.

Because the linear regression is much better supported in R (as in other statistical packages), it is usually preferable to consider if the function you want to fit can be made linear in the parameters by transforming y or x (e.g. by regressing against log(y), logit(y)). Note, however, that transforming the response changes the error structure — for example, regressing against log(y) with normal residuals implies multiplicative (log-normal) errors on the original scale, which is a different model from an NLS fit of the untransformed function, which assumes additive normal errors.

However, there are a number of functions that cannot be expressed as a linear regression. This usually means that the parameters of your model are somehow inside a more complex function that cannot be made linear in the parameters. Typical examples in ecology are growth or density dependence models, or in biology certain reaction or dose-response curves. In these and other cases, you will have to run a nonlinear regression. The simplest case is a nonlinear regression with normal residuals, which is known as nls (nonlinear least squares).

5.2 Fitting a nonlinear regression using nls

Here, I will use the Michaelis-Menten model as an example. In this example, we are interested in how the reaction rate depends on the substrate concentration. The data we have are paired measurements of the reaction rate (\(v\)) and substrate concentration (\(S\))

library(nlstools)

data(vmkm)
plot(vmkm)

From the MM model, we can derive the evolution of the reaction rate (v) as a function of the concentration of substrate (S) (For details, see ?nlstools::michaelismodels ). The function we want to fit is of the form:

\[ v = \frac{S}{S + K_m} \cdot V_{max}\]

In base R, the function nls() (package stats) fits such models: you specify the model as a formula with named parameters and provide starting values for them:

mm <- nls(v ~ S/(S + Km) * Vmax,data = vmkm, 
          start = list(Km=1,Vmax=1))

Alternatively, the helpful nlstools package already provides many functions that you may want to fit, including this particular function under the name ‘michaelis’. We have already loaded this package, so you can just write

mm <- nls(michaelis,data = vmkm, start = list(Km=1,Vmax=1))

Note that unlike in an lm, the nls function requires you to provide starting values for the parameters. This helps the optimizer to find the NLS solution (which equals the MLE under i.i.d. normal residuals), which can sometimes be problematic. In this case, you may also want to consider changing the algorithm and some of its settings (via algorithm and control, see ?nls). Starting values can be read off the plot (e.g. Vmax at least the largest observed v, Km about the S at which v reaches half of Vmax); if they are far off, nls typically fails with errors such as “singular gradient”. For many standard functions, self-starting models find the starting values automatically, e.g. nls(v ~ SSmicmen(S, Vmax, Km), data = vmkm).

The model provides a standard summary table.

summary(mm)

Formula: v ~ S/(S + Km) * Vmax

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
Km     2.5861     0.3704   6.981 8.94e-07 ***
Vmax   1.5460     0.1454  10.635 1.11e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02572 on 20 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 2.477e-06

The nlstools package has a slightly extended summary table that provides a few useful additional information

overview(mm)

------
Formula: v ~ S/(S + Km) * Vmax

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
Km     2.5861     0.3704   6.981 8.94e-07 ***
Vmax   1.5460     0.1454  10.635 1.11e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02572 on 20 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 2.477e-06

------
Residual sum of squares: 0.0132 

------
t-based confidence interval:
         2.5%    97.5%
Km   1.813418 3.358860
Vmax 1.242790 1.849303

------
Correlation matrix:
            Km      Vmax
Km   1.0000000 0.9917486
Vmax 0.9917486 1.0000000

A plot can be generated with

nlstools::plotfit(mm, smooth = TRUE)

5.3 Residual checks

Residual plots can be created using the nlstools package

res <- nlsResiduals(mm)
plot(res, which = 0)

Of course, all discussions about residuals from the section on residual checks in the linear regression chapter (Section 2.4) also apply here, so we are checking for iid normal residuals.

5.4 Uncertainties and CIs on parameters

In nonlinear models, the sampling distribution of the MLE will often not be well approximated by a multivariate normal distribution, which means that approximation errors can be large when calculating CIs based on the variance-covariance matrix as done in the regression table (which assumes that the likelihood surface is multivariate normal). The nlstools package allows you to calculate CIs based the normal approximation (asymptotic, corresponds to the se in the regression table) or on profiling.

nlstools::confint2(mm)  # asymptotic (normal approximation) CI
        2.5 %   97.5 %
Km   1.813418 3.358860
Vmax 1.242790 1.849303
confint(mm)  # profile-likelihood CI, preferable for small samples
         2.5%    97.5%
Km   1.940953 3.609793
Vmax 1.290728 1.944217

We can also visualize the RSS contour (under normal errors, minimizing the RSS is equivalent to maximizing the likelihood, so this also traces the likelihood surface)

cont <- nlsContourRSS(mm)
0%100%
 RSS contour surface array returned 
plot(cont)

Alternatively (if you assume the normal approximation is bad), you can bootstrap the model, which can also be used to generate CIs on predictions (nlsBoot resamples the residuals, i.e. it assumes i.i.d. errors)

mmboot <- nlstools::nlsBoot(mm, niter = 1000)
plot(mmboot)

summary(mmboot)

------
Bootstrap statistics
     Estimate Std. error
Km   2.637783  0.3730425
Vmax 1.565863  0.1455113

------
Median of bootstrap estimates and percentile confidence intervals
       Median     2.5%    97.5%
Km   2.602323 2.049972 3.540909
Vmax 1.547969 1.330331 1.896717

In this case, the bootstrap and the profile intervals are similar to each other and slightly shifted to the right of the symmetric asymptotic interval for Km (larger upper limit), so the normal approximation is acceptable here, but somewhat optimistic for Km.

5.5 Uncertainties and CIs on predictions and derived quantities

In principle, predictions on an nls object (also for new data) can be made with the predict.nls function

predict(mm)

however, this function does not calculate confidence intervals. You can add those using the propagate package, which propagates uncertainties through the nonlinear model based on MC samples or Taylor expansion based on the normal approximation of the likelihood.

library(propagate)
out = propagate::predictNLS(mm, nsim = 10000) # nsim must be set much higher for production

For complex nonlinear models with few data, the normal approximation will not be a good one. In this case, you can resort to the bootstrap, using the general boot function described in the section on the bootstrap in the linear regression chapter (Section 2.12). This can also be applied to any derived quantities (i.e. quantities that are calculated from a function of the model parameters such as an LD50 value). However, note that wherever possible, it is easier to fit those derived quantities by reshuffling the regression formula.

The most complete solution to calculating confidence intervals and uncertainties and predictions on derived quantities would be to fit the model fully Bayesian. This can be done, for example, using the BayesianTools R package (see the package vignette about how to do this, and the chapter on posterior estimation in Hartig (2024)).

Hartig, F. (2024) “Introduction to Bayesian statistics with R.” Available at: https://florianhartig.github.io/LearningBayes/.

5.6 Interactions in NLS

What about interactions, i.e. one variable influencing the parameter value of another variable? Numeric interactions in nls will just be coded as multiplications in the formula. Factor interactions are not explicitly supported by nls but you can create them via dummy variables (i.e. introduce a variable per factor level and put a 1 if the observation corresponds to the factor level, zero otherwise). The estimated parameters then correspond to treatment contrasts (see Section 2.6.1)

Here I created an example with the previous dataset, which I randomly split between wild type and mutant observations. I want to know if the parameter Km differs between wild type and mutant.

vmkm1 = vmkm
set.seed(1)
vmkm1$type = sample(c(0,1), nrow(vmkm), replace = TRUE) # wild type and mutant
head(vmkm1)
    S    v type
1 0.3 0.17    0
2 0.3 0.15    1
3 0.4 0.21    0
4 0.4 0.23    0
5 0.5 0.26    1
6 0.5 0.23    0

Now I can fit the model

mm1 <- nls(v ~ S/(S + Km + type*DeltaMutantKm) * Vmax,data = vmkm1, start = list(Km=1,Vmax=1, DeltaMutantKm = 0)) 
summary(mm1)

Formula: v ~ S/(S + Km + type * DeltaMutantKm) * Vmax

Parameters:
              Estimate Std. Error t value Pr(>|t|)    
Km              2.5647     0.3750    6.84 1.58e-06 ***
Vmax            1.5456     0.1477   10.47 2.51e-09 ***
DeltaMutantKm   0.0532     0.1004    0.53    0.603    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.0262 on 19 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 5.254e-06

We see here that there is no significant difference between wild type and mutant, which is reasonable as I just randomly put the observations in one or the other type. Of course, you can extend this procedure to additional variables.

There is no automatic plotting function for this, so you will have to plot by hand:

dat = data.frame(S = seq(0,2,length.out = 50), type = 0)
plot(dat$S, predict(mm1, newdata = dat), type = "l", ylim = c(0, 0.75), xlab = "S", ylab = "v")
dat$type = 1
lines(dat$S, predict(mm1, newdata = dat), type = "l", col = "red")
points(vmkm1$S, vmkm1$v, col = vmkm1$type +1)

5.7 Model selection

You can use all model selection methods described in much more detail in the chapter on model selection (Chapter 8). Just very briefly, if you want to compare, for example, if the model with the interaction as fitted in Section 5.6 is better than the null model that doesn’t distinguish between types, we can either do a nested model test with anova() (an F-test, i.e. a variant of the likelihood ratio test; for nls, its reference distribution is only approximate, see Tip 3.1)

anova(mm, mm1)
Analysis of Variance Table

Model 1: v ~ S/(S + Km) * Vmax
Model 2: v ~ S/(S + Km + type * DeltaMutantKm) * Vmax
  Res.Df Res.Sum Sq Df    Sum Sq F value Pr(>F)
1     20   0.013234                            
2     19   0.013039  1 0.0001953  0.2846 0.5999

or an AIC comparison

AIC(mm)
[1] -94.71828
AIC(mm1)
[1] -93.04535

Both of which tell us the same thing our previous result in Section 5.6 told us, which is that the model that fits separate Km parameters for wild type and mutant is not better than the model that assumes that they have the same functional relationship.

Note that, similar to what is discussed in the ANOVA chapter, the model selection procedure is often helpful to decide whether there is a difference between wild type and mutant AT ALL, whereas the procedure via calculating treatment contrasts on the interactions (previous subsection) is useful to decide for each parameter separately if there are significant differences, and if so, how strong.

5.8 Mixed effects

If you want to add mixed effects to your model (which can often make sense, e.g. consider you would have data as above, but from 3 different experiments), we can use the nlme function in the package nlme. To create a simple example, I add a variable group to the data we used before, assuming that the data came out of 3 different experiments

library(nlme)
vmkm2 = vmkm
vmkm2$group = sample(c(1,2,3), nrow(vmkm), replace = T)
vmkm2 <- groupedData(v ~ 1|group, vmkm2) # transferring to a grouped data object recommended for nlme
head(vmkm2)
Grouped Data: v ~ 1 | group
    S    v group
1 0.3 0.17     1
2 0.3 0.15     1
3 0.4 0.21     1
4 0.4 0.23     1
5 0.5 0.26     2
6 0.5 0.23     1

We now want to assume that the parameter Km may be different in the 3 experiments. In this case, we use nlme, provide the data, formula and start values as before, and we can decide which of the parameters should be treated as fixed and which as random effects. In nlme, random effects are placed on parameters of the nonlinear function: random = Km ~ 1 | group lets Km vary between groups, with deviations drawn from a normal distribution.

mm2 = nlme(model = michaelis,data = vmkm2, 
           fixed = Vmax + Km ~ 1, random = Km ~ 1|group, 
          start = c(Km=1, Vmax=1))

summary(mm2)
Nonlinear mixed-effects model fit by maximum likelihood
  Model: michaelis 
  Data: vmkm2 
        AIC       BIC   logLik
  -92.71828 -88.35411 50.35914

Random effects:
 Formula: Km ~ 1 | group
                  Km   Residual
StdDev: 3.176079e-06 0.02452676

Fixed effects:  Vmax + Km ~ 1 
        Value Std.Error DF   t-value p-value
Vmax 1.546013 0.1453726 18 10.634833       0
Km   2.586054 0.3704204 18  6.981403       0
 Correlation: 
   Vmax 
Km 0.992

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.84557904 -0.38103071 -0.01901336  0.57185719  2.39956861 

Number of Observations: 22
Number of Groups: 3

The estimated SD of Km between groups is close to zero, as expected for randomly assigned groups. Variance and correlation structures can be added in the same way, via gnls() or the weights and correlation arguments of nlme(), see Chapter 10 and Chapter 11.

5.9 Take-home messages

Tip
  1. A regression is nonlinear if the model is not linear in the parameters, not if the curve is not straight: polynomials, log(x) and splines are still linear regressions. If the function can be made linear by transforming x or y, consider lm (better supported in R), but remember that transforming the response changes the error structure.
  2. nls() needs starting values for the parameters. If the fit fails or converges badly, change the starting values, the algorithm or the control settings.
  3. nls is the maximum-likelihood estimate only for i.i.d. normal residuals. Check the residuals exactly as in a linear regression.
  4. CIs from the regression table rely on a multivariate normal approximation, which can be poor in nonlinear models. Use profile CIs (confint()) or the bootstrap, and treat large differences to the table CIs as a sign that the approximation fails. For derived quantities (e.g. an LD50), use the bootstrap or reparameterize the formula so that the quantity is a parameter.
  5. To compare groups, code factors as dummy variables inside the formula. The effects on the parameters show which parameter differs; model selection (LRT or AIC) decides whether there is a difference at all.