Statistics for Laboratory Scientists ( 140.615 )
Non-linear Regression
Example from class: Michaelis-Menten equations
The data.
library(SPH.140.615)
puro## conc vel state
## 1 0.02 76 treated
## 2 0.02 47 treated
## 3 0.06 97 treated
## 4 0.06 107 treated
## 5 0.11 123 treated
## 6 0.11 139 treated
## 7 0.22 159 treated
## 8 0.22 152 treated
## 9 0.56 191 treated
## 10 0.56 201 treated
## 11 1.10 207 treated
## 12 1.10 200 treated
## 13 0.02 67 untreated
## 14 0.02 51 untreated
## 15 0.06 84 untreated
## 16 0.06 86 untreated
## 17 0.11 98 untreated
## 18 0.11 115 untreated
## 19 0.22 131 untreated
## 20 0.22 124 untreated
## 21 0.56 144 untreated
## 22 0.56 158 untreated
## 23 1.10 160 untreatedPlot the data, just for the untreated cells.
plot(vel ~ conc, data=puro, subset=(state=="untreated"), xlab="concentration", ylab="initial velocity", lwd=2)
The non-linear least squres fit.
nls.out <- nls(vel ~ (Vm*conc)/(K+conc), data=puro, subset=(state=="untreated"), start = c(Vm=143,K=0.031))
summary(nls.out)$param## Estimate Std. Error t value Pr(>|t|)
## Vm 160.28001082 6.480242918 24.733642 1.384612e-09
## K 0.04770813 0.007781868 6.130678 1.727055e-04Plot the data with the fitted curve.
plot(vel ~ conc, data=puro, subset=(state=="untreated"), xlab="concentration", ylab="initial velocity", lwd=2)
u <- par("usr")
x <- seq(u[1],u[2],length=250)
y <- predict(nls.out, data.frame(conc=x))
lines(x, y, lwd=2)
Non-linear regression estimation for just the treated cells.
nls.outB <- nls(vel ~ (Vm*conc)/(K+conc), data=puro, subset=(puro$state=="treated"), start=c(Vm=160,K=0.048))
summary(nls.outB)$param## Estimate Std. Error t value Pr(>|t|)
## Vm 212.68377326 6.947156912 30.614505 3.241167e-11
## K 0.06412133 0.008280955 7.743229 1.565132e-05Fit both regressions all in one.
puro$x <- ifelse(puro$state=="treated",1,0)
nls.outC <- nls(vel ~ ((Vm+dV*x)*conc)/(K+dK*x+conc), data=puro, start=c(Vm=160,K=0.048,dV=0,dK=0))
summary(nls.outC)$param## Estimate Std. Error t value Pr(>|t|)
## Vm 160.28012885 6.896020314 23.242410 2.040954e-15
## K 0.04770831 0.008281173 5.761057 1.496419e-05
## dV 52.40364391 9.551021939 5.486705 2.712599e-05
## dK 0.01641302 0.011428986 1.436087 1.672387e-01Maximum velocity for the untreated cells.
summary(nls.outC)$param[1]## [1] 160.2801Maximum velocity for the treated cells.
summary(nls.outC)$param[1]+summary(nls.outC)$param[3]## [1] 212.6838Rate constant for the untreated cells.
summary(nls.outC)$param[2]## [1] 0.04770831Rate constant for the treated cells.
summary(nls.outC)$param[2]+summary(nls.outC)$param[4]## [1] 0.06412133Plot data and the fits.
plot(vel ~ conc, data=puro, xlab="concentration", ylab="initial velocity", lwd=2,
pch=as.numeric(puro$state)-1, col=ifelse(puro$state=="treated","blue","red"))
u <- par("usr")
x <- seq(u[1],u[2],len=250)
y1 <- predict(nls.out, data.frame(conc=x))
y2 <- predict(nls.outB, data.frame(conc=x))
lines(x, y1, lwd=2, col="red", lty=2)
lines(x, y2, lwd=2, col="blue", lty=2)
Fit with the K’s constrained to be equal.
nls.outD <- nls(vel ~ ((Vm+dV*x)*conc)/(K+conc), data=puro, start=c(Vm=160,K=0.048,dV=0))
summary(nls.outD)$param## Estimate Std. Error t value Pr(>|t|)
## Vm 166.60407666 5.807428495 28.688098 1.006288e-17
## K 0.05797179 0.005910172 9.808816 4.372679e-09
## dV 42.02596377 6.272141101 6.700417 1.605938e-06Test if the lines are identical.
nls.outE <- nls(vel ~ (Vm*conc)/(K+conc), data=puro, start=c(Vm=160,K=0.048))
summary(nls.outE)$param## Estimate Std. Error t value Pr(>|t|)
## Vm 190.80613580 8.7645774 21.770147 6.835199e-16
## K 0.06038853 0.0107685 5.607886 1.448696e-05anova(nls.outE,nls.outC)## Analysis of Variance Table
##
## Model 1: vel ~ (Vm * conc)/(K + conc)
## Model 2: vel ~ ((Vm + dV * x) * conc)/(K + dK * x + conc)
## Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
## 1 21 7276.5
## 2 19 2055.1 2 5221.5 24.138 6.075e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Example from class: the mud snails
The data.
sediment## time pellets
## 1 6 119
## 2 6 111
## 3 6 87
## 4 12 348
## 5 12 210
## 6 12 272
## 7 18 438
## 8 18 778
## 9 18 548
## 10 24 774
## 11 24 896
## 12 24 670
## 13 42 896
## 14 42 729
## 15 42 886
## 16 78 763
## 17 78 954
## 18 78 954Fit a “linear spline” model.
f <- function(x,b0,b1,x0) ifelse(x<x0, b0+b1*x, b0+b1*x0)
nls.out <- nls(pellets ~ f(time,b0,b1,x0), data=sediment, start=c(b0=0,b1=50,x0=30))
summary(nls.out)$param## Estimate Std. Error t value Pr(>|t|)
## b0 -146.00000 71.292691 -2.047896 5.849712e-02
## b1 38.90556 4.338735 8.967027 2.055707e-07
## x0 25.95173 1.780619 14.574561 2.911998e-10Plot the data and the fitted curve.
plot(pellets ~ time, data=sediment)
u <- par("usr")
x <- seq(u[1],u[2],length=250)
y <- predict(nls.out, data.frame(time=x))
lines(x, y, col="red", lwd=2)
Example from class: exponential decay
The data.
head(chlorine)## time chlorine
## 1 8 0.49
## 2 8 0.49
## 3 10 0.48
## 4 10 0.47
## 5 10 0.48
## 6 10 0.47Fit the model.
nls.out <- nls(chlorine ~ a+b*exp(-c*time), data=chlorine, start=c(a=0.05,b=0.49,c=0.1))
summary(nls.out)$param## Estimate Std. Error t value Pr(>|t|)
## a 0.38962769 0.005841415 66.700906 1.931064e-43
## b 0.21939846 0.031370595 6.993762 1.681635e-08
## c 0.09915567 0.018103808 5.477062 2.391030e-06Plot the data and the fitted curve.
plot(chlorine ~ time, data=chlorine, lwd=2)
u <- par("usr")
x <- seq(u[1],u[2],len=250)
y <- predict(nls.out, data.frame(time=x))
lines(x, y, lwd=2, col="blue")