[R-sig-ME] GLMM for underdispersed count data: Conway-Maxwell-Poisson and Ordinal

[Ben Bolker]

On 16-12-14 11:05 AM, Simone Santoro wrote:
> Thank you all so much for it. Really a very useful discussion, hope it
> may be so to others too.
> 
> El 14/12/2016 a las 16:06, Mollie Brooks escribió:
>> Hi Simone,
>>
>> For the glmmTMB model with the Conway-Maxwell Poisson distribution,
>> the left side of the equation should technically by fledges rather
>> than as.factor(fledges). However, it looks like glmmTMB doesn’t
>> evaluate the as.factor() command and fits the model with fledges as
>> the response anyway.

  I'd be careful with this conclusion.  I think these are *not* the same
model.  For the fake data set (where all the true effects are zero) the
results aren't that different, but what happens when you convert an
integer value to a factor is that the unique values get converted to
codes 1, 2, ...  This would potentially be disastrous.

  I got the clmm model to run with Rune's development version; it's a
little hard to see whether the results are comparable or not since it's
fitting a qualitatively different model ...


>>
>> If you end up needing zero-inflation also, it can be specified using
>> the ziformula command. See vignette("glmmTMB") or here
>> https://github.com/glmmTMB/glmmTMB/blob/master/misc/salamanders.pdf
>> for an example.
>>
>> cheers,
>> Mollie
>>
>> ———————————
>> Mollie E. Brooks, Ph.D.
>> Postdoctoral Researcher
>> National Institute of Aquatic Resources
>> Technical University of Denmark
>>
>>> On 9Dec 2016, at 19:40, Simone Santoro <santoro at ebd.csic.es
>>> <mailto:santoro at ebd.csic.es>> wrote:
>>>
>>> Hi,
>>>
>>> Thank you all very much your hints. They have been really really
>>> helpful for me. Below you may find a reproducible code to see how
>>> three approaches fit a simulated data set (clmm::ordinal,
>>> glmmTMB::glmmTMB, fitme:spaMM). Results seem to me qualitatively
>>> similar but with clmm:ordinal I cannot use the three crossed random
>>> effects because I get an error like this:
>>> Error: no. random effects (=135) >= no. observations (=100)
>>>
>>> set.seed(1234)
>>> library(ordinal)
>>> library(glmmTMB)
>>> library(spaMM)
>>> dati<- data.frame(fledges= rpois(100,10), habitatF=
>>> as.factor(rbinom(100,1,0.5)), areaPatchFath= rnorm(100), poligF01=
>>> as.factor(rbinom(100,1,0.5)),StdLayingDate= rnorm(100), ageFath1=
>>> rpois(100,3), ageMoth1= rpois(100,3), year=
>>> as.factor(rpois(100,200)), ringMoth= as.factor(rpois(100,200)),
>>> ringFath= as.factor(rpois(100,200)))
>>> str(dati)
>>>
>>> system.time(Fitclm<- clmm(as.factor(fledges) ~
>>> habitatF*(areaPatchFath+poligF01+StdLayingDate+ageFath1+ageMoth1)+(1|year)+(1|ringMoth)+(1|ringFath),data=dati,Hess=T))
>>> # this way it works...
>>> system.time(Fitclm1<- clmm(as.factor(fledges) ~
>>> habitatF*(areaPatchFath+poligF01+StdLayingDate+ageFath1+ageMoth1)+(1|year)+(1|ringFath),data=dati,Hess=T))
>>> summary(Fitclm1)
>>>
>>> system.time(FitglmmTMB<- glmmTMB(as.factor(fledges) ~
>>> habitatF*(areaPatchFath+poligF01+StdLayingDate+ageFath1+ageMoth1)+(1|year)+(1|ringMoth)+(1|ringFath),data=dati,family=
>>> "compois"))
>>> summary(FitglmmTMB)
>>>
>>> system.time(FitglmmTMB<- glmmTMB(as.factor(fledges) ~
>>> habitatF*(areaPatchFath+poligF01+StdLayingDate+ageFath1+ageMoth1)+(1|year)+(1|ringMoth)+(1|ringFath),data=dati,family=
>>> "compois"))
>>> summary(FitglmmTMB)
>>>
>>> # This lasts much more (3-4')
>>> system.time(Fitfitme<- fitme(fledges ~
>>> habitatF*(areaPatchFath+poligF01+StdLayingDate+ageFath1+ageMoth1)+(1|year)+(1|ringFath)+(1|ringMoth),data=dati,COMPoisson(),method
>>> = "ML"))
>>> summary(Fitfitme)
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> El 08/12/2016 a las 4:32, Ben Bolker escribió:
>>>>    One reference that uses ordinal regression in a similar situation
>>>> (litter size of Florida panthers) is
>>>> http://link.springer.com/article/10.1007/s00442-011-2083-0 ("Does
>>>> genetic introgression improve female reproductive performance? A test on
>>>> the endangered Florida panther")
>>>>
>>>>   Not sure about the number-of-random-effects error: a reproducible
>>>> example would probably be needed (smaller is better!)
>>>>
>>>>   Ben Bolker
>>>>
>>>>
>>>> On 16-12-06 08:41 AM, Simone Santoro wrote:
>>>>> Dear all,
>>>>>
>>>>> I am trying to find an appropriate GLMM (with temporal and individual 
>>>>> crossed random effects) to model underdispersed count data (clutch 
>>>>> size). I have found several possible ways of doing that. A good 
>>>>> distribution for data like this would seem to be the 
>>>>> Conway-Maxwell-Poisson but I have not found a way of using it within a 
>>>>> GLMM in R (I have asked here 
>>>>> <http://stats.stackexchange.com/questions/249738/how-to-define-the-nu-parameter-of-conway-maxwell-poisson-in-spamm-package> 
>>>>> and here 
>>>>> <http://stats.stackexchange.com/questions/249798/conway-maxwell-poisson-with-crossed-random-effects-in-r>).
>>>>> I have seen that Ben Bolker suggested (here 
>>>>> <https://stat.ethz.ch/pipermail/r-sig-mixed-models/2014q1/021945.html>and 
>>>>> here 
>>>>> <http://stats.stackexchange.com/questions/92156/how-to-handle-underdispersion-in-glmm-binomial-outcome-variable>) 
>>>>> to use an ordinal model in cases like this(e.g. _ordinal:clmm_). I have 
>>>>> tried this solution and the results I obtain makes (biological) sense to 
>>>>> me. However, I wonder why but I cannot put all the three crossed random 
>>>>> effects I have in the clmm model (_Error: no. random effects (=1254) >= 
>>>>> no. observations (=854)_) whereas it is not a problem for the glmer 
>>>>> model (the no. of levels of each single random effect does not exceed 854)*.
>>>>> Beyond that, and that's what I would like to ask you, *I cannot find a 
>>>>> reference to justify I used the ordinal model* to deal with 
>>>>> underdispersed count data (referee will ask it for sure).
>>>>> Best,
>>>>>
>>>>> Simone
>>>>>
>>>>> * FMglmer<- glmer(fledges ~ habitatF * (areaPatchFath + poligF01 + 
>>>>> StdLayingDate + ageFath1 + ageMoth1) + (1|year) + (1|ringMoth) + 
>>>>> (1|ringFath), data = datiDRS)
>>>>>     FMclmm<- glmer(as.factor(fledges)~ habitatF * (areaPatchFath + 
>>>>> poligF01 + StdLayingDate + ageFath1 + ageMoth1) + (1|year) + 
>>>>> (1|ringMoth) + (1|ringFath), data = datiDRS)
>>>>>
>>>>>
>>>>> 	[[alternative HTML version deleted]]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-mixed-models at r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>
>>>
>>> -- 
>>> Simone Santoro
>>> PhD
>>> Department of Ethology and Biodiversity Conservation
>>> Doñana Biological Station
>>> Calle Américo Vespucio s/n
>>> 41092 Seville - Spain
>>> Phone no. +34 954 466 700 (ext. 1213)
>>> http://www.researchgate.net/profile/Simone_Santoro
>>> http://orcid.org/0000-0003-0986-3278
>>
> 
> -- 
> Simone Santoro
> PhD
> Department of Ethology and Biodiversity Conservation
> Doñana Biological Station
> Calle Américo Vespucio s/n
> 41092 Seville - Spain
> Phone no. +34 954 466 700 (ext. 1213)
> http://www.researchgate.net/profile/Simone_Santoro
> http://orcid.org/0000-0003-0986-3278
>