On 03/03/14 17:45, Howe, Eric (MNR) wrote:
> 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.
It looks like you're ignoring any overdispersion, so you're assuming
that the variance equals the mean. For example, a mean of 900 has a
standard deviation of 30. This is smaller than the range of your data.
Because of this, the model assumes there is almost no residual error,
which forces teh effects to have a small error.
With counts this high I would start by simply square root transforming
the data and analysing that: the variance should be fairly well
stabilised and the residuals should be pretty much normally distributed.
Alternatively, you could add overdispersion by using either a negative
binomial distribution, or an extra random effect with a level for each
observation (this would be equivalent to an Poisson log-normal assumption).
Bob
> 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]]
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
--
Bob O'Hara
Biodiversity and Climate Research Centre
Senckenberganlage 25
D-60325 Frankfurt am Main,
Germany
Tel: +49 69 7542 1863
Mobile: +49 1515 888 5440
WWW: http://www.bik-f.de/root/index.php?page_id=219
Blog: http://blogs.nature.com/boboh
Journal of Negative Results - EEB: www.jnr-eeb.org
[[alternative HTML version deleted]]
_______________________________________________
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
