Good day R.sig.ecology members.

We seem to have detected highly significant fixed effects from a GLMM fit using glmer, despite having a fairly small sample size (30 observations, 4 groups in the random factor).  We are hoping that group members can either verify that we are interpreting these results correctly or tell us where we've gone wrong.  We are working in R for Windows version 3.0.2 and using the latest version of lme4.

We sought to simultaneously test for (1) effects of annual variation in natural food availability for black bears (measured using a rank bear food index; bfi) on the number of occurrences of human-bear conflict (oc), and (2) a temporal trend (across years) in oc.  oc and bfi were measured annually in different administrative districts (dist) within the larger region of interest.  Not all districts provided data in all years.

Our data look like this:
   dist   yr lyr  bfi   oc
1   ban 2004   1 2.30  852
2  pemb 2004   1 2.70  105
3    ps 2004   1 1.10  969
4   ban 2005   2 2.05  657
5  pemb 2005   2 1.86  269
6    ps 2005   2 1.69 1119
7   ban 2006   3 2.57  275
8   mid 2006   3 2.66  271
9  pemb 2006   3 2.17   57
10   ps 2006   3 2.08  429
11  ban 2007   4 1.75  530
12  mid 2007   4 2.17  606
13 pemb 2007   4 2.03  189
14   ps 2007   4 1.23 1567
15  ban 2008   5 2.86  363
16  mid 2008   5 3.00  467
17 pemb 2008   5 3.13   95
18   ps 2008   5 2.32  770
19  ban 2009   6 2.19  672
20  mid 2009   6 1.47  653
21 pemb 2009   6 2.03  471
22   ps 2009   6 1.90 1522
23  ban 2010   7 2.00  409
24  mid 2010   7 2.10  640
25 pemb 2010   7 2.10  176
26   ps 2010   7 1.90 1218
27  ban 2011   8 2.50  466
28  mid 2011   8 2.40  462
29 pemb 2011   8 2.40  117
30   ps 2011   8 2.00  963

We used a GLMM to test for effects of bfi and lyr (a dummy variable representing the # of years since the beginning of the study), with district as the only random effect.  The log-link function improved linearity.

Code for fitting the model, and the summary of the model, appear below:

> oc.sr.bfi.yr = glmer(oc ~ bfi + lyr + (1|dist), data = dat, family = poisson(link = log))
> summary(oc.sr.bfi.yr)

Generalized linear mixed model fit by maximum likelihood ['glmerMod']
 Family: poisson ( log )
Formula: oc ~ bfi + lyr + (1 | dist)
   Data: dat

      AIC       BIC    logLik  deviance
1960.3134 1965.9182 -976.1567 1952.3134

Random effects:
 Groups Name        Variance Std.Dev.
 dist   (Intercept) 0.3006   0.5483
Number of obs: 30, groups: dist, 4

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  7.040105   0.277215  25.396   <2e-16 ***
bfi         -0.486053   0.019922 -24.398   <2e-16 ***
lyr          0.035638   0.003598   9.905   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
    (Intr) bfi
bfi -0.132
lyr -0.019 -0.295

The results indicate a negative effect of bfi (more conflict when natural foods are scarce, as expected) and a slight increasing trend across years.  However, we were surprised by the reported p-values for the fixed effects, given our relatively small sample size.

We found the "Getting p-values for fitted models" help page for the lme4 package.

We fitted parameter-reduced models for comparison:
> oc.sr.null = glmer(oc ~ 1 + (1|dist), data = dat, family = poisson(link = log))
> oc.sr.bfi = glmer(oc ~ bfi + (1|dist), data = dat, family = poisson(link = log))

AIC and BIC values were greater by 100 when we excluded the effect of year, and were greater by 500 when we also excluded the effect of the bfi.  These results seem to agree with the summary of the more general model, but they also seem a bit extreme given our small sample size.

Model comparisons using PBmodcomp from the pbkrtest package also seem to support the inclusion of both fixed effects.
> mc.oc.bfi.nobfi = PBmodcomp(oc.sr.bfi, oc.sr.null, nsim = 100)
> mc.oc.yr.noyr = PBmodcomp(oc.sr.bfi.yr, oc.sr.bfi, nsim = 100)
> mc.oc.bfi.nobfi
Parametric bootstrap test; time: 26.13 sec; samples: 100 extremes: 0;
large : oc ~ bfi + (1 | dist)
small : oc ~ 1 + (1 | dist)
         stat df   p.value
LRT    513.43  1 < 2.2e-16 ***
PBtest 513.43     0.009901 **
---

> mc.oc.yr.noyr
Parametric bootstrap test; time: 39.83 sec; samples: 100 extremes: 0;
large : oc ~ bfi + lyr + (1 | dist)
small : oc ~ bfi + (1 | dist)
         stat df   p.value
LRT    98.323  1 < 2.2e-16 ***
PBtest 98.323     0.009901 **
---

We calculated bootstrap confidence intervals using "confint".  They seem to indicate that the estimated effects were unambiguous, but we were unsure whether our sample size was adequate to apply this method.

> confint (oc.sr.bfi.yr, method = c("boot"))
Computing bootstrap confidence intervals ...
                          2.5 %      97.5 %
sd_(Intercept)|dist  0.13871500  0.83969172
(Intercept)          6.54440740  7.56234905
bfi                 -0.53389761 -0.44979598
lyr                  0.02835891  0.04325712

Simple diagnostic plots (residuals vs fitted values and a normal QQ plot of residuals) of the GLMM with effects of both bfi and lyr did not reveal violations of assumptions.  We would be happy to report significant effects of food availability and a significant trend across years, but are concerned that we may not be interpreting these results correctly.  Any help or suggestions would be greatly appreciated.

Sincerely,
Eric


	[[alternative HTML version deleted]]

