Assignment 11—Solution

Question 1

I read in the data and subset things so only species 2 and trees of age 2 are used.

> water<- read.table('http://www.unc.edu/courses/2006spring/ecol/145/001/data/lab12/WaterUsageData.csv', header=TRUE, sep=',')
> dim(water)
[1] 992 8

Fig. 1   Panel display of water usage for species 2, age 2

> #subset data set
> water2<-water[water$AGE==2 & water$SPECIES==2,]
> dim(water2)
[1] 249 8
> #make tree a factor for better labeling of panels
> tree.f<-factor(water2$TREE)
> library(lattice)
> trellis.par.set(col.whitebg())
> #generate panels as instructed
> xyplot(WU~TIME|tree.f, data=water2,
panel=function(x,y) {
model1<-lm(y~x+I(x^2))
panel.xyplot(x,y)
panel.abline(lm(y~x), lty=1, col=4)
panel.loess(x,y,lty=2)
 panel.curve(coef(model1)[1]+ coef(model1)[2]*x+ coef(model1)[3]*x^2, col=2, lwd=2)
},
strip=strip.custom(par.strip.text = list(cex=0.75)),
key=list(lines=list(lty=c(1,2,1), col=c(4,1,2)), lwd=c(1,1,2),
text=list(c("linear", "lowess",
"quadratic"), cex=rep(.75,3)),
border=0,x=.7,y=.8, corner=c(0,0)))

Question 2

Our basic level-1 model is the following.

I first fit the three models of interest.

> library(nlme)
> model1<-lme(WU~TIME+I(TIME^2),random=~1|TREE,data=water,method='ML')
> model2<-lme(WU~TIME+I(TIME^2),random=~TIME-1|TREE,data=water,method='ML')
> model3<-lme(WU~TIME+I(TIME^2),random=~I(TIME^2)-1|TREE,data=water,method='ML')
> sapply(list(model1,model2,model3),AIC)
[1] 1646.792 1568.143 1650.473

Based on AIC we should choose a model in which TIME is random. This corresponds to the following two-level model.

Question 3

There are three additional models: a model in which the intercept and TIME have random effects, TIME and TIME2 have random effects, and a model in which all three have random effects. I fit all three and compare them to the model selected best in Question 2.

> model2.1<-lme(WU~TIME+I(TIME^2),random=~TIME|TREE,data=water,method='ML')
> model2.2<-lme(WU~TIME+I(TIME^2),random=~TIME+I(TIME^2)-1|TREE,data=water,method='ML')
> model2.3<-lme(WU~TIME+I(TIME^2),random=~TIME+I(TIME^2)|TREE,data=water,method='ML')
> sapply(list(model2,model2.1,model2.2,model2.3),AIC)
[1] 1568.143 1422.660 1424.577 1426.828

The model with random effects for the intercept and TIME ranks best.

Question 4

I fit the continuous autoregressive model below.

> model3<-lme(WU~TIME+I(TIME^2), random=~TIME|TREE,data=water, method='ML', correlation=corCAR1(form=~TIME|TREE))
> sapply(list(model2.1,model3),AIC)
[1] 1422.660 1384.721
> anova(model2.1,model3)
         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
model2.1     1  7 1422.660 1456.958 -704.3299
model3       2  8 1384.721 1423.919 -684.3607 1 vs 2 39.93849  <.0001

Both AIC and the likelihood ratio test argue for including the correlation structure. From the the summary statistics we see that the estimate of phi is 0.75

> summary(model3)
Linear mixed-effects model fit by maximum likelihood
Data: water
     AIC      BIC    logLik
1384.721 1423.919 -684.3607

Random effects:
Formula: ~TIME | TREE
Structure: General positive-definite, Log-Cholesky parametrization
            StdDev       Corr
(Intercept) 1.143061567 (Intr)
TIME        0.006010614 -0.808
Residual    0.443403365

Correlation Structure: Continuous AR(1)
Formula: ~TIME | TREE
Parameter estimate(s):
     Phi
0.750483
Fixed effects: WU ~ TIME + I(TIME^2)
                 Value Std.Error  DF   t-value p-value
(Intercept) -10.948345 0.6286666 950 -17.41519       0
TIME          0.119643 0.0055107 950  21.71092       0
I(TIME^2)    -0.000263 0.0000119 950 -22.07771       0
Correlation:
        ( Intr)   TIME
TIME     -0.980
I(TIME^2) 0.945 -0.982

Standardized Within-Group Residuals:
        Min          Q1        Med         Q3        Max
-4.61619705 -0.51105104 0.01349535 0.55878782 3.50418583

Number of Observations: 992
Number of Groups: 40

 

Course Home Page

Jack Weiss
Phone: (919) 962-5930
E-Mail: jack_weiss@unc.edu
Address: Curriculum in Ecology, Box 3275, University of North Carolina, Chapel Hill, 27516
Copyright © 2006
Last Revised--April 24, 2006
URL: http://www.unc.edu/courses/2006spring/ecol/145/001/docs/solutions/assign11.htm