使用lme4模型和缺失值引导

我正在通过Aguinis, Gottfredson, & Culpepper (2013)的示例进行工作。他们提供了一些R代码在R中执行自举程序来估计斜率变异的置信区间。这是他们的原始R代码：使用lme4模型和缺失值引导

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

来源

2015-11-02 vanao veneri

你能提供一个有效的引用你引用的论文吗？ – Tim

你可以从这里下载[链接]（http://mypage.iu.edu/~haguinis/JOMR.html） –

如果你想通过参数引导获得置信区间，'confint（M03b，method = “boot”）'为你工作？（我认为这些方法可能是新的或更好的发展，因为这篇论文写得很好......） –

（升级到一个答案评论）

如果你想通过参数自举得到置信区间，会为你confint(M03b,method="boot")工作？（我认为这些方法可能是新的或更好的发展，因为这篇论文写得很好......）

来源

2016-08-03 03:00:30

使用lme4模型和缺失值引导

回答

相关问题