At the end of June I posted a question about using nls() with a
factor in the model. As always, I?m amazed and thankful for the
speed and helpfulness of the responses. Not surprsingly I was
redirected down a different path (why documented help could never
replace this list - you have to be able to approximately spell a
word to use a dictionary). Summary and final model follows (only
excerpts are given for brevity).
Problem stated in the original post
As blood in a tube is centrifuged, the blood separates with
serum on top. The experiment measured the percent serum
during centrifuge for 10 donors, with measurements taken
initially and every 10 seconds for 10 minutes (10 donors x 61
measurements per donor).
Graphing the data shows that the curve family
serum ~ bo - b1 * exp( -b2 * time)
is a reasonable model. The asymptote, defined by b0, varies
widely by donor (35 to 55%), but the time till the curve
flattens outs looks consistent (somewhere between 4 and 5
minutes).
The goal is to determine the how long until the curve flattens
out to create a criterion for evaluating new systems. (Say
flattening is defined by the time until the serum percentage
achieves 98% of its final value.)
To answer this question I would like to fit a curve family
where b0 and possibly b1 vary per donor, but b2 is general to
all donors.
So far I can fit separate models for each donor, but have been
unable to make a general model.
There is an example in the S+ 4.5 manual using the Puromycin
data set where the variable 'state' is a factor, but the
factor only has two levels and is converted to 0's and 1's.
Using nls
nls(serum ~ bo[Run] - b1[Run] * exp(-b2 * time), data =
blood, start= list(bo=rep(??, 10), b1=rep(??, 10), b2=??))
where `Run' is a column giving the number of the donor.
suppose you have starting values for b0 and b1 for all
donors. Put them into separate vectors, say b00 and b10.
Let theta0 be the starting value for theta. You can call
nls with a formula as follows:
> fit <- nls(Serum ~ b0[Donor] + b1[Donor]*exp(-theta *
Time), MyData, start = list(b0 = b00, b1 = b11, theta =
theta0), trace = T)
Notice how the starting values must be given as a list.
This solution applies to linear or non-linear parameters.
Using plinear.
If b2 were fixed, then this would be a linear model. That
makes the actual model a partially linear model, and there
are special methods in nls to fit those (the 'plinear'
method).
The main point is that this kind of model is extremely easy
to fit: you can treat it (almost) like a 1 parameter model,
since once you know the value of b2, you can estimate the
other parameters using a fast method like lm(). You could
even draw a graph of the residual sum of squares from these
fits and estimate b2 by eye, by finding the minimum of the
graph.
<snip> Suppose you have a data frame, say `MyData' with
columns Serum, Time and Donor (a factor).
Firstly you could use the `plinear' algorithm in this case
since you only want to let the linear parameters vary with
donor.
Dx <- model.matrix(~Donor - 1, MyData) # creates a binary
matrix
> fit <- nls(Serum ~ cbind(Dx, -Dx*exp(- theta * Time)),
MyData, start = c(theta = 0.04), trace = T)
This will fit different b0 and b1 parameters for each donor
but a common theta. If you wanted different b0 parameters
but the same b1 and theta it would simply be
> fit <- nls(Serum ~ cbind(Dx, -exp(- theta * Time)),
MyData, start = c(theta = 0.04), trace = T)
You need only supply starting values for the non-linear
parameters, in this case theta only.
Using nlme
The general consensus was I should be using nlme as I wasn?t
interested in specific donors and should treat them as a random
effect. The model form also changed as I wanted to restrict the
percent serum to zero when time was zero:
serum ~ b0 (1 - exp( -b1 * time)).
This fit the data well and the time until the serum percentage
achieves 98% of its final value is strictly a function of b1. :
) The code for the model used was:
gelkin.m7 <- nlme(serum ~ b0 * ( 1 - exp(-1 * b1 * time)),
fixed=list(b0~.,b1~.),
random=list(b0~.),
cluster=~donor,
data=gel.kin,
start=list(fixed=c(0.5, 1)))
References
Check out the information on the library and methods at:
http://franz.stat.wisc.edu/pub/NLME
I believe chapter 10 of Venables and Ripley (Modern Applied
Statistics with S-PLUS 2nd ed.) should also prove useful.
There is an analysis of a very similar experiment on OME in
the
on-line complements to Venables & Ripley (Exercises),
including fitting
a common slope for all subjects in nls and using nlme, and
the text has an experiment on blood pressure in Rabbits with
an nls fit using factors for subject on page 318.
Lindstrom, M.J. and Bates, D.M. (1990). Nonlinear Mixed
Effects Models for Repeated Measures Data. Biometrics 46,
673-687.
Thanks again,
KJJ
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
|