[R-sig-ME] Subject: Re: About glmer with a quasi-Poisson
Nicholas Lewin-Koh
nikko at hailmail.net
Mon Apr 27 06:28:39 CEST 2009
Hi Jarrod,
Comments below
On Sun, 26 Apr 2009 21:05 +0100, "Jarrod Hadfield" <j.hadfield at ed.ac.uk>
wrote:
> Hi Nicholas,
>
> Specifying rcov=~idh(trait):units does as you suggest, with the
> covariance explicitly set to zero.
Ah, I wasn't sure what the idh() function did
> I'm not sure there's technically
> anything wrong with specifying a fully parameterised covariance matrix
> (us(trait):units) with the covariance information coming from the
> prior, as long as it's recognised as such.
What worries me here is that the prior on the covariance can have
effects on the other model parameters, and since
the unspecified covariance model is not identified, ie no information
from the data, weird things can happen. Clearly one would not pay
attention to the covariance in the posterior, as it is just the prior,
but what is that doing to everything else? Maybe it is just me, but that
would
make me very uncomfortable about trusting anything else in the model. At
least explicitly setting the covariance, you can understand what is
going on.
Nicholas
>
> Cheers,
>
> Jarrod
>
>
>
>
>
> Quoting Nicholas Lewin-Koh <nikko at hailmail.net>:
>
> > 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
> >> **************************************************
> >
> >
>
>
>
> --
> 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