[R-sig-ME] Interpreting Growth curves
Dieter Menne
dieter.menne at menne-biomed.de
Fri Mar 18 09:12:10 CET 2011
Robert Kushler suggested off-list to use different random structure:
> lme2a = lme(bone~(treat-1):I(time),
+ random = list(animal=pdDiag(~(treat-1))), data=gp,
+ weight = varIdent(form=~1|treat))
....
Parameter estimates:
DrugB DrugA DrugC
1.0000000 0.2192068 0.2753703
Fixed effects: bone ~ (treat - 1):I(time)
Value Std.Error DF t-value p-value
treatDrugA:I(time) 0.0657069 0.01381929 167 4.754724 0.0000
treatDrugB:I(time) 1.5580798 0.21538497 167 7.233930 0.0000
treatDrugC:I(time) 0.0968065 0.03694134 167 2.620546 0.0096
Thanks Robert, that's exactly what I was looking for.
He also questioned the use of a fixed (=0) intercept. I felt not comfortable
with this choice, even if it was sanctioned by definition. However, in my
original formulation the solution was even worse without that assumption. With
his suggestion, using an (Interepts) works fine, and give a less degenerate
estimate for DrugA diagonal structure variance.
David Duffy suggested not to use the growth curve at all, but only the net
outcome after the longest time. Researchers usually use different times after
implant to get some impression of the form of the growth curve; if they were
only interested in the net result, one would follow up all animals as long as
possible.
David also suggested to use a stabilizing sqrt transform. This what I normally
do, but the linear solution showed the problem better. Or use varPower as in the
"final" below.
Thanks to both of you.
Dieter
library(lattice)
library(nlme)
options(digits=2)
gp = read.csv("http://www.menne-biomed.de/uni/bonegrowth.csv",header=TRUE)
lme2b = lme(bone~treat*time,
random = list(animal=pdDiag(~(treat-1))), data=gp,
weight = varPower(fixed=0.5)
)
summary(lme2a)
plot(lme2a)
More information about the R-sig-mixed-models
mailing list