[R-sig-ME] Generalized Poisson Regression
Highland Statistics Ltd
highstat at highstat.com
Wed Jan 20 12:48:40 CET 2016
> ------------------------------
>
> Message: 2
> Date: Tue, 19 Jan 2016 15:32:24 -0500
> From: Ben Bolker <bbolker at gmail.com>
> To: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] Generalized Poisson Regression
> (Underdispersion)
> Message-ID: <569E9D58.2030104 at gmail.com>
> Content-Type: text/plain; charset=windows-1252; format=flowed
>
>
> Depends whether you mean in a standard (extended) GLM with only fixed
> effects (in which case the question is a little bit off-topic here ...)
> or in a GLMM.
>
> In the GLM case you can use the quasipoisson distribution in glm(); I
> think you have lots of other choices (including the ComPoissonReg
> package). library("sos"); ?findFn is your friend.
>
> In the GLMM case your options (that I know of) are a lot more
> limited. You could use MASS::glmmPQL with a quasi-Poisson response (but
> I would definitely test it on simulated data with known characteristic
> first!). You could use a GLMM with a Poisson response and adjust the
> standard errors yourself. Other distributions could be built into any
> of the 'toolbox' packages (glmmADMB, glmmTMB, brms), but you'd have to
> do it yourself or convince someone else ...
> You might get somewhere with the hglm (hierarchical GLM) package as
> well, although I don't immediately see an underdispersed Poisson family.
> Ordinal response models (e.g. the ordinal package) are another
> alternative ...
>
> Looking forward to hearing other people's answers
There is a bit of that in Ntzoufras (2010). He explains how you can do
the GP in MCMC. But unfortunately he only explains it for dealing with
overdispersion (and the majority of papers on this topic do the same).
It is easy to program in JAGS with the zero trick. And quite often it
performs better than the negative binomial GLM. This code should work
for a GP model to deal with overdispersion.
model{
#1. Priors beta
for (i in 1:K) { beta[i] ~ dnorm(0, 0.0001)}
#1C. Prior for theta parameter of GP distribution
theta ~ dunif(0, 1)
#2. Likelihood using the zero trick
C <- 10000
for (i in 1:N){
Zeros[i] ~ dpois(Zeros.mean[i])
Zeros.mean[i] <- -L[i] + C
l1[i] <- log(1 - theta)
l2[i] <- log(lambda[i])
l3[i] <- (Y[i] - 1) * log((1-theta) * lambda[i] + theta * Y[i])
l4[i] <- - ((1 - theta) * lambda[i] + theta * Y[i])
l5[i] <- -loggam(Y[i] + 1)
L[i] <- l1[i] + l2[i] + l3[i] + l4[i] + l5[i]
eta[i] <- inprod(beta[], X[i,])
lambda[i] <- exp(eta[i])
}
#Mean and variance
for (i in 1:N){
ExpY[i] <- lambda[i]
VarY[i] <- lambda[i] / ((1 - theta)^2)
Pres[i] <- (Y[i] - ExpY[i]) / sqrt(VarY[i])
}
AIC <- -2 * sum(L[1:N]) + 2 * (K + 1)
}
For underdispersion it gets more ugly. The theta can be between -1 and
0. But you need to include a restriction when sampling new parameters
within the MCMC algorithm. Can't remember exactly what it
was...something with the sum of two squares and being positive (See:
Modeling Underdispersed Count Data with Generalized Poisson Regression.
Harris et. (2011)).
In normal words: You need to write your own MCMC algorithm...... Maybe
someone has done this already? But once you have it running, you can
easily add random effects, and other correlation structures.
You can use rLGP from RMKdiscrete to simulate GP data. Gives you an
impression how data from such a distribution looks like.
There is a whole book about the GP distribution....by Consul (1989).
By the way..it is also called the log linear Lagrangian Poisson model.
Just in case you want to impress your colleagues.
Alain
--
Dr. Alain F. Zuur
First author of:
1. Beginner's Guide to GAMM with R (2014).
2. Beginner's Guide to GLM and GLMM with R (2013).
3. Beginner's Guide to GAM with R (2012).
4. Zero Inflated Models and GLMM with R (2012).
5. A Beginner's Guide to R (2009).
6. Mixed effects models and extensions in ecology with R (2009).
7. Analysing Ecological Data (2007).
Highland Statistics Ltd.
9 St Clair Wynd
UK - AB41 6DZ Newburgh
Tel: 0044 1358 788177
Email: highstat at highstat.com
URL: www.highstat.com
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list