[R-sig-ME] AIC comparisons between lmer and glm

[Ben Bolker]

   Yes, but: if treating the random effects as a nuisance parameter or
making inferences at the population level (see Vaida and Blanchard
2005), treating the random effect as 1 parameter is probably
approximately OK, although it is likely to be conservative because of
boundary issues (Greven et al 2008) ...

  Some proof that the deviance calculated with variance set equal to
zero in the glmer fit is not the same as the deviance from the glm --
don't know if this is a mismatch in which constants etc. are included,
or something about the limit of the penalized likelihood ...


### function to evaluate likelihood etc. with a new value of the
### random-effects variance
### DANGER DANGER DANGER DANGER
###  at least with some versions of lme4, this code will modify
### the original model.

update_dev <- function(mm,sd) {
  .Call("mer_ST_setPars", mm, sd, PACKAGE = "lme4")
  .Call("mer_update_L", mm, PACKAGE = "lme4")
  res <- try(.Call("mer_update_RX", mm, PACKAGE = "lme4"), silent = TRUE)
  if (inherits(res, "try-error")) {
    val <- NA
  } else {
    .Call("mer_update_ranef", mm, PACKAGE = "lme4")
    .Call("mer_update_dev", mm, PACKAGE = "lme4") ## added for glmmML
    val <- mm at deviance
  }
  val
}

library(lme4)
(gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
              family = binomial, data = cbpp))

gm1 at deviance["ML"] ## deviance from original model (100.0959)


dd = update_dev(gm1,sd=0)
dd["ML"]  ## deviance with random effects variance set to zero
## (114.9793)
(gm0 <- glm(cbind(incidence, size - incidence) ~ period,
              family = binomial, data = cbpp))

deviance(gm0) ## deviance from glm fit
## (114.1017)


Peter Dixon wrote:
> I think one of the issues, discussed elsewhere on this list if I  
> recall, is that the status of the random effects estimates is unclear  
> – they are not parameters of the model in the same sense as the fixed  
> effects are. Because AIC depends on counting parameters, comparing AIC  
> values with different random-effects structures is suspect.
> 
> - Pete
> 
> On Jul 29, 2009, at 10:11 AM, cacabelos at uvigo.es wrote:
> 
>> Hi R-users,
>>
>> I was unfortunately not able to find a solution to my problem about  
>> model selection. I have fixed (e.g. Height) and random (e.g.Time)  
>> factors, and I created these kind of model structures:
>>
>> Model1<-lmer(H~Height+(1|Time)+Biomass)+(Biomass| 
>> Time),data=dat,family=gaussian)
>> Model2<-glm(H~Height+Biomass,data=dat,family=gaussian)
>>
>> I am using the criterion ?the best model is the one that has the  
>> lowest AIC?, but are these AIC from different procedures (i.e. glm  
>> and lmer) comparable?
>>
>> I wonder if anyone can help me... Thank you!
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
> 
> 
> 
> ---
> Peter Dixon
> peter.dixon at ualberta.ca
> 
> 
> 
> 
> 	[[alternative HTML version deleted]]
> 
> 


-- 
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / www.zoology.ufl.edu/bolker
GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc