library(mice) set.seed(123321) pms <- c(0.05, 0.1, 0.2, 0.4) mps <- c(1, 2) # MCAR and MAR n <- 50 nrep <- 1000 m <- 5 FUN <- function(datainput) { result <- t.test(y ~ x, data = datainput, var.equal = TRUE) means <- result\$estimate diff <- means[1] - means[2] se <- diff / result\$statistic c(diff, se) } ovresult <- NULL for(k in 1:length(pms)) { pm <- pms[k] for(j in 1:length(mps)) { mp <- mps[j] result <- matrix(NA, nrep, 6) for(i in 1:nrep) { y1 <- rnorm(n, 50, 10) y2 <- rnorm(n, 45, 10) y <- c(y1, y2) dx <- rep(c(1, 0), each = n) x <- rep(c(1, 2), each = n) z <- dx + rnorm(n * 2, 0, 1) my <- rep(FALSE, n * 2) while(all(!my)) { if(mp == 1) { my <- as.logical(rbinom(n * 2, 1, pm)) } else { expintcept <- log(pm/(1 - pm)) - (0.5 * 0.5) zhat <- expintcept + 0.5 * z phat <- 1 / (1 + exp(-zhat)) my <- as.logical(rbinom(n, 1, phat)) } } compdat <- data.frame(x, y) result[i, 1:2] <- FUN(compdat) y[my] <- NA dat <- data.frame(x, y, z) # Listwise result[i, 3:4] <- FUN(dat) # mice, MI imp <- mice(dat, m = m, print = FALSE) dat.l <- NULL for(w in 1:m) dat.l[[w]] <- complete(imp, action = w) statmi <- sapply(dat.l, FUN) estmi <- mean(statmi[1,]) bvarmi <- var(statmi[1,]) wvarmi <- sum(statmi[2,]^2)/m tvarmi <- wvarmi + ((m+1)/m)*bvarmi semi <- sqrt(tvarmi) result[i, 5:6] <- c(estmi, semi) } result <- cbind(pm, mp, result) ovresult <- rbind(ovresult, result) } } colnames(ovresult) <- c("pm", "mp", "d", "sed", "dlist", "sedlist", "dmi", "sedmi") ave <- aggregate(cbind(d, sed, dlist, sedlist, dmi, sedmi) ~ pm + mp, data = ovresult, FUN = mean) empse <- aggregate(cbind(d, dlist, dmi) ~ pm + mp, data = ovresult, FUN = sd) conds <- ave[,1:2] relbiasest <- cbind(conds, dlist = (ave[,"dlist"] - 5)/5, dmi = (ave[,"dmi"] - 5)/5) relbiasest2 <- cbind(conds, dlist = (ave[,"dlist"] - ave[,"d"])/ave[,"d"], dmi = (ave[,"dmi"] - ave[,"d"])/ave[,"d"]) round(relbiasest, 4) round(relbiasest2, 4) relbiasse <- cbind(conds, dlist = (ave[,"sedlist"] - empse[,"dlist"])/empse[,"dlist"], dmi = (ave[,"sedmi"] - empse[,"dmi"])/empse[,"dmi"]) relbiasse2 <- cbind(conds, dlist = (ave[,"sedlist"] - empse[,"d"])/empse[,"d"], dmi = (ave[,"sedmi"] - empse[,"d"])/empse[,"d"]) round(relbiasse, 4) round(relbiasse2, 4)