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