5
$\begingroup$

I'm considering a mixed-effects model to try to understand factors that influence the number of ticks sampled on wild rodents. My data is nested so that I have one tick count per rodent, multiple rodents per site and multiple sites per year (sites are repeated across years but not every site is present in all years).

So far, without fixed effects, my model looks something like this:

glmer(Ticks~(1|Year)+(1|Site), family=poisson, data=tickdata)

I understand that this will account for random variation between sites and between years. My main concern is that I think the effect of specific sites will change between years, e.g. some sites may recieve more rainfall than others, but the identity of the sites with the highest rainfall differs between years.

So do I need some sort of interaction term? Maybe something like this:

glmer(Ticks~(1|Year)+(1|Site)+(1|Year:Site), family=poisson, data=tickdata) 
$\endgroup$
  • $\begingroup$ I wonder if it makes more sense to nest year at the bottom of the hierarchy: observation in year nested in ticks, nested in rodents, nested in sites? $\endgroup$ – Alexis May 6 '14 at 14:28
  • $\begingroup$ @Alexis, you don't need to nest in this case (and probably shouldn't) $\endgroup$ – Ben Bolker May 6 '14 at 21:35
  • $\begingroup$ Because...? Are there no year-specific determinants of tick prevalence? Are there no site-specific determinants of tick prevalence? Rodent-specific? Ignoring such structural heterogeneity of variation surely guarantees biased estimates of the standard errors? $\endgroup$ – Alexis May 6 '14 at 21:49
  • 1
    $\begingroup$ I mean that you can cross random effects rather than nesting: exactly as the OP said, there are site-specific effects, year-specific effects, and site-by-year effects. (You could also consider rodent-specific effects, which would be observation-level effects/characterizing overdispersion if each rodent is measured once. Otherwise (if rodents are recaptured/recounted) you could also include observation-level effects). $\endgroup$ – Ben Bolker May 7 '14 at 1:43
7
$\begingroup$

Have you tried it? That sounds like it should be fine.

set.seed(101)
## generate fully crossed design:
d <- expand.grid(Year=2000:2010,Site=1:30)
## sample 70% of the site/year comb to induce lack of balance
d <- d[sample(1:nrow(d),size=round(0.7*nrow(d))),]
## now get Poisson-distributed number of obs per site/year
library(plyr)
d <- ddply(d,c("Site","Year"),transform,rep=seq(rpois(1,lambda=10)))
library(lme4)
d$ticks <- simulate(~1+(1|Year)+(1|Site)+(1|Year:Site),
                    family=poisson,newdata=d,
                    newparams=list(beta=2, ## mean(log(ticks))=2
                               theta=c(1,1,1)))[[1]]
mm <- glmer(ticks~1+(1|Year)+(1|Site)+(1|Year:Site),
                    family=poisson,data=d)

We get out approximately what we put in -- equal variances at each level:

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: ticks ~ 1 + (1 | Year) + (1 | Site) + (1 | Year:Site)
##    Data: d
## 
##      AIC      BIC   logLik deviance df.resid 
##  12487.3  12510.2  -6239.7  12479.3     2267 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9944 -0.6842 -0.0726  0.6010  3.8532 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Year:Site (Intercept) 1.0818   1.0401  
##  Site      (Intercept) 1.0490   1.0242  
##  Year      (Intercept) 0.9787   0.9893  
## Number of obs: 2271, groups:  Year:Site, 231 Site, 30 Year, 11
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.1952     0.3593   6.109    1e-09 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

You may want to include an observation-level random effect to allow for overdispersion (see the "grouse ticks" example in http://rpubs.com/bbolker/glmmchapter)

$\endgroup$
  • $\begingroup$ Trying your code here, I get an error: Error in UseMethod("simulate") : no applicable method for 'simulate' applied to an object of class "formula". Any thoughts on what that may be about? $\endgroup$ – gung Jan 13 '15 at 17:57
  • $\begingroup$ are you sure (1) you loaded lme4 (2) you have a recent version of lme4 ? $\endgroup$ – Ben Bolker Jan 13 '15 at 18:11
  • $\begingroup$ Seems I had an older version of lme4. Sorry to bother you. $\endgroup$ – gung Jan 13 '15 at 18:19
  • $\begingroup$ Why did we need to sample 70% of the site/year comb to induce lack of balance here? $\endgroup$ – Joshua Rosenberg Apr 8 '17 at 13:36
  • $\begingroup$ I no longer remember why I did that. $\endgroup$ – Ben Bolker Apr 8 '17 at 14:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.