Hi Philip,

Thanks very much for this, i was completely unaware. I  have read various papers using lmer to calculate the AIC statistic and none have mentioned this?

I have just run through a random section of my models with this correction, however the AIC / BIC values are the same with the REML=F in and out?

Chris

MODEL WITH REML=F

> model0 <- lmer(threatornot~1+(1|order/family) + seasonality + pollendispersal +  breedingsystem*fruit + habit + lifeform +  woodyness, family=binomial,REML=F )


Generalized linear mixed model fit by the Laplace approximation 
Formula: threatornot ~ 1 + (1 | order/family) + seasonality + pollendispersal +      breedingsystem * fruit + habit + lifeform + woodyness 
  AIC  BIC logLik deviance
 1399 1479 -683.6     1367
Random effects:
 Groups       Name        Variance Std.Dev.
 family:order (Intercept) 0.27526  0.52466 
 order        (Intercept) 0.00000  0.00000 
Number of obs: 1116, groups: family:order, 43; order, 9

Fixed effects:
                        Estimate Std. Error z value Pr(>|z|)   
(Intercept)             0.384574   0.734960   0.523  0.60079   
seasonality2           -1.127996   0.353013  -3.195  0.00140 **
pollendispersal2        0.693255   0.314600   2.204  0.02755 * 
breedingsystem2         0.761067   0.493404   1.542  0.12296   
breedingsystem3         1.226269   0.557236   2.201  0.02776 * 
fruit2                  1.047648   0.616723   1.699  0.08937 . 
habit2                 -1.146334   0.551682  -2.078  0.03772 * 
habit3                 -0.731207   0.872805  -0.838  0.40216   
habit4                 -0.190900   0.551427  -0.346  0.72920   
lifeform2              -0.295342   0.182667  -1.617  0.10592   
lifeform3              -0.376204   0.501825  -0.750  0.45345   
woodyness2              0.006274   0.390241   0.016  0.98717   
breedingsystem2:fruit2 -1.273811   0.651011  -1.957  0.05039 . 
breedingsystem3:fruit2 -1.633424   0.744563  -2.194  0.02825 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

MODEL WITHOUT REML=F

model126 <- lmer(threatornot~1+(1|order/family) + seasonality + pollendispersal +  breedingsystem*fruit + habit + lifeform +  woodyness, family=binomial)

Generalized linear mixed model fit by the Laplace approximation 
Formula: threatornot ~ 1 + (1 | order/family) + seasonality + pollendispersal +      breedingsystem * fruit + habit + lifeform + woodyness 
  AIC  BIC logLik deviance
 1399 1479 -683.6     1367
Random effects:
 Groups       Name        Variance Std.Dev.
 family:order (Intercept) 0.27526  0.52466 
 order        (Intercept) 0.00000  0.00000 
Number of obs: 1116, groups: family:order, 43; order, 9

Fixed effects:
                        Estimate Std. Error z value Pr(>|z|)   
(Intercept)             0.384574   0.734960   0.523  0.60079   
seasonality2           -1.127996   0.353013  -3.195  0.00140 **
pollendispersal2        0.693255   0.314600   2.204  0.02755 * 
breedingsystem2         0.761067   0.493404   1.542  0.12296   
breedingsystem3         1.226269   0.557236   2.201  0.02776 * 
fruit2                  1.047648   0.616723   1.699  0.08937 . 
habit2                 -1.146334   0.551682  -2.078  0.03772 * 
habit3                 -0.731207   0.872805  -0.838  0.40216   
habit4                 -0.190900   0.551427  -0.346  0.72920   
lifeform2              -0.295342   0.182667  -1.617  0.10592   
lifeform3              -0.376204   0.501825  -0.750  0.45345   
woodyness2              0.006274   0.390241   0.016  0.98717   
breedingsystem2:fruit2 -1.273811   0.651011  -1.957  0.05039 . 
breedingsystem3:fruit2 -1.633424   0.744563  -2.194  0.02825 * 

On 5 Aug 2010, at 12:08, Philip Dixon wrote:

I believe the recent discussion about AIC or p-values has missed a crucial practical issue.

The AIC statistic reported by a default call to lmer() has NOTHING to do with the choice of fixed effects.  lmer() uses reml to define the fit of a model.  REML (residual ML or restricted ML) is a log likelihood for a variance-covariance structure based on the residuals after fitting the fixed effects.  This has two consequences:
1) the AIC calculated from the reml lnL only informs you about the fit of the random effects model.
2) the reml AIC can only be compared between models with the same fixed effects.  Otherwise, the lnL is being calculated from different data (because different residuals with different X's).

This is a very common issue/problem.

use REML=F to force lmer() to calculate the full ML lnL.

Best wishes,
Philip Dixon

_______________________________________________
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology


	[[alternative HTML version deleted]]

