[R-sig-ME] lme4 and calculating QAICc

[Ben Bolker]

  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