[R-sig-ME] Fixed effects only model with lme4

[Ben Bolker]

 A quick hack to compute a likelihood profile, for changing
random-effects variances.  A slight special case is to compute the
likelihood of a model with variance of the random effect set to zero.
This works only for a single-random-effect model, and is pulled from
code in Doug Bates's vignette.  In the past this would have dangerously
modified your original model, but I think that's fixed now.

   Other glmm hacks are available at
http://glmm.wikidot.com/local--files/reef-fish/glmmfuns.R

  I would be happy to hear about extensions or corrections.

  cheers
    Ben Bolker

## extract likelihood based on zero variance for a single random
##  effect
zerodev <- function(mm) {
  varprof(mm,0,0,1)[["ML"]]
}

varprof <- function(mm,lower=0,upper=20,n=101) {
  sg <- seq(lower, upper, len = n)
  orig.sd <- attr(VarCorr(mm)[[1]],"stddev")
  dev <- mm at deviance
  nc <- length(dev)
  nms <- names(dev)
  vals <- matrix(0, nrow = length(sg), ncol = nc, dimnames = list(NULL,
nms))
  update_dev <- function(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
  }
  for (i in seq(along = sg)) {
      vals[i,] <- update_dev(sg[i])
    }
  update_dev(orig.sd) ## hack! to restore original sd
  data.frame(sd=sg,vals)
}


Emmanuel Charpentier wrote:
> Dear Pr Bates,
> 
> Le samedi 06 juin 2009 à 12:49 -0500, Douglas Bates a écrit :
>> On Sat, Jun 6, 2009 at 12:03 PM, Jeroen Ooms<jeroenooms at gmail.com> wrote:
>>> For my GUI, I would like the user to be able to compare a fixed effects
>>> model with a random effects model. For example:
>>> fm0 <- lmer(Reaction ~ Days, sleepstudy);
>>> fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy);
>>> anova(fm0,fm1);
>>>
>>> However, this returns the obvious "No random effects terms specified in
>>> formula" error for the first model. I've also tried fitting the fixed
>>> effects model with lm:
>>>
>>> fm0 <- lm(Reaction ~ Days, sleepstudy);
>>> fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy);
>>> anova(fm0,fm1);
>> Try listing them the other way around
>>
>> anova(fm1, fm0)
>>
>> If the first model in the call to anova is of class "lm" then the
>> method for that class is the one chosen and that method doesn't know
>> about models created by lmer.  You must list them so that the lmer
>> model comes first.
> 
> You could do that ... with lme objects. As explained by the OP, this no
> longer works with lmer objects.
> 
> One may always extract deviances and substract, and make rash hypotheses
> in difference in numerator DFs, but, while this should work (= give
> expected results) with lm and gaussian objects, nothing guarantees that
> glm and lmer use the same parametrizations for non-gaussian models (the
> GLM chapter in V&R4 states that deviances are computed up to an additive
> constant. I tried to follow lmer code, but quickly got lost in C
> code...).
> 
> May this join your (already well-grown) wishlist (along with
> deviance/likelihood under arbitrary parameters, maybe :-) ?
> 
> Sincerely,
> 
> 					Emmanuel Charpentier
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


-- 
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