使用lme4模型和缺失值进行自举

2023年5月16日 319次阅读

我正在研究
Aguinis, Gottfredson, & Culpepper (2013)的一个例子.他们提供了一些R代码来在R中执行自举程序来估计斜率方差的置信区间.这是他们原来的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))

我试图调整代码以适应我自己的数据和模型.到目前为止,这是无用的,因为(a)我不完全理解所有代码行和(b)我在一个预测器中缺少数据点.这是我到目前为止：

#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))

最佳答案 (将评论升级为答案)

如果您试图通过参数自举获得置信区间,那么confint(M03b,method =“boot”)是否适合您？ (我认为这些方法可能是新的或更好的开发,因为那篇论文是写的…)

点赞