[R-sig-ME] lme4 and calculating QAICc

[Ben Bolker]

Marcus Rowcliffe wrote:
> Dear all
> I found the following useful exchange in the archive, but there are a
> couple of crucial things confusing me:
> 1.	In the Bolker et al supplement, it appears that his QAICc
> function (essentially the same as that outlined below) is run on quasi
> models, so the log-likelihood used is from quasi models, not straight
> poisson, as recommended below. In fact, looking at simple model fits,
> log-likelihoods for poisson and quasipoisson appear to be the same
> anyway, so why is it necessary to fit both types of model as suggested
> by Bolker below?

     Some confusion (on my part) between different kinds of models: for
more traditional generalized quasilikelihood models (e.g. fitted by
glm() or glmmPQL()), R refuses (perhaps rightly) to return a
log-likelihood, hence one has to fit the non-quasi- version to get the
equivalent likelihood value.  Not really necessary for glmer, although I
can imagine this will change in the future ... ?

  **NB** someone (alas I can no longer remember who nor quickly find the
e-mail) pointed out that the "scale" estimate for QAICc should be
lme4:::sigma(model)^2 [i.e. SQUARED].  I've been meaning to send in a
correction to Elsevier but haven't gotten around to it yet.

> 2.	Comparing poisson and quasipoisson models, I note that the
> former has no residual random variance, and hence one fewer degrees of
> freedom than the quasipoisson. Why is this? It obviously has important
> implications for the QAIC calculation.

  I find this a bit difficult myself.  Technically, it's another
'estimated parameter', but it's a nuisance parameter that does *not*
affect the predictions at all (it is simply calculated from the
residuals of the model), so I'm not really sure whether it should affect
the estimate of expected Kullback-Leibler distance or not ... however, I
would quibble a little bit with "important implications for the QAIC
calculation" -- it's just 1 degree of freedom!  To quote _Numerical
Recipes in C_ 2d ed. p. 611 (Press et al): "We might also comment that
if the difference between N and N-1 ever matters to you, then you are
probably up to no good anyway -- e.g., trying to substantiate a
questionable hypothesis with marginal data.")

> Apologies if the latter should be obvious and I've missed it somewhere.
> Marcus
> 
> 
> 
>>> 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="quasipoi
> sson",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
> 
> 
> The Zoological Society of London is incorporated by Royal Charter
> Principal Office England. Company Number RC000749
> Registered address: 
> Regent's Park, London, England NW1 4RY
> Registered Charity in England and Wales no. 208728 
> 
> _________________________________________________________________________
> This e-mail has been sent in confidence to the named a...{{dropped:13}}