Lecture 44 (lab 11b)—April 5 , 2006

What was covered?

• R topics covered
  • fixef for extracting fixed effects
• Entering level-2 predictors into a mixed effects model
• Calculating pseudo-R² statistics
• Centering to achieve numerical stability

R functions and commands demonstrated

• fixef (from nlme) extracts fixed effect parameter estimates from an lme model object

R packages used

• nlme contains the functions lme, VarCorr, and fixef for mixed effect models

Summary of the last session

• We continue with exploring the larval development time data set that we studied yesterday.

> #obtain 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"

• I fit the models from last time and then compare them using AIC. For this it's important that the mixed effect models are fit using method="ML".

> library(nlme)
> modelOLS<-lm(lnPLD~lntemp, data=inverts)
> model0<-lme(lnPLD~1,random=~1|species, data=inverts, method='ML')
> model1<-lme(lnPLD~lntemp, random=~1|species, data=inverts, method='ML')
> model2<-lme(lnPLD~lntemp, random=~lntemp|species, data=inverts, method='ML')
> mod.names<-c('OLS', 'uncond. mean', 'random ints', 'random ints/slopes')
> mod.aic<-sapply(list(modelOLS, model0, model1, model2), AIC)
> data.frame(mod.names, mod.aic)
           mod.names  mod.aic
1                OLS 553.0685
2       uncond. mean 451.2372
3        random ints 239.5597
4 random ints/slopes 210.4372

• As noted previously any model that accounts for the structure of the data set is a huge improvement over the unstructured model. It's worth noting that we've accounted for the structure by estimating only three additional parameters (in the random slopes and intercepts model), the two variances and and one covariance that define the distribution of the random effects.
• Graphs of the four models we've fit thus far are shown below. In the figures,
  • and refer to the variability in the intercepts of the subject-specific regression lines about the intercept of population-averaged line,
  • is the variability of the individual slopes about the population line slope, and
  • is the average variability of a species values about its own individual regression line.

The remaining parameter is not shown.

Fig. 1 The OLS model and the three basic multilevel models for the larval development data set.

Adding level-2 predictors

• The climate.3 variable in the data set is a level-2 predictor. It varies only at the species level and is a categorical variable with three levels.

> table(inverts$climate.3)

polar temperate tropical
   29       152       37

• This last display is not overly informative since it is counting the multiple observations from the same species. To get a count of how many species there are from each geographic region do the following.

> table(tapply(inverts$climate.3, inverts$species, function(x) x[1]))

1  2  3
8 53 13

• In the tapply function call we group the climate observations by species but the generic function then only takes the first one. The result is that we get one observation per species and hence a correct count of the totals when the table function is applied.
• Observe that the data are very unbalanced. The preponderance of species in the data set are classified as temperate.
• The climate.3 variable is a factor variable in R. By default dummy variables have been constructed that contrast temperate and tropical species with polar species.

> contrasts(inverts$climate.3)
          temperate tropical
polar             0        0
temperate         1        0
tropical          0        1

• To include climate.3 in the model we have to decide in what way we think it alters the basic level-1 relationship between lnPLD and lntemp. Does it modify the intercept, the slope, or both? In general theory should guide us in answering this question. Since we as outsiders don't have theory to guide us, I'll try all three possibilities.

Does climate.3 affect the population intercept?

• If climate.3 affects the intercept it should be entered into the level-2 equation for the intercept.

• The corresponding composite model is shown below.

• To estimate this model in R we start with the random slopes and intercepts formulation and just add climate to the fixed part of the model.

> model3<-lme(lnPLD~lntemp+climate.3, random=~lntemp|species, data=inverts, method='ML')
Error in solve.default(pdMatrix(a, fact = TRUE)) :
        system is computationally singular: reciprocal condition number = 1.26964e-053

• On a PC, although not on a Macintosh, running the model yields the error message shown above. The problem arises because of the high degree of correlation between the random effects. We can reduce this correlation and get convergence of the estimation routine if we center the predictor lntemp. I begin by centering lntemp about its mean.

> model3<-lme(lnPLD~I(lntemp-mean(lntemp))+climate.3, random=~I(lntemp-mean(lntemp))|species, data=inverts, method='ML')

The model now converges.

• Using the mean for this data set is not particularly appealing. Since the data are not a random sample of any sort (except perhaps of the literature), the mean isn't likely to have any natural interpretation. In particular it does not represent the average temperature at which these marine species are likely to be found. Given this I elect to use a whole number that is near the mean to center the data.

> mean(inverts$lntemp)
[1] 2.809248
> exp(mean(inverts$lntemp))
[1] 16.59743

• The mean of lntemp is equal to the logarithm of the geometric mean of temperature. Since the geometric mean is 16.59, I instead center the data about a temperature of 15 degrees Celsius. This works as the model once again converges.

> model3<-lme(lnPLD~I(lntemp-log(15))+climate.3,random=~I(lntemp-log(15))|species,data=inverts, method='ML')
> summary(model3)
Linear mixed-effects model fit by maximum likelihood
 Data: inverts
       AIC      BIC    logLik
  213.8831 240.9590 -98.94154

Random effects:
Formula: ~I(lntemp - log(15)) | species
Structure: General positive-definite, Log-Cholesky parametrization
                    StdDev    Corr
(Intercept)         0.9099084 (Intr)
I(lntemp - log(15)) 0.5313982 -0.456
Residual            0.1519110

Fixed effects: lnPLD ~ I(lntemp - log(15)) + climate.3
                         Value Std.Error  DF    t-value p-value
(Intercept)          2.9544791 0.3136574 143   9.419446  0.0000
I(lntemp - log(15)) -1.4529806 0.0852620 143 -17.041368  0.0000
climate.3temperate   0.1940039 0.3313066  71   0.585572  0.5600
climate.3tropical    0.3081079 0.3860982  71   0.798004  0.4275
Correlation:
                   (Intr)   I(-l(1 clmt.3tm
I(lntemp - log(15)) -0.160
climate.3temperate  -0.925  0.019
climate.3tropical   -0.791 -0.004 0.749

Standardized Within-Group Residuals:
        Min          Q1         Med         Q3        Max
-2.50718336 -0.36428146 -0.01874458 0.40493293 3.66881543

Number of Observations: 218
Number of Groups: 74

• Because climate.3 is a categorical variable with more than two categories, the Wald tests shown in the summary output do not provide a test of the construct climate. For that we need to carry out a likelihood ratio test.

> anova(model2,model3)
       Model df      AIC      BIC    logLik   Test  L.Ratio p-value
model2     1  6 210.4372 230.7441 -99.21859
model3     2  8 213.8831 240.9590 -98.94154 1 vs 2 0.554096   0.758

• We can further quantify the relationship between climate and the intercept of the population-averaged model by computing a pseudo-R² statistic. In order to do this I first need to estimate a centered version of the random slopes and intercepts model.

> model2.5<-lme(lnPLD~I(lntemp-log(15)), random=~I(lntemp-log(15))|species, data=inverts, method='ML')
> VarCorr(model2.5)
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(model3)
species = pdLogChol(I(lntemp - log(15)))
                    Variance   StdDev    Corr
(Intercept)         0.82793327 0.9099084 (Intr)
I(lntemp - log(15)) 0.28238407 0.5313982 -0.456
Residual            0.02307696 0.1519110

• Since climate.3 was entered in the level-2 intercept equation we should look to the intercept variance component, , for its effect. From the output we can see the effect has been minimal. Formally we compute the statistic as follows.

> (as.numeric(VarCorr(model2.5)[1,1]) - as.numeric(VarCorr(model3)[1,1])) / as.numeric(VarCorr(model2.5)[1,1])
[1] 0.03360387

So only 3% of the species-level variation in the intercepts is explained by climate. This is consistent with the fact that adding climate.3 to the level-2 equation for the intercept did not significantly improve the model.

• Note: centering has no effect on the likelihood ratio. We get exactly the same result as was obtained above using the uncentered version of the unconditional growth model.

> anova(model2.5,model3)
         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
model2.5     1  6 210.4372 230.7441 -99.21859
model3       2  8 213.8831 240.9590 -98.94154 1 vs 2 0.554096   0.758

Does climate.3 affect the population slope?

• To see if climate.3 affects the slope we enter it into the level-2 equation for the slope.

or in composite form

• In R we specify this model by entering an interaction term for the interaction between climate.3 and lntemp. Thus we need to include the term I(lntemp-log(15)):climate.3 in the fixed part of the model expression.

> model4<-lme(lnPLD~I(lntemp-log(15)) + I(lntemp-log(15)):climate.3, random=~I(lntemp-log(15))|species, data=inverts, method='ML')
> anova(model4,model2)
       Model df      AIC      BIC    logLik   Test  L.Ratio p-value
model4     1  8 206.7896 233.8655 -95.39479
model2     2  6 210.4372 230.7441 -99.21859 1 vs 2 7.647598  0.0218

• The likelihood ratio test finds the effect of climate on the slope to be statistically significant. This is confirmed by looking at the AIC. We can quantify the effect with a pseudo-R² statistic. This time we compare the change in between the current model and random slopes and intercepts model.

> VarCorr(model4)
species = pdLogChol(I(lntemp - log(15)))
                    Variance   StdDev    Corr
(Intercept)         0.86149947 0.9281700 (Intr)
I(lntemp - log(15)) 0.20631864 0.4542231 -0.453
Residual            0.02314854 0.1521464
> VarCorr(model2.5)
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(model2.5)[2,1]
[1] "0.28109623"
> (as.numeric(VarCorr(model2.5)[2,1]) - as.numeric(VarCorr(model4)[2,1])) / as.numeric(VarCorr(model2.5)[2,1])
[1] 0.2660213

So 26.6% of the variability in slopes between species is explained by climate.3.

Does climate.3 affect the population slopes and intercepts?

• If climate.3 affects both the slope and intercept it should be entered into both of the level-2 equations.

or in composite form

• From the composite form the fixed part of the model should include lntemp, climate.3, and the interaction between the two. We can obtain this in R as follows.

fixed=lnPLD~I(lntemp-log(15))+climate.3+I(lntemp-log(15)):climate.3

or using R's shortcut notation

fixed=lnPLD~I(lntemp-log(15))*climate.3

where the * notation automatically generates all three terms (plus the intercept).

> model5<-lme(lnPLD~I(lntemp-log(15))*climate.3, random=~I(lntemp-log(15))|species, data=inverts, method='ML')
> anova(model5,model2)
       Model df      AIC      BIC    logLik   Test  L.Ratio p-value
model5     1 10 206.4266 240.2715 -93.21330
model2     2  6 210.4372 230.7441 -99.21859 1 vs 2 12.01058  0.0173

• So compared to a model without climate, we find evidence for a significant climate effect on the slopes and intercepts. We can obtain a more focused test by comparing the current model with our previous one, a model in which climate is only allowed to affect the slopes.

> anova(model5,model4)
       Model df      AIC      BIC    logLik   Test  L.Ratio p-value
model5     1 10 206.4266 240.2715 -93.21330
model4     2  8 206.7896 233.8655 -95.39479 1 vs 2 4.362981  0.1129

• So it appears we should allow the slopes to be affected by climate, but not the intercepts. AIC more or less supports this proposition. The AIC for the last two models we've fit are barely distinguishable although they are both quite a bit lower than that of a model in which climate is absent and a model in which climate only affects the intercepts.

> sapply(list(model2,model3,model4,model5),AIC)
[1] 210.4372 213.8831 206.7896 206.4266

• To examine the estimates for the fixed part of the model we can use the summary function, or we can apply the function fixef to the model.

> fixef(model4)
                           (Intercept)                   I(lntemp - log(15))
                             3.1724449                            -1.0277108
I(lntemp - log(15)):climate.3temperate I(lntemp - log(15)):climate.3tropical
                            -0.4736086                            -0.7131552

Observe that as we move from polar to temperate to tropical the slope becomes more negative.

• We can plot these results in the usual way by constructing functions that generate the appropriate slope terms based on the contrasts that are coded into the dummy variables. Since climate consists of two dummy variables, our slope model actually is the following.

where and are dummy variables created by R. Using the contrasts defined for these variables we can write the slope equations for the three regions.

• Based on this the following three functions create the population-averaged regression lines for the three regions.

> polar.func<-function(x) fixef(model4)[1] + fixef(model4)[2]*(x-log(15))
> temp.func<-function(x) fixef(model4)[1] + (fixef(model4)[2] + fixef(model4)[3])*(x-log(15))
> trop.func<-function(x) fixef(model4)[1] + (fixef(model4)[2] + fixef(model4)[4])*(x-log(15))

• I plot the three regression lines in Fig. 2

> plot(inverts$lntemp,inverts$lnPLD, type='n', xlab='log(temperature) ', ylab='log(PLD) ')
> lines(log(seq(2,42,10)), polar.func(log(seq(2,42,10))), col=1)
> lines(log(seq(2,42,10)), temp.func(log(seq(2,42,10))), col=2)
> lines(log(seq(2,42,10)), trop.func(log(seq(2,42,10))), col=3)
> legend(1,1.5,c('polar','temperate','tropic'), lty=c(1,1,1), col=c(1,2,3), bty='n', cex=c(.8,.8,.8))
> mtext('Marginal model: climate affects the slope', side=3, line=.5)

Fig. 2  Effect of climate on population-averaged model

• Observe the common intersection point of the lines. This occurs at log(15) and follows from the fact that in the model we're using only the slopes were assumed to be affected by climate.

Course Home Page