[R-sig-ME] About glmer with a quasi-Poisson
Jarrod Hadfield
j.hadfield at ed.ac.uk
Sat Apr 25 20:45:58 CEST 2009
Hi,
The zero-inflated Poisson is specified as a bivariate model, so the
residual and random effect models usually follow a bivariate structure
such as us(trait):site or idh(trait):units. These specify 2x2
(co)variance matrices, in the first a covariance is estimated and in
the second it is set to zero. For ZIP models I would use the idh
function because the covariance between the Poisson process and the
zero-inflation process cannot be estimated from the data. The prior
term would have something like:
prior=list(R=list(V=diag(2), n=2, fix=2), G=list(G1=list(V=diag(2), n=2),
G2=list(V=diag(2), n=2), G3=list(V=diag(2), n=2)))
or if you don't want to have a random effect for the zero-inflation
process have
something like:
prior=list(R=list(V=diag(2), n=2, fix=2), G=list(G1=list(V=diag(c(1,0.000001))
, n=2, fix=2), G2=list(V=diag(c(1,0.000001))
, n=2, fix=2), G3=list(V=diag(c(1,0.000001))
, n=2, fix=2))
Don't fix the residual variance to 0.000001 because the chain will not
then mix.
As long as n is greater then or equal to the dimension of the
covariance matrix the prior is proper and the covariance matrix should
not become ill-conditioned.
Bear in mind that I just made these priors up - you need to pick
something sensible, although the residual variance for the
zero-inflation process (trait 2) is immaterial because there is no
info in the data regarding this variance as with standard binary
variables.
Cheers,
Jarrod
Quoting Strubbe Diederik <diederik.strubbe at ua.ac.be>:
> Dear all,
>
> Thank you for the clarifications on the correct use and coding of
> nested random effects. I decided to try the MCMCglmm package for my
> analysis, but encountered some other problems.
>
> First: dataset of 57 points (each point has a unique ID (a,b,c,...)
> and is nested in a certain site (26 sites, Site1,Site2,...). Points
> have been visited 5 times (v1,v2,...,v5). For simplicity, I want to
> relate the number of birds counted to one environmental variable,
> called pca1. Using most default settings, the model should look like:
>
> MCMCglmm(abundance~data$pca1,family="zipoisson",random=~site+point+visit,data=data,verbose=TRUE)
>
> However, this results in the following error:"Error in
> MCMCglmm(abundance ~ data$pca1, family = "zipoisson", random = ~site
> + : R-structure does not define unique residual for each data
> point". [[ Note that if I specify another distribution, e.g. an
> inappropriate Poisson, this error does not occur]]. As far as I
> understand, this error relates to the rcov part of the model and I
> updated the syntax to:
>
> MCMCglmm(abundance~data$pca1,family="zipoisson",rcov=~us(trait):units,random=~site+point+visit,data=data,verbose=TRUE)
>
> However, this results in another error:
> MCMC iteration = 0
> E[MH acceptance ratio] = 0.000431
> MCMC iteration = 1000
> E[MH acceptance ratio] = 0.439460
> MCMC iteration = 2000
> E[MH acceptance ratio] = 0.447993
> Error in MCMCglmm(abundance ~ data$pca1, family = "poisson", random
> = ~site + :
> ill-conditioned G/R structure: use proper priors if you haven't or
> rescale data if you have
>
> From this I understand that I should define appropriate priors.
> Based on the Plodia tutorial, I tried 3 priors:
> halfV = var(abundance)/2
> prior=list(R=list(n=1,V=halfV),G=list(G1=list(n=1,V=halfV)))
> ##halfV = 0.08639657##
> prior = list(R=list(V=1,n=1e+06),G=list(G1=list(V=1,n=0)))
> prior = list(R=list(V=1,n=0,fix=1),G=list(G1=list(V=1, + n=0)))
>
> However, all of them result in the following error:
> Error in MCMCglmm(abundance ~ data$pca1, family = "zipoisson", rcov
> = ~us(trait):units, :
> ill-conditioned G/R structure: use proper priors if you haven't or
> rescale data if you have
>
> Could anyone offer some advice on what my mistakes are? Or suggested
> papers to read, as this is my first encounter with Bayesian based
> analysis.
>
> Many thanks,
>
> Diederik
>
>
> Diederik Strubbe
> Evolutionary Ecology Group
> Department of Biology, University of Antwerp
> Universiteitsplein 1
> B-2610 Antwerp, Belgium
> http://webhost.ua.ac.be/deco
> tel : 32 3 820 23 85
>
>
>
> -----Original Message-----
> From: Jarrod Hadfield [mailto:jhadfiel at staffmail.ed.ac.uk]
> Sent: Thu 23-4-2009 9:33
> To: Strubbe Diederik
> Cc: r-sig-mixed-models at r-project.org; boristmr at yahoo.com
> Subject: Re: [R-sig-ME] About glmer with a quasi-Poisson
>
>
>
> Diederik Strubbe
> Evolutionary Ecology Group
> Department of Biology, University of Antwerp
> Universiteitsplein 1
> B-2610 Antwerp, Belgium
> http://webhost.ua.ac.be/deco
> tel : 32 3 820 23 85
>
>
>
> -----Original Message-----
> From: Jarrod Hadfield [mailto:jhadfiel at staffmail.ed.ac.uk]
> Sent: Thu 23-4-2009 9:33
> To: Strubbe Diederik
> Cc: r-sig-mixed-models at r-project.org; boristmr at yahoo.com
> Subject: Re: [R-sig-ME] About glmer with a quasi-Poisson
>
> Hi,
>
> This model can be fitted using MCMCglmm (and others I expect) if it is
> coded correctly. As Doug mentioned in an earlier email, calling
> different points in different sites the same name is not really
> necessary any more. If they are called different names then:
>
> random=~site+point
>
> is the same as the glmer model you use. If they are called the same
> thing then
>
> random=~site+site:point
>
> is the same as the glmer model below, as can be seen in the glmer summary.
>
> Coded this way many packages can probably fit the type of model you envisage.
>
> Cheers,
>
> Jarrod
>
>
>
>
>
> Quoting Strubbe Diederik <diederik.strubbe at ua.ac.be>:
>
>> Hey Mike,
>>
>> I am working on a similar dataset as you. I have point counts of
>> birds in 26 sites, in total 57 points [number of point counts per
>> site varying from 1 to 5]. These points were visited 5 times ( in
>> the same year). Most of the counts are 0 (> 90%), and the dataset is
>> thus zero-inflated - which, as far as I know is a special case of
>> overdispersion. There are 3 continuous explanatory variables under
>> consideration. I have been looking at different packages to fit a
>> model with
>> 1)a zero-inflated poison distribution
>> 2)a nested random effect to account for the fact that I have points
>> located within sites
>> 3)a random effect to account for the 5 visits (repeated measures)
>>
>> However, none of the packages that I tried is able to do this ( I
>> tried zeroinfl {pscl}, fmr {gnlm}, MCMCglmm {MCMCglmm} and
>> glmm.ADMB{glmmADMB}). The problem seems to be the nested random
>> effect (1|site/point in my case).
>>
>> The only function that I found able to handle this nested random
>> design and a distribution that comes close to be 'good' for my data
>> is glmer with a quasi-poisson distribution. My model syntax is:
>>
>> test <-
>> glmer(data$abundance~data$pca1+data$pca2+data$Gynoxis+(1|site/point)
>> + (1|visit),family=quasipoisson(link="log"))
>>
>> and yields the folowing output
>>
>> Generalized linear mixed model fit by the Laplace approximation
>> Formula: data$abundance ~ data$pca1 + data$pca2 + data$Gynoxis + (1
>> | site/point) + (1 | visit)
>> AIC BIC logLik deviance
>> 158.1 187.1 -71.05 142.1
>> Random effects:
>> Groups Name Variance Std.Dev.
>> point:site (Intercept) 0.0468101 0.216356
>> site (Intercept) 0.0000000 0.000000
>> visit (Intercept) 0.0012396 0.035208
>> Residual 0.0187386 0.136889
>> Number of obs: 276, groups: point:site, 57; site, 26; visit, 5
>>
>> Fixed effects:
>> Estimate Std. Error t value
>> (Intercept) -3.20340 0.10701 -29.936
>> data$pca1 -0.02629 0.03100 -0.848
>> data$pca2 -0.15055 0.02907 -5.180
>> data$Gynoxis -0.02246 0.01320 -1.702
>>
>> Correlation of Fixed Effects:
>> (Intr) dt$pc1 dt$pc2
>> data$pca1 -0.239
>> data$pca2 -0.098 -0.071
>> data$Gynoxs -0.878 0.330 0.205
>>
>> I was not aware of the comments of Ben Bolker on glmer and
>> quasi-Poisson, and am thus also interested in any responses on this.
>> Am I right to conclude that, ecept glmer, none of the packages
>> mentioned above is capable to handle nested random effects and
>> zero-inflated Poisson data?
>>
>> Cheers and many thanks,
>>
>> Diederik
>>
>>
>>
>> Diederik Strubbe
>> Evolutionary Ecology Group
>> Department of Biology, University of Antwerp
>> Universiteitsplein 1
>> B-2610 Antwerp, Belgium
>> http://webhost.ua.ac.be/deco
>> tel : 32 3 820 23 85
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>
>
>
> --
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
>
>
>
>
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list