Dear R users,
I am having problems trying to fit quasipoisson and negative binomials glm.
My data set
contains abundance (counts) of a species under different management regimens.
First, I tried to fit a poisson glm:
> summary(model.p<-glm(abund~mgmtcat,poisson))
Call:
glm(formula = abund ~ mgmtcat, family = poisson)
.
.
.
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 1904.7 on 19 degrees of freedom
Residual deviance: 1154.3 on 16 degrees of freedom
AIC: 1275.4
Number of Fisher Scoring iterations: 4
Wich suggests the existence of STRONG overdispersion, so I tried:
> summary(model.qp<-glm(abund~mgmtcat,quasipoisson))
Call:
glm(formula = abund ~ mgmtcat, family = quasipoisson)
.
.
.
(Dispersion parameter for quasipoisson family taken to be 73.51596)
Null deviance: 1904.7 on 19 degrees of freedom
Residual deviance: 1154.3 on 16 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 4
Here I found the first problem: AIC is not available.
I know that count data for the studied species usually show aggregation.
So, I fitted
a negative binomial glm with the glm.nb in MASS:
> summary.negbin(model.nb<-glm.nb(abund~mgmtcat))
Call: glm.nb(formula = abund ~ mgmtcat, init.theta =
1.23560100958978, link = log)
.
.
.
(Dispersion parameter for Negative Binomial(1.2356) family taken to
be 1)
Null deviance: 33.173 on 19 degrees of freedom
Residual deviance: 22.316 on 16 degrees of freedom
AIC: -15948
Number of Fisher Scoring iterations: 1
Correlation of Coefficients:
(Intercept) mgmtcat1 mgmtcat2
mgmtcat1 -0.7052
mgmtcat2 -0.7053 0.4974
mgmtcat3 -0.7005 0.4940 0.494
Theta: 1.236
Std. Err.: 0.362
2 x log-likelihood: -211.079
And now, I am getting a negative AIC value! I have seen that this problem
have been discused in the S-news list.
Much of the discussion there is far beyond my statistical and R knowledge.
One of the solutions proposed there
was adding - lgamma(y +1) to the internal function loglik in glm.nb, but I
have seen that the current version of
MASS contains that term.
My problem is that I want to compare the quasipoisson and negative binomial
models, and I have a NA value and a negative one.
Can I obtain an AIC for the quasipoisson model? What about the negative
AIC? Can I use it or do you think that anything is wrong?
Thanks in advance,
--
Vicente Piorno
Departamento de Ecologia y Biologia Animal - Universidad de Vigo
EUIT Forestal - Campus Universitario
36005 Pontevedra SPAIN
Vicente Piorno wrote:> > Dear R users, > I am having problems trying to fit quasipoisson and negative binomials glm. > My data set > contains abundance (counts) of a species under different management regimens. > First, I tried to fit a poisson glm: > > > summary(model.p<-glm(abund~mgmtcat,poisson)) > > Call: > glm(formula = abund ~ mgmtcat, family = poisson) > . > . > . > (Dispersion parameter for poisson family taken to be 1) > > Null deviance: 1904.7 on 19 degrees of freedom > Residual deviance: 1154.3 on 16 degrees of freedom > AIC: 1275.4 > Number of Fisher Scoring iterations: 4 > > Wich suggests the existence of STRONG overdispersion, so I tried: > > > summary(model.qp<-glm(abund~mgmtcat,quasipoisson)) > > Call: > glm(formula = abund ~ mgmtcat, family = quasipoisson) > . > . > . > (Dispersion parameter for quasipoisson family taken to be 73.51596) > > Null deviance: 1904.7 on 19 degrees of freedom > Residual deviance: 1154.3 on 16 degrees of freedom > AIC: NA > Number of Fisher Scoring iterations: 4 > > Here I found the first problem: AIC is not available. > > I know that count data for the studied species usually show aggregation. > So, I fitted > a negative binomial glm with the glm.nb in MASS: > > > summary.negbin(model.nb<-glm.nb(abund~mgmtcat)) > > Call: glm.nb(formula = abund ~ mgmtcat, init.theta > 1.23560100958978, link = log) > . > . > . > (Dispersion parameter for Negative Binomial(1.2356) family taken to > be 1) > > Null deviance: 33.173 on 19 degrees of freedom > Residual deviance: 22.316 on 16 degrees of freedom > AIC: -15948 > Number of Fisher Scoring iterations: 1 > > Correlation of Coefficients: > (Intercept) mgmtcat1 mgmtcat2 > mgmtcat1 -0.7052 > mgmtcat2 -0.7053 0.4974 > mgmtcat3 -0.7005 0.4940 0.494 > > Theta: 1.236 > Std. Err.: 0.362 > 2 x log-likelihood: -211.079 > > And now, I am getting a negative AIC value! I have seen that this problem > have been discused in the S-news list. > Much of the discussion there is far beyond my statistical and R knowledge. > One of the solutions proposed there > was adding - lgamma(y +1) to the internal function loglik in glm.nb, but I > have seen that the current version of > MASS contains that term. > > My problem is that I want to compare the quasipoisson and negative binomial > models, and I have a NA value and a negative one. > Can I obtain an AIC for the quasipoisson model? What about the negative > AIC? Can I use it or do you think that anything is wrong? > > Thanks in advance, > > -- > Vicente Piorno > Departamento de Ecologia y Biologia Animal - Universidad de Vigo > EUIT Forestal - Campus Universitario > 36005 Pontevedra SPAIN > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-helpWhat's wrong with a negative AIC? I would vote against comparing different model classes using AIC. Instead, it's a better idea to think about which model class makes sense in a given context. Uwe Ligges
1) A quasi- mode does not have a likelihood and so does not have an AIC,
by definition.
2) There is nothing wrong with a negative AIC, but you don't have one.
The value of AIC in your problem is 211.079 + 2*16 + 2. It is summary.glm
that is giving you misleading results. Try the following
logLik.negbin <- function(object, ...)
{
if (length(list(...)))
warning("extra arguments discarded")
p <- object$rank + 1 # for theta
val <- object$twologlik/2
attr(val, "df") <- p
class(val) <- "logLik"
val
}
then AIC(fit) will work.
3) Don't call methods like summary.negbin directly: they will not be
user-visible.
On Wed, 12 Mar 2003, Vicente Piorno wrote:
> Dear R users,
> I am having problems trying to fit quasipoisson and negative binomials glm.
> My data set
> contains abundance (counts) of a species under different management
regimens.
> First, I tried to fit a poisson glm:
>
> > summary(model.p<-glm(abund~mgmtcat,poisson))
>
> Call:
> glm(formula = abund ~ mgmtcat, family = poisson)
> .
> .
> .
> (Dispersion parameter for poisson family taken to be 1)
>
> Null deviance: 1904.7 on 19 degrees of freedom
> Residual deviance: 1154.3 on 16 degrees of freedom
> AIC: 1275.4
> Number of Fisher Scoring iterations: 4
>
> Wich suggests the existence of STRONG overdispersion, so I tried:
>
> > summary(model.qp<-glm(abund~mgmtcat,quasipoisson))
>
> Call:
> glm(formula = abund ~ mgmtcat, family = quasipoisson)
> .
> .
> .
> (Dispersion parameter for quasipoisson family taken to be 73.51596)
>
> Null deviance: 1904.7 on 19 degrees of freedom
> Residual deviance: 1154.3 on 16 degrees of freedom
> AIC: NA
> Number of Fisher Scoring iterations: 4
>
> Here I found the first problem: AIC is not available.
>
> I know that count data for the studied species usually show aggregation.
> So, I fitted
> a negative binomial glm with the glm.nb in MASS:
>
> > summary.negbin(model.nb<-glm.nb(abund~mgmtcat))
>
> Call: glm.nb(formula = abund ~ mgmtcat, init.theta =
> 1.23560100958978, link = log)
> .
> .
> .
> (Dispersion parameter for Negative Binomial(1.2356) family taken to
> be 1)
>
> Null deviance: 33.173 on 19 degrees of freedom
> Residual deviance: 22.316 on 16 degrees of freedom
> AIC: -15948
> Number of Fisher Scoring iterations: 1
>
> Correlation of Coefficients:
> (Intercept) mgmtcat1 mgmtcat2
> mgmtcat1 -0.7052
> mgmtcat2 -0.7053 0.4974
> mgmtcat3 -0.7005 0.4940 0.494
>
> Theta: 1.236
> Std. Err.: 0.362
> 2 x log-likelihood: -211.079
>
> And now, I am getting a negative AIC value! I have seen that this problem
> have been discused in the S-news list.
> Much of the discussion there is far beyond my statistical and R knowledge.
> One of the solutions proposed there
> was adding - lgamma(y +1) to the internal function loglik in glm.nb, but I
> have seen that the current version of
> MASS contains that term.
>
> My problem is that I want to compare the quasipoisson and negative binomial
> models, and I have a NA value and a negative one.
> Can I obtain an AIC for the quasipoisson model? What about the negative
> AIC? Can I use it or do you think that anything is wrong?
>
> Thanks in advance,
>
>
> --
> Vicente Piorno
> Departamento de Ecologia y Biologia Animal - Universidad de Vigo
> EUIT Forestal - Campus Universitario
> 36005 Pontevedra SPAIN
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595