On the P-value of the variance of random intercept in glmer model

Question

I'm using a logistic mixed-effect model with random intercept through glmer function from lme4 package. I want to test the significance of the variance of the random intercept to decide if I continue with the same model or switch to a glm. I read about the boundary effect when we test the variance of the random effect. I want to know if it's correct to use for the glmer class the same procedure explained in (Verbenke & Molenberghs, 2000) by dividing the p_value of the Likelihood ratio test from anova by 2.

Answer 1

Based on Molenberghs and Verbeke (2007) (ref below) I believe the "divide by 2" rule still applies (I haven't fully absorbed/carefully read the paper). In any case, that rule applies asymptotically, so I would double-check by doing a parametric bootstrap:

1. fit glm() model (== null model)
2. (for loop x 1000)
   • use simulate() to generate a new set of response variables for the glm()
   • fit both glm and glmer models to this data set
   • compute logLik(glmer)-logLik(glm)
   • store results in delta_nll_vec[i]
3. compute obsval <- logLik(glmer)-logLik(glm) for your observed data
4. the PB p-value is sum(obsval>=delta_nll_vec)/1001

(notes: (a) this is a one-tailed test because log-likelihood of glmer will always be higher than LL of the nested glm unless something goes wrong; (b) dividing by 1001 rather than 1000 because we want to count the observed deviation in the ensemble along with the PB results)

---

Molenberghs, Geert, and Geert Verbeke. “Likelihood Ratio, Score, and Wald Tests in a Constrained Parameter Space.” The American Statistician 61, no. 1 (February 1, 2007): 22–27. https://doi.org/10.1198/000313007X171322.

Answer 2

@BenBolker, thank you very much for your suggestions. When I applied your method: 1/ I got different results between the parametric bootstrap and the method of Molenberghs and Verbeke (devising the p_value by 2): PB_PAV= 0.06193806 where P_val_corrected= 0.4999999. 2/ I got absolutely the same values between the loglik of the GLM and the GLMER inside the loop. Knowing that I observed a difference in loglik between the two models for the observed data. I don’t know if I made something wrong but this is my code :

#First method (P_value correction)
      L_1 <- -2*(logLik(GLM_MODEL)-logLik(GLMM_MODEL)) #The observed LRT 
      p_value_corr <-0.5 * (1 - pchisq(L_1, 1)) #The corrected p_val (p_val/2)
      Data_new <- Data[, -2] # remove the  dependent variable vector 
#Second method parametric boostrap
      Delta <- c(). # LRT vector under H0
      LLGLM <- c() #Loglik vector under H0 for GLM
      LLGLMM <- c()#Loglik vector under H0 for GLMER
      for (p in 1: 1000){
        sim <- simulate(GLM_MODEL, nsim =1)
        names(sim) <- "P_A"
        Data_new <- cbind(Data_new, sim)
        GLM_new <- glm(as.formula(MODEL_GLM), data = Data_new, family=binomial(link="logit"))
        GLMM_new <- glmer(as.formula(MODEL_GLMM),  data = Data_new, family=binomial(link="logit"))
        L_2 <- -2*(logLik(GLM_new)-logLik(GLMM_new))
        LLGLM <- c(LLGLM,logLik(GLM_new))
        LLGLMM <- c(LLGLMM,logLik(GLMM_new))
        Delta <- c(Delta, L_2)
        Data_new <- Data_new[, -ncol(Data_new)] #Remove the dependent variables for the next iteration
      }
      PB_p_value <- sum(Delta>=L_1)/1001
      LL <- cbind(LLGLM,LLGLMM)