OXITEC: Parameter inferrence - nonlinear regression [R]

With nonlinear regression there is an iterative process, we provide the initial guess; and the algorithm does then adjusts and improves the function step by step, until it hopefully converges on a least square solution. From my experience, providing initial guesses may be quite troublesome. Therefore, it is very convenient to use R functions including the appropriate self-starting routine, which can greatly simplify the fitting process.

The packages that you can look for self-starter formulas can be found in the following packages:

• drc - used with the drm() nonlinear regression facility
• nlme - used with nls(), nlsList(), nlme() nonlinear regression facilities
• aomisc - package developed by the author

if(!require('drc')){install.packages('drc');library('drc')}else{library('drc')}
if(!require('nlme')){install.packages('nlme');library('nlme')}else{library('nlme')}
if(!require('aomisc')){install.packages('aomisc');library('aomisc')}else{library('aomisc')}

Shapes

Polynomials

They arehey are actually the most flexible tool to descibe biological processes. They are simple and, although *curvilinear, they are linear in the parameters and can be fitted by using linear regression. One disadvantage is that they cannot describe asymptotic processes, which are very common in biology. Polynomialas are also prone to overfitting, we can really get a nice fit to the data which often does not reflect a sort of biological realism.

• Nowdays (becaus the wide availability of nonlinear regression algorithms) the use of polynomials has decreased
• Linear and quadratic polynomials are sitll used, but higher order polynomials are very rarely seen in practice

Polynomials: Linear equations

\[Y = b_0 + b_1 X\]

Polynomials: Quadratic equations

The equation is:

\[Y = b_0 + b_1X + b_2X^2\]

It is interesting to see that the first derivative is \(dY/dx\):

D(expression(a + b*X + c*X^2), "X")

b + c * (2 * X)

Polynomials: Polynomial fitting in R

Polynomials in R are fit by using the linear model function lm(). Although this is not efficient, in a couple of cases I found myself in the need of fitting a polynomial by using the nls() or drm() functions.

Concave curves

Exponential

The exponential equation describes an increasing/decreasing trend, with constant relative rate. The most common parameterisation is: \[Y = ae^{kx}\] Other possible parameterisations are \[Y = ab^X=e^{d+kX}\] First derivative of this function is \(dY/dX\):

D( expression(a * exp(k * X)), "X")

a * (exp(k * X) * k)

The exponential curve is normally used to describe the growth of a population in unlimiting environmental conditions, or to describe the degradation of xenobiotics in the environment (first-order degradation of kinetic).

The exponential function is nonlinear in k and needs to be fitted by using nls(), and drm().

Let’s see an example of fitting using the degradation dataset

library(aomisc)

Loading required package: drc
Loading required package: MASS
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

'drc' has been loaded.

Please cite R and 'drc' if used for a publication,
for references type 'citation()' and 'citation('drc')'.


Attaching package: ‘drc’

The following objects are masked from ‘package:stats’:

    gaussian, getInitial

Loading required package: plyr
Loading required package: car
Loading required package: carData
Loading required package: multcompView

data(degradation)
head(degradation)

summary(degradation)

      Time           Conc        
 Min.   : 0.0   Min.   :  0.200  
 1st Qu.:17.5   1st Qu.:  3.075  
 Median :35.0   Median : 11.715  
 Mean   :35.0   Mean   : 25.536  
 3rd Qu.:52.5   3rd Qu.: 32.665  
 Max.   :70.0   Max.   :102.300  

plot(degradation$Conc, degradation$Timex)

model <- drm(Conc ~ Time, fct = DRC.expoDecay(), data = degradation)
summary(model)

Model fitted: Exponential Decay Model (2 parms)

Parameter estimates:

                   Estimate Std. Error t-value   p-value    
init:(Intercept) 99.6349312  1.4646680  68.026 < 2.2e-16 ***
k:(Intercept)     0.0670391  0.0019089  35.120 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error:

 2.621386 (22 degrees of freedom)

plot(model, log = "")

This is the graphics plot function. The log argument is a character string which contains “x” if the x axis is to be logarithmic, “y” if the y axis is to be logarithmic and “xy” or “yx” if both axes are to be logarithmic.

Asymptotic regression model

This is a function often used to describe limited growth in ecological models, Monomolecular growth, Mitscherlich law, and Von Bertalanffy law.

The most widespread parameterisation in ecological models is:

\[Y = a - (a-b) \cdot exp(-cX)\]

where a is the maximum attainable Y, b is Y at x = 0, and c is proportional to the relative rate of Y increase while X increases. The first derivative \(dY/dx\) is:

D(expression(a - (a - b) * exp (- c * X)), "X")

(a - b) * (exp(-c * X) * c)

that is:

\[Y' = c \cdot (a-Y)\]

what is interesting here is that we see that the rate of growth is actually max when Y is 0 and it decreases as Y increase.

Let’s examine an example

X <- c(1, 3, 5, 7, 9, 11, 13, 20)
Ye <- asymReg.fun(predictor = X, init = 0.3, m = 5, plateau = 20)
epsilon <- rnorm(8, 0, 0.005)
Y <- Ye + epsilon
model <- drm(Y ~ X, fct = DRC.asymReg())
plot(model,)# log="", main = "Asymptotic regression")

Negative exponential equation

If we adapt the asymptotic equation and add a constraint that \(b=0\), we get the following equation, that is often known as a negative exponential equation.

\[Y = a[1-exp(-cX)]\]

similar to asymptotic but for \(X = 0\), Y is also 0, meaning that the curve passes through the origin. It is often used to model the absorbed photosynthetically active radiation.