lme4 and calculating QAICc
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