Assignment 10—Solution

Setting up the data

I read in the data.

> inverts<- read.table('http://www.unc.edu/courses/2006spring/ecol/145/001/data/lab11/74species.csv', sep=',', header=TRUE)
> names(inverts)
 [1] "species"      "meantemp" "temp"  "PLD" "taxon"   "climate.3"
 [7] "feeding.type" "lntemp"   "lnPLD" "n"   "PLD.15c" "lnPLD.15c"
[13] "X1000.k"      "lograte"

Question 1

I fit the two models discussed in the question.

> model2<-lme(lnPLD~lntemp, random=~lntemp|species, data=inverts, method="ML")
> model2.5<-lme(lnPLD~I(lntemp-log(15)), random=~I(lntemp-log(15))|species, data=inverts, method="ML")
> summary(model2)$tTable
                Value  Std.Error  DF   t-value      p-value
(Intercept)  7.088667 0.29554220 143  23.98529 5.647789e-52
lntemp      -1.454417 0.08465556 143 -17.18040 1.371520e-36
> summary(model2.5)$tTable
                        Value  Std.Error  DF   t-value      p-value
(Intercept)          3.150033 0.11072433 143  28.44933 9.559019e-61
I(lntemp - log(15)) -1.454417 0.08465556 143 -17.18041 1.371507e-36

The fixed effect portion of model 2 is where x is lntemp. The fixed effect portion of model 2.5 is the following.

Thus if the estimates are truly the same it must be the case that

Plugging in the coefficient estimates from model 2.5 into these formulas we obtain the following for the intercept and slope.

> #intercept
> fixef(model2.5)[1]-fixef(model2.5)[2]*log(15)
(Intercept)
7.088667
> #slope
> fixef(model2.5)[2]
I(lntemp - log(15))
-1.454417

which we see are exactly the values estimated in model 2 that are shown in the summary table above.

Question 2

I first fit all the models of interest.

Model 3a—feeding type affects only the intercept

The multilevel model we're fitting is the following.

Its composite form is shown below.

The composite form gives us the fixed and random expressions required by the lme function of R. To minimize estimation problems due to excessive correlation between the random effects that we've encountered previously, I center the log temperature variable around log(15).

> model3a<-lme(lnPLD~I(lntemp-log(15))+feeding.type, random=~I(lntemp-log(15))|species, data=inverts, method='ML')

Model 3b—feeding type affects only the slope

The multilevel model we're fitting is the following.

which in composite form is the following.

The composite form gives us the fixed and random expressions required by the lme function of R. Once again I center the log temperature variable around log(15).

> model3b<-lme(lnPLD~I(lntemp-log(15)) + I(lntemp-log(15)):feeding.type, random=~I(lntemp-log(15))|species, data=inverts, method='ML')

Model 3c—feeding type affects both the intercept and slope

The multilevel model we're fitting is the following.

which in composite form is the following.

The composite form gives us the fixed and random expressions required by the lme function of R. As explained in Lecture 44 I can use R's shortcut notation to obtain the two main effects and the interaction term. Once again I center the log temperature variable around log(15).

> model3c<-lme(lnPLD~I(lntemp-log(15))*feeding.type, random=~I(lntemp-log(15))|species, data=inverts, method='ML')

Comparing the models

For completeness I also fit a model in which feeding.type is omitted altogether.

> model3<-lme(lnPLD~I(lntemp-log(15)), random=~I(lntemp-log(15))|species, data=inverts, method='ML')

I compare the four models using AIC.

> sapply(list(model3,model3a,model3b,model3c),AIC)
[1] 210.4372 203.2484 211.8966 204.7071

Model 3a, the model in which feeding type affects only the intercept, ranks best. The only model that comes close is one in which feeding type affects both the slope and intercept. But notice that the decrease in the loglikelihood obtained with this model is not even enough to offset the increase in the penalty term.

Since they are nested we can use likelihood ratio tests to test model3a against each of model3 and model3c.

> anova(model3,model3a)
        Model df      AIC      BIC    logLik   Test  L.Ratio p-value
model3      1  6 210.4372 230.7441 -99.21859
model3a     2  7 203.2484 226.9399 -94.62421 1 vs 2 9.188761  0.0024
> anova(model3a,model3c)
        Model df      AIC      BIC    logLik   Test  L.Ratio p-value
model3a     1  7 203.2484 226.9399 -94.62421
model3c     2  8 204.7071 231.7830 -94.35353 1 vs 2 0.541343  0.4619

From the first test we see that adding feeding type to the intercept equation is a significant improvement over a model without feeding type. On the other hand, the second test tells us that there is no evidence that the magnitude of the relationship between larval development time and temperature differs for larva of different feeding types. Therefore feeding type does not modify the slope.

Question 3

The two-stage and composite versions of model 3a were shown in Question 2 and are repeated below.

Part 1

Two stage model

Part 2

Composite model

Part 3

Population-averaged (marginal) model

Subject-specific (conditional) model

Part 4.  R has set up the contrasts so that feeding.type = 0 for lecithotrophic larvae and is equal to 1 for planktotrophic larva.

> contrasts(inverts$feeding.type)
  P
L 0
P 1

Marginal models

Conditional models

Question 4

To obtain a pseudo-R² we compare the variance component of the random slopes and intercepts model (model3) with its corresponding value in model 3a.

> VarCorr(model3)
species = pdLogChol(I(lntemp - log(15)))
                    Variance   StdDev    Corr
(Intercept)         0.85672246 0.9255930 (Intr)
I(lntemp - log(15)) 0.28109623 0.5301851 -0.49
Residual            0.02310587 0.1520062
> VarCorr(model3a)
species = pdLogChol(I(lntemp - log(15)))
                    Variance   StdDev    Corr
(Intercept)         0.75255594 0.8674998 (Intr)
I(lntemp - log(15)) 0.28135382 0.5304280 -0.479
Residual            0.02303037 0.1517576

Observe that except for , labeled (Intercept) in the output, the remaining components have changed hardly at all. This is good because it means the pseudo-R² can be interpreted more like a real R². The formula is

> (as.numeric(VarCorr(model3)[1,1]) - as.numeric(VarCorr(model3a)[1,1])) / as.numeric(VarCorr(model3)[1,1])
[1] 0.1215872

So roughly 12% of the variability in intercepts between species is explained by their feeding type.

Question 5

From the marginal and conditional equations obtained in Question 3, we see that the two feeding types differ only in their intercept. From a printout of the fixed effects we observe the estimate for β₂ is positive.

> fixef(model3a)
(Intercept) I(lntemp - log(15)) feeding.typeP
  2.6374682          -1.4489413     0.6773972

Thus at any fixed temperature planktotrophic larvae have a larval duration time that is longer than that of lecithotrophic larvae. A plot of the population-averaged models for each feeding type is shown below superimposed on the raw data. I use the curve function to simplify the plotting. With the curve function it is not necessary to specify a range for the x-variable; it is determined automatically from the plot function when you use the add=TRUE option of curve.

> plot(inverts$lntemp, inverts$lnPLD, pch= (as.numeric(inverts$feeding.type) -1)*18 + 2, xlab='log(temperature)', ylab='log(PLD)')
> curve(fixef(model3a)[1] + fixef(model3a)[2] * (x-log(15)), col=2, add=TRUE)
> curve(fixef(model3a)[1] + fixef(model3a)[3] + fixef(model3a)[2] * (x-log(15)), col=4, add=TRUE)
> legend(1,.8, c('lecithotrophic', 'planktotrophic'), col=c(2,4), lty=c(1,1), cex=c(.9,.9), bty='n')

Fig. 1 is a bit dishonest in that it suggests we are fitting lines to the raw data. That's not correct. Instead the displayed lines are the population averaged lines for two families of subject-specific lines. A more honest picture is one in which the subject-specific lines are also displayed. The next bit of code carries this out.

> plot(inverts$lntemp, inverts$lnPLD, pch=(as.numeric(inverts$feeding.type)-1)*18+2, xlab='log(temperature)', ylab='log(PLD)')
#get list of species names for each feeding type
> unique(inverts[inverts$feeding.type=='L', "species"])->lecitho
> unique(inverts[inverts$feeding.type=='P', "species"])->plankto
> random.effs<-ranef(model3a)
#separate random effects by feeding type
> lecitho.ran<-random.effs[rownames(random.effs)%in%lecitho,]
> plankto.ran<-random.effs[rownames(random.effs)%in%plankto,]
#plot subject specific curves--planktotrophic
> apply(cbind(plankto.ran[,1]+fixef(model3a)[1]+fixef(model3a)[3],
plankto.ran[,2]+fixef(model3a)[2]),1,
function(y) curve(y[1]+y[2]*(x-log(15)), col='skyblue', add=TRUE, lty=2))
#plot subject specific curves--lecithotrophic
> apply(cbind(lecitho.ran[,1]+fixef(model3a)[1],
lecitho.ran[,2]+fixef(model3a)[2]),1,
function(y) curve(y[1]+y[2]*(x-log(15)), col='pink', add=TRUE))
#plot population-averaged curves
> curve(fixef(model3a)[1] + fixef(model3a)[2] * (x-log(15)), col=2, add=TRUE, lwd=2)
> curve(fixef(model3a)[1] + fixef(model3a)[3] + fixef(model3a)[2] * (x-log(15)), col=4, lty=2, add=TRUE, lwd=2)
> legend(1,.8, c('lecithotrophic', 'planktotrophic'), col=c(2,4), lty=c(1,2), cex=c(.9,.9), bty='n')

Fig. 2  Conditional and marginal models for feeding types

Observe that there is one clearly deviant lecithotrophic species (top right corner) and perhaps three or four truly deviant planktotrophic species (these have subject-specific lines that are far below the lecithotrophic population-averaged line)

Course Home Page