Dear all,
I am trying to calculate QAICc using to compare two Poisson models. Unfortunately all I seem to get as printed values is NaN.
Is there something I'm missing? Even though I am able to generate model output, I do receive "convergence errors". Would these warning message have anything to do with this?
Incidentally, I'm using the methodology extracted from Bolker et al (2009)
Here is the code I have used.
######
library(lme4)
######
mp1=lmer(abundance~year+controlA+(year|groupC:sitecode),family="poisson",data=testData)
######
mq1=lmer(abundance~year+controlA+(year|groupC:sitecode),family="quasipoisson",data=testData)
######
QAICc <- function(mod, scale, QAICc = TRUE) {
LL <- logLik(mod)
ll <- as.numeric(LL)
df <- attr(LL, "df")
n <- length(mod at y)
if (QAICc)
qaic = as.numeric(-2 * ll/scale + 2 * df + 2 * df * (df +
1)/(n - df - 1))
else qaic = as.numeric(-2 * ll/scale + 2 * df)
qaic
}
#######
QAICc(mq1,scale=phi)
....and this is what I generate...
[1] NaN
Thank you.
Andrew
Andrew Close
Research Associate
Institute for Research on Environment and Sustainability (IRES)
School of Biology
4th Floor Devonshire Building
Newcastle University
NE1 7RU
+44 (0)191 2464840
lme4 and calculating QAICc
4 messages · Andrew Close, Christoph Scherber, David Atkins +1 more
Dear Andrew, You might want to check if you can extract a log-Likelihood from your models logLik(mq1) If you get an NaN here, there may be something wrong with your model(s). Best wishes Christoph
Dear all,
I am trying to calculate QAICc using to compare two Poisson models.
Unfortunately all I seem to get as printed values is NaN.
Is there something I'm missing? Even though I am able to generate model
output, I do receive "convergence errors". Would these warning message
have anything to do with this?
Incidentally, I'm using the methodology extracted from Bolker et al (2009)
Here is the code I have used.
######
library(lme4)
######
mp1=lmer(abundance~year+controlA+(year|groupC:sitecode),family="poisson",data=testData)
######
mq1=lmer(abundance~year+controlA+(year|groupC:sitecode),family="quasipoisson",data=testData)
######
QAICc <- function(mod, scale, QAICc = TRUE) {
LL <- logLik(mod)
ll <- as.numeric(LL)
df <- attr(LL, "df")
n <- length(mod at y)
if (QAICc)
qaic = as.numeric(-2 * ll/scale + 2 * df + 2 * df * (df +
1)/(n - df - 1))
else qaic = as.numeric(-2 * ll/scale + 2 * df)
qaic
}
#######
QAICc(mq1,scale=phi)
....and this is what I generate...
[1] NaN
Thank you.
Andrew
Andrew Close
Research Associate
Institute for Research on Environment and Sustainability (IRES)
School of Biology
4th Floor Devonshire Building
Newcastle University
NE1 7RU
+44 (0)191 2464840
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
Andrew-- Your code worked fine for me using an lmer() fit on some of my own data, and setting scale = 1. Is phi defined somewhere that we can't see, or might that be causing your problems? I think it should work if scale is passed a numeric value. cheers, Dave
Dave Atkins, PhD
Research Associate Professor
Center for the Study of Health and Risk Behaviors
Department of Psychiatry and Behavioral Science
University of Washington
1100 NE 45th Street, Suite 300
Seattle, WA 98105
206-616-3879
datkins at u.washington.edu
Dear all,
I am trying to calculate QAICc using to compare two Poisson models.
Unfortunately all I seem to get as printed values is NaN.
Is there something I'm missing? Even though I am able to generate model
output, I do receive "convergence errors". Would these warning message
have anything to do with this?
Incidentally, I'm using the methodology extracted from Bolker et al (2009)
Here is the code I have used.
######
library(lme4)
######
mp1=lmer(abundance~year+controlA+(year|groupC:sitecode),family="poisson",data=testData)
######
mq1=lmer(abundance~year+controlA+(year|groupC:sitecode),family="quasipoisson",data=testData)
######
QAICc <- function(mod, scale, QAICc = TRUE) {
LL <- logLik(mod)
ll <- as.numeric(LL)
df <- attr(LL, "df")
n <- length(mod at y)
if (QAICc)
qaic = as.numeric(-2 * ll/scale + 2 * df + 2 * df * (df +
1)/(n - df - 1))
else qaic = as.numeric(-2 * ll/scale + 2 * df)
qaic
}
#######
QAICc(mq1,scale=phi)
....and this is what I generate...
[1] NaN
Thank you.
Andrew
Andrew Close
Research Associate
Institute for Research on Environment and Sustainability (IRES)
School of Biology
4th Floor Devonshire Building
Newcastle University
NE1 7RU
+44 (0)191 2464840
Checking in a bit late here ... been busy the last couple of days.
Dr. Christoph Scherber wrote:
Dear Andrew, You might want to check if you can extract a log-Likelihood from your models logLik(mq1) If you get an NaN here, there may be something wrong with your model(s). Best wishes Christoph
Dear all, I am trying to calculate QAICc using to compare two Poisson models. Unfortunately all I seem to get as printed values is NaN. Is there something I'm missing? Even though I am able to generate model output, I do receive "convergence errors". Would these warning message have anything to do with this? Incidentally, I'm using the methodology extracted from Bolker et al (2009) Here is the code I have used. ###### library(lme4) ###### mp1=lmer(abundance~year+controlA+(year|groupC:sitecode),family="poisson",data=testData) ###### mq1=lmer(abundance~year+controlA+(year|groupC:sitecode),family="quasipoisson",data=testData)
OR mq1 = update(mp1,family="quasipoisson") (as in Bolker et al supplement)
######
QAICc <- function(mod, scale, QAICc = TRUE) {
LL <- logLik(mod)
ll <- as.numeric(LL)
df <- attr(LL, "df")
n <- length(mod at y)
if (QAICc)
qaic = as.numeric(-2 * ll/scale + 2 * df + 2 * df * (df +
1)/(n - df - 1))
else qaic = as.numeric(-2 * ll/scale + 2 * df)
qaic
}
#######
QAICc(mq1,scale=phi)
In order to use QAICc you have to do the following:
phi = lme4:::sigma(mq1)
QAICc(mp1,scale=phi)
(see p. 8 of the Bolker et al supplement)
The basic problem is that quasi- models don't return likelihoods,
and non-quasi- models don't return estimates of scale parameters,
so you have fit both and combine the information.
good luck,
Ben Bolker