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
