[R-sig-ME] Subject: Re: About glmer with a quasi-Poisson
Jarrod Hadfield
j.hadfield at ed.ac.uk
Sun Apr 26 22:05:05 CEST 2009
Hi Nicholas,
Specifying rcov=~idh(trait):units does as you suggest, with the
covariance explicitly set to zero. 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.
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