[R-sig-ME] Generalized Poisson Regression

[Highland Statistics Ltd]

> ------------------------------
>
> 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]]