使用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模型和缺失值引導

回答

相關問題