[R-sig-ME] inference for random effects

[Jeff Evans]

Hi Ben,

I've just tryied to run the simulation approach you posted on glmm.wikidot,
but all the values returned in sim1 are NaN. I think my recoding parallels
your example, but maybe you'd best have a look. You ran your simulation from
the glm object, then refit the full glmer model from the simulation right? 

Also, the file where the function zerodev is said to be (glmmprof.r) on the
wiki page isn't there, though I found it in glmm_Stier.Rnw.

Here's the code I've run so far:

library(lme4)
library(emdbook)
source("C:/Users/Jeff/Documents/Jeffs Work/Analytic
Tools/R_Resources/Bolker/glmmfuns.r")

# Full Model
fm1 = glmer(cbind(SdlFinal, SdlMax-SdlFinal)~ State * Year  +
    (1|Site),
    data = dat.gm,
    family = "binomial")

# Reduced Model
rm1 = glm(cbind(SdlFinal, SdlMax-SdlFinal) ~ State * State,
    data = dat.gm,
    family = "binomial")

zerodev <- function(mm) {
  varprof(mm,0,0,1)[["ML"]]
}


svals <- simulate(rm1, nsim = 1000)
t0 <- system.time(sim1 <- apply(svals, 2, function(x) {
  r <- refit(fm1, x)
  zerodev(r) - deviance(r)
  }))
save("svals", "t0", "sim1", file = "glmersim1.RData")

> str(sim1)
 Named num [1:1000] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
 - attr(*, "names")= chr [1:1000] "V1" "V2" "V3" "V4" ... 


Jeff

#############################################################


-----Original Message-----
From: Ben Bolker [mailto:bolker at ufl.edu] 
Sent: Thursday, February 05, 2009 2:44 PM
To: Juan Pedro Steibel
Cc: Jeff Evans; r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] inference for random effects



Juan Pedro Steibel wrote:
> Jeff,
> Why not use the model without the random effect as the null model?
> JP

  Because the model without the random effect can't be fitted by lmer
(and there is no equivalent of nlme::gls in the lme4 package, which
handles the same syntax but without random effects), and as it turns
out glm calculates the likelihood differently so that they are not
comparable, so that

glm(cbind(successes, total-successes) ~ A * B  data=dat,
           family="binomial")

 would not give an appropriate comparison.

  Jeff's approach is clever, if I (or someone else) gets a
chance it would be nice to compare it against the approach
discussed in http://glmm.wikidot.com/reef-fish , which uses
a modification of some of the code in one of the lme4 vignettes
to compute a likelihood profile for the random effects variance
(including at V=0, which is the likelihood we're interested in
here).

  Ben Bolker

> 
> Jeff Evans wrote:
>> I'm sure this must have been discussed before, but in searching the
archives
>> I haven't found an answer yet. 
>>
>> Simple question:
>>
>> In lme4 can I evaluate the significance of a random effect in a model by
>> substituting an uninformative dummy variable for it and comparing it to
the
>> model with the "real" random effect using anova? 
>>
>> M1 = lmer(cbind(successes, total-successes) ~ A * B + (1|C), data=dat,
>> family="binomial")
>>
>> M2 = lmer(cbind(successes, total-successes) ~ A * B + (1|Cdummy) ,
data=dat,
>> family="binomial")
>>
>> anova(M1,M2)
>>
>> Where A, B, and C are factors, and Cdummy is a column with the word
"dummy"
>> in every row.
>>
>> Then compare the AIC, subtracting 2 from the M2 AIC score since it
"falsely"
>> estimated a parameter for the random effect. When I do this, I get delta
AIC
>> of about 600 favoring the more informative M1. Is this approach
>> fundamentally wrong? 
>>
>>
>> Thanks,
>>
>> Jeff Evans
>> Michigan State University
>>
>> _______________________________________________
>> 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