[R-sig-ME] Interpreting Growth curves
Dieter Menne
dieter.menne at menne-biomed.de
Thu Mar 17 17:15:54 CET 2011
Friends of nlme and lmer,
I have growth curves of new bone material under 3 treatments a, b, c. By
definition, at t=0 we have no new bone, so fitting a linear slope without
intercept is a reasonable model. For each animal, data are available only at
one point in time, but several samples are taken at each location with large
variation.
Treatment a and c are within animal (different locations in body).
Treatment b is a much more effective that could not be paired with a and c;
it has a much higher variance.
I used varIdent to handle heteroscedasticity. Things are fine when only the
paired data are analyzed:
DrugA DrugC
1.0 1.3
Fixed effects: bone ~ (treat - 1):I(time)
Value Std.Error DF t-value p-value
treatDrugA:I(time) 0.082 0.031 100 2.7 0.0088
treatDrugC:I(time) 0.097 0.032 100 3.0 0.0031
....
Standard error of 30% are reasonable. However, when I add DrugB, I get
standard errors of 100% for A and C, at similar estimates for the slopes.
DrugB DrugA DrugC
1.00 0.22 0.28
Fixed effects: bone ~ (treat - 1):I(time)
Value Std.Error DF t-value p-value
treatDrugA:I(time) 0.09 0.11 167 0.8 0.43
treatDrugB:I(time) 1.62 0.15 167 10.9 0.00
treatDrugC:I(time) 0.10 0.11 167 0.9 0.38
The way out would be to analyze the data separately, but I don't like the
approach. Any alternatives?
Full code below, data can be loaded from the net.
Dieter
# -------------------------------------------------------------------------
library(lattice)
library(nlme)
options(digits=2)
gp = read.csv("http://www.menne-biomed.de/uni/bonegrowth.csv",header=TRUE)
# By definition, bone = 0 for time = 0
# Without DrugB (2 treatments within on animal)
# Using varIdent to model heteroscedasticity
lme1 = lme(bone~(treat-1):I(time),
subset=treat != "DrugB",random = ~1|animal, data=gp,
weight = varIdent(form=~1|treat))
summary(lme1)
# With Drug B (parallel group, different animals)
lme2 = lme(bone~(treat-1):I(time),
random = ~1|animal, data=gp,
weight = varIdent(form=~1|treat))
summary(lme2)
doplot = function(lmeRes,main){
# Plot jittered data an fits
lmcoef = summary(lmeRes)$tTable
xyplot(bone~jitter|treat,data=gp,
panel = function(x,y,subscripts,groups,...){
panel.superpose(x,y,subscripts,groups,...)
treat = gp[subscripts[1],"treat"]
treati = grep(treat,rownames(lmcoef))[1]
xp = seq(0,max(x),length.out=50)
yp = xp*lmcoef[treati,1]
panel.xyplot(xp,yp,type="l", lwd=3,col="black")
yp = xp*(lmcoef[treati,1]+lmcoef[treati,2])
panel.xyplot(xp,yp,type="l", lwd=3,col="gray")
yp = xp*(lmcoef[treati,1]-lmcoef[treati,2])
panel.xyplot(xp,yp,type="l", lwd=3,col="gray")
},
layout = c(1,3),
main = main,
groups=animal,cex=0.6,pch=16,xlim=c(0,18),
scales=list(x=list(alternating=1,at=unique(gp$Time))))
}
print(doplot(lme1, "Drug B not fitted"),split=c(1,1,2,1),more=TRUE)
print(doplot(lme2, "Drug B fitted"),split=c(2,1,2,1))
More information about the R-sig-mixed-models
mailing list