[R-sig-ME] Subject: Re: About glmer with a quasi-Poisson
Nicholas Lewin-Koh
nikko at hailmail.net
Sun Apr 26 21:41:10 CEST 2009
Hi Jarrod,
Since there is no information in the data for the covariance between the
Poisson process and the
zero-inflation process, in essence the model is not identified. If I am
understanding the structure
correctly, isn't it quite dangerous to even put a prior on the
covariance? Wouldn't it be better to
write the model in a way that the covariance structure is explicit? For
instance could one
not assume that the Poisson process and the zero inflation are
separable, can be written as
a product of the two distributions? I am shooting from the hip here, as
usual, but maybe I just need to write out the model to see it.
Nicholas
> Message: 3
> Date: Sat, 25 Apr 2009 19:45:58 +0100
> From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
> Subject: Re: [R-sig-ME] About glmer with a quasi-Poisson
> To: Strubbe Diederik <diederik.strubbe at ua.ac.be>
> Cc: r-sig-mixed-models at r-project.org, boristmr at yahoo.com
> Message-ID: <20090425194558.ri00kmkf4k0o0c4k at www.staffmail.ed.ac.uk>
> Content-Type: text/plain; charset=ISO-8859-1; DelSp="Yes";
> format="flowed"
>
> 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.
>
>
>
> ------------------------------
>
> _______________________________________________
> R-sig-mixed-models mailing list
> R-sig-mixed-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
> End of R-sig-mixed-models Digest, Vol 28, Issue 29
> **************************************************
More information about the R-sig-mixed-models
mailing list