Linear Mix Effects Models - A primer
We will use an example from Zuur et al. (2009). Data include the marine benthic data from nine inter-tidal areas along the Dutch coast (Zuur et al. 2007). Five samples were taken at 9 sites and the macro-fauna and abiotic variables were measured, including species richness (the number of different species) and NAP (the height of a sampling station compared to mean tidal level).
setwd("C:/Liem/GEOG515/Spring15/Labs/")
# Read in data
btdata <- read.table("RIKZ.txt", header=TRUE, sep="\t")
btdata## Sample Richness Exposure NAP Beach
## 1 1 11 10 0.045 1
## 2 2 10 10 -1.036 1
## 3 3 13 10 -1.336 1
## 4 4 11 10 0.616 1
## 5 5 10 10 -0.684 1
## 6 6 8 8 1.190 2
## 7 7 9 8 0.820 2
## 8 8 8 8 0.635 2
## 9 9 19 8 0.061 2
## 10 10 17 8 -1.334 2
## 11 11 6 11 -0.976 3
## 12 12 1 11 1.494 3
## 13 13 4 11 -0.201 3
## 14 14 3 11 -0.482 3
## 15 15 3 11 0.167 3
## 16 16 1 11 1.768 4
## 17 17 3 11 -0.030 4
## 18 18 3 11 0.460 4
## 19 19 1 11 1.367 4
## 20 20 4 11 -0.811 4
## 21 21 3 10 1.117 5
## 22 22 22 10 -0.503 5
## 23 23 6 10 0.729 5
## 24 24 0 10 1.627 5
## 25 25 6 10 0.054 5
## 26 26 5 11 -0.578 6
## 27 27 4 11 -0.348 6
## 28 28 1 11 2.222 6
## 29 29 6 11 -0.893 6
## 30 30 4 11 0.766 6
## 31 31 2 11 0.883 7
## 32 32 1 11 1.786 7
## 33 33 1 11 1.375 7
## 34 34 3 11 -0.060 7
## 35 35 4 11 0.367 7
## 36 36 3 10 1.671 8
## 37 37 5 10 -0.375 8
## 38 38 7 10 -1.005 8
## 39 39 5 10 0.170 8
## 40 40 0 10 2.052 8
## 41 41 7 10 -0.356 9
## 42 42 11 10 0.094 9
## 43 43 3 10 -0.002 9
## 44 44 0 10 2.255 9
## 45 45 2 10 0.865 9Assume we want to develop separate regression models for each beach (i.e., site)
Beta <- vector(length = 9)
for (i in 1:9){
Model_i <- summary(lm(Richness ~ NAP,
subset = (Beach==i), data=btdata))
Beta[i] <- Model_i$coefficients[2]
}
Beta## [1] -0.3718279 -4.1752712 -1.7553529 -1.2485766 -8.9001779 -1.3885120
## [7] -1.5176126 -1.8930665 -2.9675304# Instead of the loop in the code above, you can also use the lmList command
# from the nlme package for the same results.Now we put two variables, L.AREA and fGRAZE, into one model
# Create a factor for exposure which has two categories
fExposure9 <- factor(c(0, 0, 1, 1, 0, 1, 1, 0, 0))
# Regression model with slopes of Model_i as depedent variable and Exposure as indepedent variable.
# Actually this is an ANOVA test
Model1 <- lm(Beta ~ fExposure9)
summary(Model1) ##
## Call:
## lm(formula = Beta ~ fExposure9)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.2386 -0.2778 0.0890 0.6940 3.2897
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.662 1.099 -3.332 0.0126 *
## fExposure91 2.184 1.649 1.325 0.2268
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.458 on 7 degrees of freedom
## Multiple R-squared: 0.2005, Adjusted R-squared: 0.08625
## F-statistic: 1.755 on 1 and 7 DF, p-value: 0.2268The Random Intercept Model
library(nlme)
btdata$fBeach <- factor(btdata$Beach)
Mlme1 <- lme(Richness ~ NAP, random = ~1 | fBeach,data = btdata)
summary(Mlme1)## Linear mixed-effects model fit by REML
## Data: btdata
## AIC BIC logLik
## 247.4802 254.525 -119.7401
##
## Random effects:
## Formula: ~1 | fBeach
## (Intercept) Residual
## StdDev: 2.944065 3.05977
##
## Fixed effects: Richness ~ NAP
## Value Std.Error DF t-value p-value
## (Intercept) 6.581893 1.0957618 35 6.006682 0
## NAP -2.568400 0.4947246 35 -5.191574 0
## Correlation:
## (Intr)
## NAP -0.157
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.4227495 -0.4848006 -0.1576462 0.2518966 3.9793918
##
## Number of Observations: 45
## Number of Groups: 9coef(Mlme1)## (Intercept) NAP
## 1 9.203412 -2.5684
## 2 11.781501 -2.5684
## 3 3.966113 -2.5684
## 4 4.306275 -2.5684
## 5 8.532072 -2.5684
## 6 4.952491 -2.5684
## 7 4.816416 -2.5684
## 8 5.520228 -2.5684
## 9 6.158529 -2.5684Model2<- lm(Richness ~ NAP + fBeach, data=btdata)
summary(Model2)##
## Call:
## lm(formula = Richness ~ NAP + fBeach, data = btdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8518 -1.5188 -0.1376 0.7905 11.8384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.8059 1.3895 7.057 3.22e-08 ***
## NAP -2.4928 0.5023 -4.963 1.79e-05 ***
## fBeach2 3.0781 1.9720 1.561 0.12755
## fBeach3 -6.4049 1.9503 -3.284 0.00233 **
## fBeach4 -6.0329 2.0033 -3.011 0.00480 **
## fBeach5 -0.8983 2.0105 -0.447 0.65778
## fBeach6 -5.2231 1.9682 -2.654 0.01189 *
## fBeach7 -5.4367 2.0506 -2.651 0.01196 *
## fBeach8 -4.5530 1.9972 -2.280 0.02883 *
## fBeach9 -3.7820 2.0060 -1.885 0.06770 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.06 on 35 degrees of freedom
## Multiple R-squared: 0.7025, Adjusted R-squared: 0.626
## F-statistic: 9.183 on 9 and 35 DF, p-value: 5.645e-07Now we make plot of Mlme1
F0 <- fitted(Mlme1, level = 0)
F1 <- fitted(Mlme1, level = 1)
I <- order(btdata$NAP)
NAPs <- sort(btdata$NAP)
plot(NAPs, F0[I], lwd = 2, type = "l", ylim = c(0, 22), ylab = "Richness", xlab = "NAP")
for (i in 1:9){
x1 <- btdata$NAP[btdata$Beach == i]
y1 <- F1[btdata$Beach == i]
K <- order(x1)
lines(sort(x1), y1[K])
}
text(btdata$NAP, btdata$Richness, btdata$Beach, cex = 0.5)
The Random Intercept and Slope Model
Mlme2 <- lme(Richness ~ NAP,random = ~ 1 + NAP | fBeach, data = btdata)
summary(Mlme2)## Linear mixed-effects model fit by REML
## Data: btdata
## AIC BIC logLik
## 244.3839 254.9511 -116.1919
##
## Random effects:
## Formula: ~1 + NAP | fBeach
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 3.549064 (Intr)
## NAP 1.714963 -0.99
## Residual 2.702820
##
## Fixed effects: Richness ~ NAP
## Value Std.Error DF t-value p-value
## (Intercept) 6.588706 1.264761 35 5.209448 0e+00
## NAP -2.830028 0.722940 35 -3.914610 4e-04
## Correlation:
## (Intr)
## NAP -0.819
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.8213326 -0.3411043 -0.1674617 0.1921216 3.0397049
##
## Number of Observations: 45
## Number of Groups: 9Now we make plot of Mlme2
F0 <- fitted(Mlme2, level = 0)
F1 <- fitted(Mlme2, level = 1)
I <- order(btdata$NAP); NAPs <- sort(btdata$NAP)
plot(NAPs, F0[I], lwd = 4, type = "l", ylim = c(0, 22), ylab = "Richness", xlab = "NAP")
for (i in 1:9){
x1 <- btdata$NAP[btdata$Beach == i]
y1 <- F1[btdata$Beach == i]
K <- order(x1)
lines(sort(x1), y1[K])
}
text(btdata$NAP, btdata$Richness, btdata$Beach, cex = 0.5)
Random Effects Model
Mlme3 <- lme(Richness ~ 1, random = ~ 1 | fBeach,data = btdata)
summary(Mlme3)## Linear mixed-effects model fit by REML
## Data: btdata
## AIC BIC logLik
## 267.1142 272.4668 -130.5571
##
## Random effects:
## Formula: ~1 | fBeach
## (Intercept) Residual
## StdDev: 3.237112 3.938415
##
## Fixed effects: Richness ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 5.688889 1.228419 36 4.631066 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.77968689 -0.50704111 -0.09795286 0.25468670 3.80631705
##
## Number of Observations: 45
## Number of Groups: 9