3
$\begingroup$

I'm sorry if this is a duplicate question; I searched around for an answer for some time, but couldn't find anything.

I want to build a model in R, with the proportional number of individuals ('count'/'allCount') as a response variable. 'count' is the number of individuals within a specific group, and 'allCount' the number of individuals with all groups combined. Here is the data, which is a subset for this specific group (a foraging guild of birds, within a larger bird community):

"count","allCount","site","treatment","elevation","field.season"
3,20,"CM","fragmented",800,"2016"
4,19,"CM","fragmented",800,"2016"
0,12,"SM","continuous",800,"2016"
1,21,"CA","fragmented",1200,"2016"
3,11,"SA","continuous",1200,"2016"
3,29,"SA","continuous",1200,"2016"
0,16,"SB","continuous",300,"2016"
12,38,"CB","fragmented",300,"2016"
3,22,"SM","continuous",800,"2016"
2,19,"CM","fragmented",800,"2016"
0,22,"SA","continuous",1200,"2016"
3,28,"CA","fragmented",1200,"2016"
15,39,"CB","fragmented",300,"2016"
3,19,"SB","continuous",300,"2016"
15,38,"CA","fragmented",1200,"2017"
12,22,"CA","fragmented",1200,"2017"
3,10,"CM","fragmented",800,"2017"
1,10,"SM","continuous",800,"2017"
1,11,"SM","continuous",800,"2017"
1,17,"CM","fragmented",800,"2017"
5,17,"SA","continuous",1200,"2017"
4,10,"SA","continuous",1200,"2017"
1,21,"SB","continuous",300,"2017"
1,9,"SB","continuous",300,"2017"
1,14,"CB","fragmented",300,"2017"
2,19,"CB","fragmented",300,"2017"
0,6,"CA","fragmented",1200,"2017"
4,14,"SA","continuous",1200,"2017"
2,12,"CM","fragmented",800,"2017"
1,12,"SM","continuous",800,"2017"
2,11,"SB","continuous",300,"2017"

I have tried three different ways of fitting the data.


1. Binomial with lme4

I started with a binomial family in the package lme4:

model1 <- glmer(count/allCount~ (1|site) + treatment + elevation + field.season + 
treatment*elevation + treatment*field.season, family = "binomial", data=dat)

Curiously, this first model failed:

Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  : 
Downdated VtV is not positive definite

2. Binomial with glmmadmb

I then tried fitting the same model in the package glmmadmb, this time without the error message:

model2 <- glmmadmb(count/allCount~ (1|site) + treatment + elevation + field.season + 
treatment*elevation + treatment*field.season, family = "binomial", data=dat)

After checking overdispersion, it seemed like something was wrong with this model, too:

> overdisp.glmer(model2)
 Residual deviance: 1848789055.502 on 24 degrees of freedom (ratio: 77032877.313)

Here the model fit: enter image description here


3. Negative binomial

Finally, I tried fitting the model with a negative binomial distribution. This time, it converged without problems, and the residuals seemed well dispersed:

model3 <- glmmadmb(count/allCount*100~ (1|site) + treatment + elevation + field.season + 
          treatment*elevation + treatment*field.season, family = "nbinom", data=dat)
> overdisp.glmer(model3)
  Residual deviance: 21.87 on 24 degrees of freedom (ratio: 0.911)

enter image description here


My Question:

  1. Is there something conceptually wrong with choosing the negative binomial distribution in this case?
  2. Why did my first model fail in the first place / why did the binomial model look so bad?
$\endgroup$
  • 1
    $\begingroup$ The error message has nothing to do with validity of your model. $\endgroup$ – Tim May 27 '17 at 19:08
  • $\begingroup$ Thanks, that's good to know! Do you know what could have caused the error? $\endgroup$ – Joris May 27 '17 at 19:16
6
$\begingroup$

There is one big conceptual error here and one procedural/programming error (at least).

  1. (conceptual) in general, the negative binomial distribution is not appropriate for proportion data, or count data with an upper bound. Binomial (or beta-binomial, if you want to allow for additional variation) are appropriate for proportions expressed as counts; Beta is appropriate for proportions with an unknown denominator. You can use count models (Poisson or negative binomial) with an offset if the probability of a successful outcome is low, i.e. the counts are always much less than the maximum possible.
  2. (procedural) if expressing a binomial response as a proportion in glm or glmer, you must include the total number via the weights argument, e.g. model1 <- glmer(count/allCount~ (1|site) + treatment*(elevation + field.season) + family = "binomial", weights=allCount, data=dat)
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Thank you for answering both of my questions! After adding the weights argument, the model runs without producing an error. $\endgroup$ – Joris May 28 '17 at 8:22
  • $\begingroup$ Hi @Ben Bolker - I apologise for responding to an old answer, but do you know if there is a good reference for the conceptual concerns? I've come up against a reviewer saying something similar and need to justify my use of quasibinomial as opposed to negative binomial $\endgroup$ – NatWH Jun 27 '19 at 18:33
  • $\begingroup$ Sorry, don't know of anything offhand. (You could ask a new question referring to this one and specifically looking for literature references ...) $\endgroup$ – Ben Bolker Jun 27 '19 at 22:38

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.