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
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
