[R-sig-ME] Binomial GLMM vs GLM question

[Justin Touchon]

Hi all,
Thank you for all for you informative and helpful comments about my 
dilemma. I've pasted below a summary table of egg survival at the two 
ponds during the three years and the GLM output for a model looking for 
interactions between pond, rainfall and year. Below that is the output 
of a LMER model Ben suggested. Also, as Ben wondered earlier, the data 
are highly overdispersed. To give a little background, this wasn't an 
experiment but an observational study of what was happening in nature. 
The two ponds were selected because they are near my field site, no 
other reason. They differ in the amount of canopy shade overhead, which 
in turn is clearly going to affect egg desiccation (the primary source 
of mortality). Rainfall did vary during the three years (it was much 
heavier in 2005 than in 2003 or 2004). As a result, mortality was higher 
at Bridge Pond (little shade) than at Ocelot Pond, and it varied between 
years. These points are evident from the GLM anova table. All three 
factors, significantly affect egg mortality, pond being the most 
important. Rainfall varied across the three years, and there is a 3-way 
interaction which indicates to me simply that the effects of habitat 
(pond) and rainfall are not consistent across time. The model estimates 
make sense (e.g. rainfall has a negative slope, the estimate for Ocelot 
Pond is highly negative, etc.) In other words, the GLM output does make 
biological sense given what we know about the data and lends itself to 
fairly easy interpretation. Regarding the LMER output, it seems the 
variance of both rain and pond is obscenely huge. Is that a result of 
only have 2 ponds and 3 years? The estimate for the rainfall effect is 
negative, like in the GLM, which makes sense as well (as rain increase, 
mortality decreases). Any opinions on where to go from here? Again, 
thank you all for you helpful advice. I am very appreciative.
-Justin


Egg Survival
Year Bridge Ocelot Total
2003 0.069 0.453 0.251
2004 0.115 0.429 0.272
2005 0.270 0.505 0.376
Total 0.141 0.460 0.292


 > GLM.2 <- glm(tbl_mort ~ Pond *total_rainfall *Year , 
family=quasibinomial(logit), data=FieldData0305)

 > summary(GLM.2)

Call:
glm(formula = tbl_mort ~ Pond * total_rainfall * Year, family = 
quasibinomial(logit),
data = FieldData0305)

Deviance Residuals:
Min 1Q Median 3Q Max
-17.965 -4.838 1.062 5.197 15.225

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.68318 0.47947 5.596 4.53e-08 ***
Pond[T.Ocelot] -2.30872 0.53614 -4.306 2.18e-05 ***
total_rainfall -0.11465 0.02333 -4.915 1.39e-06 ***
Year[T.2004] -0.97218 0.63583 -1.529 0.127198
Year[T.2005] -3.92557 0.76808 -5.111 5.38e-07 ***
Pond[T.Ocelot]:total_rainfall 0.09799 0.02609 3.756 0.000204 ***
Pond[T.Ocelot]:Year[T.2004] 0.48230 0.77189 0.625 0.532506
Pond[T.Ocelot]:Year[T.2005] 2.15914 1.22853 1.757 0.079739 .
total_rainfall:Year[T.2004] 0.08475 0.02717 3.120 0.001967 **
total_rainfall:Year[T.2005] 0.12082 0.02440 4.952 1.16e-06 ***
Pond[T.Ocelot]:total_rainfall:Year[T.2004] -0.05511 0.03585 -1.537 0.125162
Pond[T.Ocelot]:total_rainfall:Year[T.2005] -0.11371 0.02948 -3.857 
0.000137 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 38.23861)

Null deviance: 20762 on 349 degrees of freedom
Residual deviance: 15252 on 338 degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

 > Anova(GLM.2, test="F")
Anova Table (Type II tests)

Response: tbl_mort
SS Df F Pr(>F)
Pond 1002.7 1 26.2229 5.119e-07 ***
total_rainfall 238.2 1 6.2306 0.0130318 *
Year 496.3 2 6.4894 0.0017157 **
Pond:total_rainfall 101.7 1 2.6599 0.1038378
Pond:Year 172.7 2 2.2587 0.1060579
total_rainfall:Year 650.3 2 8.5037 0.0002493 ***
Pond:total_rainfall:Year 704.8 2 9.2159 0.0001267 ***
Residuals 12924.5 338
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


 > m1 <- lmer(tbl_mort~total_rainfall + (1|Year)+(1|Pond), family = 
quasibinomial, data=FieldData0305)

 > summary(m1)
Generalized linear mixed model fit using Laplace
Formula: tbl_mort ~ total_rainfall + (1 | Year) + (1 | Pond)
Data: FieldData0305
Family: quasibinomial(logit link)
AIC BIC logLik deviance
14657 14673 -7325 14649
Random effects:
Groups Name Variance Std.Dev.
Year (Intercept) 2.0141e+15 44878595
Pond (Intercept) 1.3427e+15 36643219
Residual 8.8116e+16 296844003
number of obs: 350, groups: Year, 3; Pond, 2

Fixed effects:
Estimate Std. Error t value
(Intercept) 7.486e-01 5.736e+13 1.305e-14
total_rainfall -2.316e-02 4.406e+11 -5.257e-14

Correlation of Fixed Effects:
(Intr)
total_rnfll -0.956






> On 17/05/2008, at 6:53 AM, Ben Bolker wrote:
>
> <snip>
> >/
> />/ | I think this is a good example of what seems a common problem for
> />/ | people using mixed models - How to decide what are the random  
> />/ factors
> />/ | and what are the fixed factors.
> />/
> />/ ~  I agree that this is difficult, and contentious (see e.g.
> />/ Andrew Gelman, "Analysis of variance--why it is more important than
> />/ ever," Annals of Statistics 33, no. 1 (2005): 1-53,
> />/ doi:http://dx.doi.org/doi:10.1214/009053604000001048)
> />/
> /Available at
> http://www.stat.columbia.edu/~gelman/research/published/AOS259.pdf <http://www.stat.columbia.edu/%7Egelman/research/published/AOS259.pdf>
>
> >/
> />/ | In your study you have taken repeated measurements on the same  
> />/ ponds,
> />/ | as such, observations from each pond across years are not  
> />/ independent
> />/ | (observations in pond 1 at t+1 will be correlated with observations
> />/ | at t).  If you were only interested in the effect of rainfall on egg
> />/ | survival I would suggest you use a mixed model with year and pond as
> />/ | random effects (based on a similar eg in Crawley's R Book p605). In
> />/ | this case I assume you have measurements on every pond in every year
> />/ | and as such these random effects are crossed (not nested).
> />/
> />/ ~  (a) I would suggest that a mixed model is NOT going to work well
> />/ in this case, even if year and pond are "philosophically" random  
> />/ effects
> />/ (i.e., you don't really care about what happens in those specific
> />/ ponds, and you may even have chosen them with a random-number  
> />/ generator
> />/ out of a list of all possible ponds -- although this is much less
> />/ likely with years ...).  The technical problem is that estimating
> />/ variances from 2 or 3 points is nasty.  This translates into
> />/ inference/philosophical terms because these few points really
> />/ don't give you the data to generalize about the population, even if
> />/ you want to.  (I think one of the confusions is that in the classical
> />/ method-of-moments world there's nothing that says you can't have 2
> />/ denominator degrees of freedom -- your power will be terrible, but
> />/ the expressions won't blow up on you [unless you get negative variance
> />/ estimates ...])
> />/ |
> />/ | m1 <- lmer(mort~rainfall + (year|pond), family = binomial, data=
> />/ | FieldData0305)
> />/
> />/ ~   (1|year)+(1|pond) might work OK, and be slightly more
> />/ parsimonious (OR ponds nested within years, (1|year)+(1|pond:year), or
> />/ vice versa -- and if you're not dead set on treating the random  
> />/ effects
> />/ as "effects of year" and "effects of pond", but were willing to  
> />/ treat it
> />/ as "effects of pond within year" that would buy you the flexibility
> />/ to use some other package that wasn't so good at crossed
> />/ effects.)
> />/
> />/
> />/ | summary(m1)
> />/ |  From the summary output you can assess the influence of the fixed
> />/ | effect to see if the estimate for rainfall (slope estimate) is
> />/ | different to zero.
> />/ |
> />/ | Even if the number of groups (years and ponds) is small it is still
> />/ | better to take account of this group variation than ignore it (see
> />/ | Gelman and Hill 200&.
> />/
> />/ ~  Yes, but ... what if all you can get out of the model is that
> />/ the estimate is nearly zero?  If you go Bayesian instead you can
> />/ deal with some problems by setting an informative prior [in theory
> />/ setting a proper prior, no matter how weak, would generally solve
> />/ the problem, but in reality I suspect that if the model is
> />/ overparameterized you're going to have nasty technical difficulties
> />/ even if the model is theoretically OK])
> />/ ~  With respect to Andrew Gelman (I can't believe I'm saying this --
> />/ sacrilege) but I think he's used to really big social-science data  
> />/ sets,
> />/ where "leave stuff in when you're not sure" is more generally
> />/ a good idea than it is in typical field ecology data sets, where
> />/ the bias-variance tradeoff bites harder.  (At least there are 350
> />/ data points here,
> />/
> />/ |  From this model you can also get a feeling for the between year and
> />/ | between pond variation by looking at the random effects estimates  
> />/ and
> />/ | variances. Obviously you will have some idea if you just plot the
> />/ | means for each year and each pond.
> />/
> />/ ~  Only if the variance estimates don't suck.
> />/
> />/ | If you are also interested in understanding how egg survival differs
> />/ | between ponds or between years and if this interacts with rainfall
> />/ | then it become less straight forward.
> />/ |
> />/ | I would suggest that you try
> />/ |
> />/ | m2 <- lmer(mort~rainfall * year * pond + (year|pond), family =
> />/ | binomial, data= FieldData0305)
> />/ | summary(m2)
> />/ |
> />/ | In the summary you will have estimates for all three factors and
> />/ | their interactions and you can ascertain if these are good
> />/ | explanatory variables for egg survival.
> />/
> />/ ~   Isn't this overparameterized?  We have a fixed effect for
> />/ each year:pond combination (and variation in the slopes of
> />/ the effect with respect to rainfall), as well as a random-effect
> />/ level for years and ponds?
> />/
> /
> This does bring up the important point as to whether the effect of  
> rainfall does vary within year and site. One advantage of the mixed  
> effects method for dealing with varying slopes is it gives a nice  
> population estimate. Unfortunately this is a lot of random effects for  
> not many groups.  It is looking very messy.
>
> I think with this sort of data, as in it has very few groups) the best  
> option is to start with the basics (probably always the best option).  
> Fit a separate model for each group and determine an appropriate  
> transform for rainfall. Do the models fit reasonably? Look at the  
> parameter estimates with 95% CI (plotting is a good idea) and see how  
> much variation there is, and there is nothing strange happening like  
> one unusual year-pond. Then decide whether the more complex models  
> will tell you anything more.
>
> Ken
> >

-- 
Justin Touchon
Doctoral Student
Boston University
Ecology, Behavior and Evolution
http://people.bu.edu/jtouchon