[R-sig-ME] Interpreting Growth curves

[Dieter Menne]

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