[R] mixed-effects model using lmer
hlynch at umd.edu
hlynch at umd.edu
Thu Jun 28 22:44:31 CEST 2007
Hello R-users,
I have been trying to fit what I think is a simple mixed-effects model using lmer (from lme4), but I've run into some difficulty that I have not been able to resolve using the existing archives or Pinheiro and Bates (2000).
I am measuring populations (of birds) which change with time at a number of different sites. These sites are grouped into regions. Sites are not measured with any regularity, so each site has a different number of observations. I want to treat the effect of year (i.e. the rate of population change) as a fixed effect which is modeled separately for each region, but I want to allow sites within regions to be random effects. I've come up with the following model, but I'm not sure its correct. Any help would be most appreciated.
fit<-lmer(log.count~site+month+year:region-1+(year:region-1|site))
What my variables mean:
log.count=log(count)
site=factor for the different sites, each site is associated with exactly one region
month=fixed effect for the month of the survey
year:region = allows me to model different year effects for each region
Each site should have a completely separate intercept, so I have chosen to suppress the intercept and use all the factors of site instead.
The difficulty is, of course, with the last term. I want to have sites vary randomly *within* a region, but I don't want them to vary randomly throughout all the regions, i.e. I want the random effects to sum to zero for each region separately.
My output looks like (using the display() function from Gelman and Hill (2007)):
*******************************
lmer(formula = log.count ~ site + month + year:region - 1 + (year:region -
1 | site))
coef.est coef.se
siteAITC 7.25 0.15
siteALMI 4.57 0.29
siteBENE 5.98 0.29
siteBOOT 7.08 0.17
siteBROW 6.47 0.14
(I cut some out here for space)
siteUSEF 6.97 0.21
siteWATE 7.38 0.15
siteYANK 8.31 0.11
monthJAN -0.15 0.07
monthNOV -0.01 0.07
year:regionNE -0.16 0.12
year:regionNW 0.04 0.01
year:regionSH 0.02 0.02
year:regionSW 0.07 0.01
Error terms:
Groups Name Std.Dev. Corr
site year:regionNE 0.12
year:regionNW 0.03 0.00
year:regionSH 0.00 0.00 0.00
year:regionSW 0.00 0.00 0.00 0.00
Residual 0.26
---
number of obs: 110, groups: site, 23
AIC = 157.8, DIC = -72.6
deviance = 3.6
********************************
The coefficients look totally reasonable given all the analysis I've already done on the system. My confusion comes about when I look at the random effects using ranef(fit) and se.ranef(fit):
********************************
>ranef(fit)
An object of class "ranef.lmer"
[[1]]
year:regionNE year:regionNW year:regionSH year:regionSW
AITC 0.0e+00 0.0e+00 5.8e-10 0.0e+00
ALMI 0.0e+00 -2.3e-04 0.0e+00 0.0e+00
BENE 0.0e+00 -2.1e-17 0.0e+00 0.0e+00
BOOT 0.0e+00 0.0e+00 0.0e+00 -2.4e-09
BROW 3.0e-15 0.0e+00 0.0e+00 0.0e+00
BRYE 0.0e+00 -1.3e-02 0.0e+00 0.0e+00
BRYS 0.0e+00 -8.9e-17 0.0e+00 0.0e+00
CUVE 0.0e+00 -1.7e-02 0.0e+00 0.0e+00
DAMO 0.0e+00 -5.0e-03 0.0e+00 0.0e+00
DANC 0.0e+00 -5.2e-03 0.0e+00 0.0e+00
DOBE 0.0e+00 8.1e-04 0.0e+00 0.0e+00
GEOR 0.0e+00 -9.3e-03 0.0e+00 0.0e+00
HANN 0.0e+00 0.0e+00 3.2e-10 0.0e+00
HERO -6.3e-16 0.0e+00 0.0e+00 0.0e+00
JOUG 0.0e+00 -2.3e-02 0.0e+00 0.0e+00
MADD 7.3e-16 0.0e+00 0.0e+00 0.0e+00
MOOT 0.0e+00 0.0e+00 0.0e+00 8.7e-11
NEKO 0.0e+00 -2.1e-03 0.0e+00 0.0e+00
PETE 0.0e+00 0.0e+00 0.0e+00 2.5e-09
PLEN 0.0e+00 0.0e+00 0.0e+00 -2.0e-10
USEF 0.0e+00 5.6e-02 0.0e+00 0.0e+00
WATE 0.0e+00 1.8e-02 0.0e+00 0.0e+00
YANK 0.0e+00 0.0e+00 -9.0e-10 0.0e+00
>se.ranef(fit)
$site
year:regionNE year:regionNW year:regionSH year:regionSW
AITC 0.120 0.0296 5.8e-06 5.8e-06
ALMI 0.120 0.0204 5.8e-06 5.8e-06
BENE 0.120 0.0196 5.8e-06 5.8e-06
BOOT 0.120 0.0296 5.8e-06 5.8e-06
BROW 0.016 0.0296 5.8e-06 5.8e-06
BRYE 0.120 0.0144 5.8e-06 5.8e-06
BRYS 0.120 0.0231 5.8e-06 5.8e-06
CUVE 0.120 0.0136 5.8e-06 5.8e-06
DAMO 0.120 0.0086 5.8e-06 5.8e-06
DANC 0.120 0.0211 5.8e-06 5.8e-06
DOBE 0.120 0.0221 5.8e-06 5.8e-06
GEOR 0.120 0.0153 5.8e-06 5.8e-06
HANN 0.120 0.0296 5.8e-06 5.8e-06
HERO 0.070 0.0296 5.8e-06 5.8e-06
JOUG 0.120 0.0086 5.8e-06 5.8e-06
MADD 0.047 0.0296 5.8e-06 5.8e-06
MOOT 0.120 0.0296 5.8e-06 5.8e-06
NEKO 0.120 0.0148 5.8e-06 5.8e-06
PETE 0.120 0.0296 5.8e-06 5.8e-06
PLEN 0.120 0.0296 5.8e-06 5.8e-06
USEF 0.120 0.0137 5.8e-06 5.8e-06
WATE 0.120 0.0226 5.8e-06 5.8e-06
YANK 0.120 0.0296 5.8e-06 5.8e-06
**************************************
I'm not sure this is correct. For one thing, I get an warning when I run the model that "Estimated variance-covariance for factor 'site' is singular". Also, although the random effects do add to zero for each region (within rounding error), I'm not sure how to interpret their standard errors. Each site is in only one of the four regions, and yet lmer seems to be estimating coefficients and standard errors as though each site was replicated in each region which is not the case. I have tried several other strategies based on suggestions from the archive, but this is the closest I have come to something reasonable.
Thanks in advance for the help,
Heather Lynch
More information about the R-help
mailing list