Lab3: Nonlinear Models and Least Squares


Getting Started

To Start R: Click on the R icon

If no R icon, go to Start and then All Programs to find R


Example1: Michaelis Menton Relationship

R function: nls

> help(nls)

What is the default algorithm?

Here is an example:

# fitting Michaelis and Menten's original data
conc <- c(0.3330, 0.1670, 0.0833, 0.0416, 0.0208, 0.0104, 0.0052)
vel <- c(3.636, 3.636, 3.236, 2.666, 2.114, 1.466, 0.866)
Micmen <- data.frame(conc=conc, vel=vel)

Let's plot the data.

> plot(conc,vel, pch="+")

You must give starting values for the optimization.

fit <- nls(vel~Vm*conc/(K+conc),Micmen,start=list(Vm=3.7,K=0.02))

What is Vm? What is K?

What's in an object of class "nls" ?
There is a summary and a plot function! What other diagnostic plots would be useful?

There are also functions resid, fitted, predict, and AIC function.


Example2: Michaelis Menton Relationship

Data from a biochemical experiment where the intial velocity of a reaction measured for different concentations of the substrate are given in the data frame Puromycin . For more information use

> help(Puromycin)

Plot velocity against concentration separately for the two levels of treatment.

> plot(Puromycin$conc, Puromycin$rate, type="n")
> text(Puromycin$conc, Puromycin$rate, ifelse(Puromycin$state== "treated", "T", "U"))

Fit nonlinear model to treated data.

> Treated <- Puromycin[Puromycin$state == "treated",]
> Purfit1 <- nls(rate ~ Vm*conc/(K + conc), Treated, list(Vm = 200, K= 0.1), trace=T)

What does the trace=T option do?

Use the summary command and make some diagonistic plots for the model fit.

What is nls.control? To see, use the help function.


Example3: Logistic Growth of Chickens

Weight versus age of chicks on different diets For more information use

> help(ChickWeight)

What does the coplot function do? Copy and paste from the (ChickWeight) help file to see.

Plot chick weight over time by chicks.

> plot(ChickWeight$Time, ChickWeight$weight, type="n")
> text(ChickWeight$Time, ChickWeight$weight, ChickWeight$Chick)

Let's pick out information for Chick #1, a representative Chick

> Chick.1 <- ChickWeight[ChickWeight$Chick == 1, ]
> lines(Chick.1$Time, Chick.1$weight)

Let's fit a logistic growth model:

> fit1 <- nls(weight~Asym/(1+exp((xmid-Time)/scal)),Chick.1,start=list(Asym=368,xmid=14,scal=6),trace=T)

What happens if you change scal=6 to scal=0.1?
What happens if you change xmid=14 to xmid=0.1?

How did I come up with starting values? There is a getInitial and SSlogis function in R.

> getInitial(weight ~ SSlogis(Time, Asym, xmid, scal), data = Chick.1)

Note: Initial values are in fact the converged values.


How would you make a function?

Make your own function myplotnls that will make appropriate summary plots for an nls object.

Here is a place to start: myplotnls1.r

Here is mine: myplotnls.r (Feel free to use this one!)