[R] Generalised poisson regression

[David Winsemius]

On Aug 22, 2015, at 8:40 AM, Babatunde Yakub wrote:

> Using object$deviance for generalised poisson regression model gives NULL as response. Code attached. please check....thanks

I'm not attacking your manhood when I say that your attachments are a problem. The only person who can see that attachment besides yourself is me because you did not copy the list and used an extension of `.r`. The mailing list requires that attachments be MIME-text and most mailers will not construct an email that lists them as such unless you use an extension of .txt. So your attachment with an .r extension was scrubbed. I'm not a user of the VGAM package, but I do read help pages. First your code which will let the rest of the audience correct my errors:

#EQUI-DISPERSION POISSON COUNT
library(MASS)
library(GPseq)
library(VGAM)
library(COMPoissonReg)
x=rnorm(50)
link=0.2+0.4*x
rpois.od<-function (n,lambda,d) {
   if (d==1)
      rpois(n, lambda)
   else
      rnbinom(n, size=(lambda/(d-1)), mu=lambda)
}
y=rpois.od(50,lambda=exp(link),d=1)

# GENERALIZED POISSON
genp= vglm(y~x, genpoisson,trace = TRUE)
gen=summary(genp)


The need for 4 different packages was no at all clear. I installed pkg:VGAM and looked at `?genpoisson` after loading only that package. Since it was there then only pkg:COMPoissonReg might be modifying the behavior of htat code, so I left it uninstalled (and unloaded initially). The ?vglm page tells us that there is no deviance component to objects returned from it and states that the recommended method for doing LR-tests is to use `lrtest` and so my w

The 'gen' object was examined with `str`:

str(gen)

#--------output truncated at beginning and end
 ..@ family          :Formal class 'vglmff' [package "VGAM"] with 18 slots
 .. .. ..@ blurb             : chr [1:8] "Generalized Poisson distribution\n\n" "Links:    " "rhobit(lambda)" ", " ...
 .. .. ..@ constraints       :  expression({     M1 <- 2     dotzero <- -1     eval(negzero.expression.VGAM) })
 .. .. ..@ deviance          :function ()  
 .. .. ..@ fini              :  expression({ })
 .. .. ..@ first             :  expression({ })
 .. .. ..@ infos             :function (...)  

#-------

So that didn't look very promising, ..... seeing that empty `deviance` function. So I continued with my hacking efforts at following the manual and built a NULL model and ran an LRT:

?lrtest

> genpnull= vglm(y~1, genpoisson,trace = TRUE)
VGLM    linear loop  1 :  loglikelihood = -64.750744
VGLM    linear loop  2 :  loglikelihood = -64.696782
VGLM    linear loop  3 :  loglikelihood = -64.696568
VGLM    linear loop  4 :  loglikelihood = -64.696567
VGLM    linear loop  5 :  loglikelihood = -64.696567
> gen0=summary(genpnull)
> lrtest(genp,genpnull)
Likelihood ratio test

Model 1: y ~ x
Model 2: y ~ 1
 #Df  LogLik Df  Chisq Pr(>Chisq)    
1  97 -59.045                         
2  98 -64.697  1 11.304  0.0007735 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


The LogLik values would allow you to construct a "deviance"-value for each model although the `lrtest` function has already calculated the difference in 2*LL between the 2 models and labeled it 'Chisq'. The deviance is not of much meaning as a number applied to a single model because it can vary drastically for different data arrangements. Only differences in nested models is meaningful.

Turned out that the actions I just performed are already done automatically if a single model is given to lrtest. Exactly the same output results. You can look at the code:

lrtest_vglm

# and see where it does that
----snippet----
   if (nmodels < 2) {
       objects <- c(objects, . ~ 1)
       nmodels <- 2
   }
#-----and then see that it is using a function named LogLik------

logLlist <- lapply(objects, logLik)

#-----------

> logLik(genp)
[1] -59.04466

So if I had really known what I was doing, I might have just typed ?logLik, except when I do that it appears the package author doesn't really intend for us to do this, since that pulls up a page for "Undocumented and Internally Used Functions and Classes". Thus endeth today's meanderings among the help pages for pkg:VGAM.

Please read the Posting Guide and also go back and read the general information webpages. You seem to have missed some important details.

-- 
David.

> 
> 
> On Friday, August 21, 2015 7:40 PM, David Winsemius <dwinsemius at comcast.net> wrote:
> 
> 
> 
> On Aug 21, 2015, at 5:22 AM, Babatunde Yakub via R-help wrote:
> 
>> I want to know how to extract or obtain the deviance for a fitted generalised poisson regression model. Thanks in advance
> 
> 
> If you post the code and some sample data, for building such a model, I'm sure someone can help you extract the deviance. If you simply mean what is returned by an ordinary glm-call with family="poisson" then deviance should be one of the elements of the glm-object.
> 
> 
> object$deviance
> 
>>    [[alternative HTML version deleted]]
> 
> When you do reply (if needed) please send in plain text rather than HTML.
> 
> -- 
> David Winsemius
> Alameda, CA, USA
> 
> 
> 
> <for mail.r>

David Winsemius
Alameda, CA, USA