Self-tuning with lme4 model and missing values

Question

I am working on an example from Aguinis, Gottfredson and Culpepper (2013) . They provided some R code to perform the boot procedure in R to estimate confidence intervals for slope deviations. This is their source code R:

library(RLRsim)

#STEP 3: Random Intercept and Random Slope model
lmm.fit3=lmer(Y ~ (Xc|l2id) + Xc + I(Wj-mean(Wj)), data=exdata, REML=F)

# Nonparametric Bootstrap Function
REMLVC=VarCorr(lmer(Y ~Xc+(Xc|l2id)+I(Wj-mean(Wj) ),data=exdata,REML=T))$l2id[1:2,1:2]
U.R=chol(REMLVC)

REbootstrap=function(Us,es,X,gs){ 
  nj=nrow(Us)
  idk=sample(1:nj,size=nj,replace=T)
  Usk=as.matrix(Us[idk,])
  esk=sample(es,size=length(es),replace=T)

  S=t(Usk)%*%Usk/nj
  U.S = chol(S)
  A=solve(U.S)%*%U.R
  Usk = Usk%*%A

  datk=expand.grid(l1id = 1:6,l2id = 1:nj)
  colnames(X)=c('one','Xc','Wjc')
  datk=cbind(datk,X)
  datk$yk = X%*%gs + Usk[datk$l2id,1]+Usk[datk$l2id,2]*X[,2]+esk
  lmm.fitk=lmer(yk ~Xc+(Xc|l2id)+Wjc,data=datk,REML=F)
  tau11k = VarCorr(lmm.fitk)$l2id[2,2]
  tau11k
}

# Implementing Bootstrap
bootks=replicate(1500,REbootstrap(Us=ranef(lmm.fit3)$l2id,es=resid(lmm.fit3),X=model.matrix(lmm.fit3),gs=fixef(lmm.fit3)))
quantile(bootks,probs=c(.025,.975))

I tried to adapt the code according to my own data and model. Until now, this has been fruitless, because (a) I do not fully understand all the lines of code and (b) I do not have data points in one of my predictors. Here is what I still have:

#reproducible code
set.seed(855)
exdf <- data.frame(
    ID= c(rep(1:105, 28)),
    content= sort(c(rep(1:28, 105))),
    PrePost= sample(0:1, 105*28, replace=TRUE),
    eyeFRF= sort(rep(rnorm(28), 105)),
    APMs= sample(0:1, 105*28, replace=TRUE),
    Gf= rep(rnorm(105), 28)
)
exdf[which(exdf$ID==62), "eyeFRF"] <- NA
RandomMissing <- sample(rownames(exdf[-which(exdf$ID==62), ]), 17)
exdf[RandomMissing, "eyeFRF"] <- NA 
View(exdf)


#model
M03b <- glmer(APMs ~ PrePost + Gf + eyeFRF + (1|content) + (eyeFRF|ID), data=exdf, family=binomial("logit"))

#own adaptation
REMLVC=VarCorr(M03b)$ID[1:2,1:2]
U.R=chol(REMLVC)

REbootstrap=function(Us, es, X, gs){ 
    #Us = random effects
    #es = residuals
    #X = design matrix
    #gs = fixed effects

    nj = nrow(Us) #104 in this case, one is excluded (#62) b/c no eye-data
    idk = sample(1:nj, size=nj, replace=TRUE) #104 IDs
    Usk = as.matrix(Us[idk,]) #104 intercepts and slopes
    esk = sample(es, size=length(es), replace=TRUE) #2895 datapoints called 'x' (errors?)

    S = t(Usk)%*%Usk/nj #?
    U.S = chol(S) #?
    A = solve(U.S)%*%U.R #?
    Usk = Usk%*%A #?

    datk = expand.grid(content=1:28, ID=1:nj)
    colnames(X) = c('one', 'PrePost', 'Gf', 'eyeFRF')
    datk = cbind(datk, X)
    datk$APMsk = X%*%gs + Usk[datk$ID,1] + Usk[datk$ID,2]*X[ ,2] + esk
    lmm.fitk = glmer(APMsk ~ PrePost + Gf + eyeFRF + (1|content) + (zb|ID), data=datk, family=binomial("logit"))
    tau11k = VarCorr(lmm.fitk)$l2id[2,2]
    tau11k
}

# Implementing Bootstrap
bootks <- replicate(1500, REbootstrap(Us=ranef(M03b)$ID, es=resid(M03b), X=model.matrix(M03b), gs=fixef(M03b)))
quantile(bootks, probs=c(.025,.975))

Answer 1

(update comment to answer)

, confint(M03b,method="boot") ? ( , , ...)