On the P-value of the variance of random intercept in glmer model 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.
 A: 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:

*

*fit glm() model (== null model)

*(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]



*compute obsval <- logLik(glmer)-logLik(glm) for your observed data

*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.
A: @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)

