Ecological mixed model example
Two more things: (1) I accidentally left out the main effects for "year" in the fixed-effect part of the lmer model. So the call should actually be:
mod.2 = lmer(y ~ year + ys + sm + (1 | site) + (0 + ys | site) )
or for direct estimates of (almost all) cell means
mod.3 = lmer(y ~ 0 + ys + sm + (1 | site) + (0 + ys | site) )
In real data, year should probably be specified with trends, e.g., as poly(year, 2). (2) I think , in principle, replacing
mod = lmer(y ~ year/season + season/month + (1 + year/season | site) )
with
mod = lmer(y ~ year/season + season/month + (1 | site) + (0 + year/season | site) )
should also work for the sigmas in your script. I did run into a convergence problem with this specification. Also, the fixed effects were not very transparent to me in this specification. Best Reinhold Kliegl On Thu, May 14, 2009 at 11:38 PM, Reinhold Kliegl
<reinhold.kliegl at gmail.com> wrote:
The following parameterization recovers both the fixed effects and
also sigma_mu, sigma_beta, and sigma. The correlation does not look
bad either. There are probably also ways of forcing some of the
irrelevant correlations to zero.
Reinhold Kliegl
############################################
....
P = rmvnorm(nsite,mean=theta,sigma=V, method="chol")
...
# specifying seasons nested within year
ys ?<- factor(paste(year, season, sep="_"))
cmat.ys <- matrix(c( ? ? ? ? ?-1, ?1, rep(0,8),
? ? ? ? ? ? ? ? ? ? ? ? ? ?rep(0,2), -1, ?1, rep(0,6),
? ? ? ? ? ? ? ? ? ? ? ? ? ?rep(0,4), -1, ?1, rep(0,4),
? ? ? ? ? ? ? ? ? ? ? ? ? ?rep(0,6), -1, ?1, rep(0,2),
? ? ? ? ? ? ? ? ? ? ? ? ? ?rep(0,8), -1, ?1 ? ? ? ? ? ? ), ?10, ?5)
colnames(cmat.ys) <- c(".s2-s1|y1", ".s2-s1|y2", ".s2-s1|y3",
".s2-s1|y4",".s2-s1|y5")
contrasts(ys, 5) <- cmat.ys
# specifying months nested within season
# contrasting months 2 vs 1 and 3 vs 2 within season
sm ?<- factor(paste(season, month, sep="_"))
cmat.sm <- matrix(c(-2/3, ?1/3, ?1/3, ? ?0, ? ? 0, ? ? 0,
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?-1/3, -1/3, ?2/3, ? ?0, ? ? 0, ? ? 0,
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?0, ? ? 0, ? ? 0, -2/3, ?1/3, ?1/3,
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?0, ? ? 0, ? ? 0, -1/3, -1/3, ?2/3 ?), ?6, ?4)
colnames(cmat.sm) <- c(".m2-m1|s1", ".m3-m2|s1", ".m2-m1|s2", ".m3-m2|s2")
contrasts(sm, 4) <- cmat.sm
mod.2 = lmer(y ~ ys + sm + (1 | site) + (0 + ys | site) )
print(summary(mod.2), cor=FALSE)
Linear mixed model fit by REML
Formula: y ~ ys + sm + (1 | site) + (0 + ys | site)
?AIC ?BIC logLik deviance REMLdev
?1377 1745 -621.4 ? ? 1200 ? ?1243
Random effects:
?Groups ? Name ?Variance ?Std.Dev. Corr
?site ? (Intercept) 7.894963 2.80980
?site ? ys1_1 ? ? ? ?7.401678 2.72060
? ? ? ? ?ys1_2 ? ? ? ?4.987042 2.23317 ? 0.691
? ? ? ? ?ys2_1 ? ? ? ?3.228583 1.79683 ? 0.126 ?0.256
? ? ? ? ?ys2_2 ? ? ? ?5.913934 2.43186 ?-0.177 -0.205 ?0.419
? ? ? ? ?ys3_1 ? ? ?11.414706 3.37857 ?-0.239 -0.242 ?0.018 ?0.651
? ? ? ? ?ys3_2 ? ? ?10.336715 3.21508 ? 0.109 ?0.070 ?0.244 ?0.426
0.661
? ? ? ? ?ys4_1 ? ? ? ?4.113868 2.02827 ? 0.392 ?0.180 -0.115 -0.159
-0.015 ?0.521
? ? ? ? ?ys4_2 ? ? ? ?4.801364 2.19120 ? 0.391 ?0.108 -0.268 -0.215
-0.061 ?0.217 ?0.503
? ? ? ? ?ys5_1 ? ? ? ?2.867584 1.69339 ? 0.186 -0.073 -0.356 -0.190
0.107 ?0.373 ?0.097 ?0.268
? ? ? ? ?ys5_2 ? ? ? ?3.745493 1.93533 ?-0.002 -0.107 -0.049 ?0.075
0.174 ?0.224 -0.069 ?0.046 ?0.241
?Residual ? ? ? ? ? ? 0.010422 0.10209
Number of obs: 1800, groups: site, 60
Fixed effects:
? ? ? ? ? ? Estimate Std. Error t value
(Intercept) 13.892119 ? 0.375663 ? 36.98
ys.s2-s1|y1 -4.447045 ? 0.122316 ?-36.36
ys.s2-s1|y2 -4.164777 ? 0.120216 ?-34.64
ys.s2-s1|y3 -3.161133 ? 0.155749 ?-20.30
ys.s2-s1|y4 -3.047036 ? 0.129606 ?-23.51
ys.s2-s1|y5 -3.863464 ? 0.136543 ?-28.29
sm.m2-m1|s1 ?1.988407 ? 0.008335 ?238.55
sm.m3-m2|s1 -0.981789 ? 0.008335 -117.79
sm.m2-m1|s2 -1.990586 ? 0.008335 -238.81
sm.m3-m2|s2 -0.008253 ? 0.008335 ? -0.99
###########################################
On Wed, May 13, 2009 at 6:19 AM, Murray Jorgensen
<maj at stats.waikato.ac.nz> wrote:
I am trying to help a friend towards an appropriate lmer() analysis of her
data and help myself understand how to use lmer syntax in the right way.
She says "...here is what I think is a very standard situation in
ecology....
A normally distributed measure (say bird density) is collected for:
?5 years of data (fixed factor),
?2 seasons within each year (fixed factor),
?3 months within each season (fixed factor)..."
The response is measured at the same six sites at every visit. What follows
is my interpretation of what she might want, which may not be what she
really wants, but this may not matter for the purpose of creating an
example.
I will simulate some data from a known model, then use lmer() to fit the
model and compare the estimated parameters with those used in the
simulation.
I would be interested in comments about
a) ?whether I have simulated the data from my model correctly
and
b) ?whether my lmer syntax is appropriate for my model.
The systematic part of my model is:
?\mu + \beta_{is} + \gamma_{sm}
where i varies over years, s over seasons, and m over months.
(I assume that months are nested within seasons but crossed with years.)
Now let me decide how these coefficients should vary with site. I am
going to assume that the constant term and the year-season term vary
with site but that month-within-season does not.
The random part of the model needs to be generated from a 1 + 5*2 =
11-dimensional multivariate normal. I will assume 1 correlation
parameter, that between two seasons in the same year. All other
correlations are zero.
Values for fixed effects:
Suppose ?\mu = 10
? the beta_{i1} are ? 7, 8, 9, 6, 7.5
? the beta_{i2} are ? 1, 3, 4, 2, 2.5
? the gamma_{1m} are ?0, 2, 1
? the gamma_{2m} are ?0, -2, -2
The s.d. parameters are ?sigma_mu = 3, sigma_beta = 2, sigma = 0.1
The correlation parameter is 0.4.
Instead of 6 sites I will assume nsite = 60.
I have written the following to generate data and fit the model:
library(mvtnorm)
library(lme4)
V = diag(11)
for ( i in 2:10) {
?V[i,i+1] = 0.4
?V[i+1,i] = 0.4
?}
sig = diag(c(3,rep(2,10)))
V = sig%*%V%*%sig
nsite = 60 ? ? ? ? ? #Why not?
sigma = 0.1
theta = c(10,7, 1, 8, 3, 9, 4, 6,2, 7.5,2.5)
P = rmvnorm(nsite,mean=theta,sigma=V)
a = ?rep(1:10,rep(3,10))
TH = ?cbind(rep(1,30),model.matrix(~factor(a)-1))
gg = rep(c(0,2,1,0,-2,-2),5)
y = as.vector(TH%*%t(P)) + rep(gg,nsite) + rnorm(30*nsite,0,sigma)
year = gl(5,2*3,30*nsite)
season = gl(2,3,30*nsite)
month = gl(3,1,30*nsite)
site = gl(nsite,30,30*nsite)
mod = lmer(y ~ year/season + season/month + (1 + year/season | site) )
summary(mod)
The output I get is
library(mvtnorm) library(lme4)
V = diag(11)
for ( i in 2:10) {
+ ? ?V[i,i+1] = 0.4 + ? ?V[i+1,i] = 0.4 + ? ?}
sig = diag(c(3,rep(2,10))) V = sig%*%V%*%sig
nsite = 60 sigma = 0.1 theta = c(10,7, 1, 8, 3, 9, 4, 6,2, 7.5,2.5) P = rmvnorm(nsite,mean=theta,sigma=V) a = ?rep(1:10,rep(3,10)) TH = ?cbind(rep(1,30),model.matrix(~factor(a)-1)) gg = rep(c(0,2,1,0,-2,-2),5) y = as.vector(TH%*%t(P)) + rep(gg,nsite) + rnorm(30*nsite,0,sigma) year = gl(5,2*3,30*nsite) season = gl(2,3,30*nsite) month = gl(3,1,30*nsite) site = gl(nsite,30,30*nsite) mod = lmer(y ~ year/season + season/month + (1 + year/season | site) ) summary(mod)
Linear mixed model fit by REML Formula: y ~ year/season + season/month + (1 + year/season | site) ?AIC ?BIC logLik deviance REMLdev 1332 1717 -596.1 ? ? 1153 ? ?1192 Random effects: Groups ? Name ? ? ? ? ?Variance ?Std.Dev. Corr site ? ? (Intercept) ? 14.874181 3.85671 ? ? ? ? year2 ? ? ? ? ?8.955943 2.99265 -0.234 ? ? ? ? year3 ? ? ? ? ?8.751238 2.95825 ?-0.343 0.504 ? ? ? ? year4 ? ? ? ? ?8.671103 2.94467 ?-0.208 ?0.450 0.422 ? ? ? ? year5 ? ? ? ? ?7.388752 2.71823 ?-0.324 ?0.402 ?0.516 0.462 ? ? ? ? year1:season2 ?4.283068 2.06956 ?-0.156 ?0.638 ?0.322 ?0.517 0.358 ? ? ? ? year2:season2 ?5.046195 2.24637 ?-0.149 -0.469 ?0.223 -0.053 0.015 -0.308 ? ? ? ? year3:season2 ?5.490231 2.34312 ? 0.078 -0.139 -0.418 ?0.243 -0.007 ?0.019 -0.294 ? ? ? ? year4:season2 ?5.353655 2.31380 ? 0.135 -0.204 -0.056 -0.583 0.086 -0.254 ?0.009 -0.318 ? ? ? ? year5:season2 ?5.400134 2.32382 ?-0.053 ?0.047 ?0.068 -0.044 -0.472 -0.068 ?0.074 -0.018 Residual ? ? ? ? ? ? ? ?0.010340 0.10169 -0.301 Number of obs: 1800, groups: site, 60 Fixed effects: ? ? ? ? ? ? ? Estimate Std. Error t value (Intercept) ? ?17.473505 ? 0.497980 ? 35.09 year2 ? ? ? ? ? 0.731252 ? 0.386499 ? ?1.89 year3 ? ? ? ? ? 1.965976 ? 0.382058 ? ?5.15 year4 ? ? ? ? ?-1.014649 ? 0.380309 ? -2.67 year5 ? ? ? ? ? 0.688820 ? 0.351081 ? ?1.96 season2 ? ? ? ?-6.255910 ? 0.267479 ?-23.39 year2:season2 ? 1.889409 ? 0.451063 ? ?4.19 year3:season2 ? 1.117068 ? 0.400037 ? ?2.79 year4:season2 ? 2.564384 ? 0.448766 ? ?5.71 year5:season2 ? 1.055936 ? 0.415366 ? ?2.54 season1:month2 ?1.993234 ? 0.008303 ?240.07 season2:month2 -2.029299 ? 0.008303 -244.41 season1:month3 ?1.001748 ? 0.008303 ?120.65 season2:month3 -2.027852 ? 0.008303 -244.24 Correlation of Fixed Effects: ? ? ? ? ? (Intr) year2 ?year3 ?year4 ?year5 ?seasn2 yr2:s2 yr3:s2 yr4:s2 yr5:s2 ssn1:2 ssn2:2 ssn1:3 year2 -0.235 year3 ? ? ? -0.343 0.504 year4 ? ? ? -0.208 ?0.450 0.422 year5 ? ? ? -0.324 ?0.402 ?0.516 0.462 season2 ? ? -0.157 ?0.637 ?0.322 ?0.516 0.358 year2:sesn2 -0.003 -0.680 -0.047 -0.341 -0.203 -0.790 year3:sesn2 ?0.163 -0.531 -0.531 -0.162 -0.245 -0.654 0.377 year4:sesn2 ?0.183 -0.515 -0.229 -0.696 -0.156 -0.765 ?0.575 0.343 year5:sesn2 ?0.062 -0.376 -0.159 -0.364 -0.572 -0.693 ?0.572 ?0.444 0.377 sesn1:mnth2 -0.008 ?0.000 ?0.000 ?0.000 ?0.000 ?0.016 ?0.000 ?0.000 0.000 ?0.000 sesn2:mnth2 ?0.000 ?0.000 ?0.000 ?0.000 ?0.000 -0.016 ?0.000 ?0.000 0.000 ?0.000 ?0.000 sesn1:mnth3 -0.008 ?0.000 ?0.000 ?0.000 ?0.000 ?0.016 ?0.000 ?0.000 0.000 ?0.000 ?0.500 ?0.000 sesn2:mnth3 ?0.000 ?0.000 ?0.000 ?0.000 ?0.000 -0.016 ?0.000 ?0.000 0.000 ?0.000 ?0.000 ?0.500 ?0.000
The variance estimates for the random effects look large, but maybe because of the way I have parameterized the model I am estimating some combination of sigma_mu and sigma_beta? Cheers, Murray Jorgensen -- Dr Murray Jorgensen ? ? ?http://www.stats.waikato.ac.nz/Staff/maj.html Department of Statistics, University of Waikato, Hamilton, New Zealand Email: maj at waikato.ac.nz ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?Fax 7 838 4155 Phone ?+64 7 838 4773 wk ? ?Home +64 7 825 0441 ? Mobile 021 0200 8350
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models