[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)

# 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,"treat"]
treati = grep(treat,rownames(lmcoef))
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))