[R-sig-ME] Bootstrapping glmer random effects

[Joe King]

Dear all

I am attempting to obtain a bootstrap confidence interval for the random effect in a simple (random intercept) model using glmer.

The problem I have is that the interval I obtain consistently does not contain the value I am trying to get an interval for ! For example I get the following output when I run glmer on the full data:

Generalized linear mixed model fit by the Laplace approximation
Formula: wg~ (1 | city)
   Data: dt
   AIC   BIC logLik deviance
 10115 10131  -5056    10111
Random effects:
 Groups   Name        Variance Std.Dev.
 city(Intercept)       0.14155  0.37623
Number of obs: 19318, groups: city, 218

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -2.58566    0.04045  -63.93   <2e-16 ***

So I am trying to obtain the confidence interval for random effect variance : 0.14155.  Yet, the confidence interval I got was  0.2839343 , 0.3534999. Moreover, the value in every one of the bootstrap replicates is greater than 0.14155. For example, the output from glmer in the last replicate the last bootstrap replicate was

Generalized linear mixed model fit by the Laplace approximation
Formula: wg~ (1 | city)
   Data: sam
   AIC   BIC logLik deviance
 10480 10496  -5238    10476
Random effects:
 Groups   Name        Variance Std.Dev.
 city(Intercept)     0.32769  0.57245
Number of obs: 19318, groups: city, 217

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept) -2.58779    0.05142  -50.33   <2e-16 ***

There are no missing data. This is the code I have used to obtain the interval:

for (i in 1:k) {
    sam <- dt[sample(nrow(dt), replace=T, size=nrow(dt)), ]
    m1<- glmer(wg~(1|city), data=sam, family=binomial)   
    bs[i] <- VarCorr(m1)$city[1]
}
quantile(bs,c(0.025,0.975))

Could anyone suggest why this is happening, and what I might be able to do about it ?

Thank you
JK