md.pattern(selfreport[, c("age", "sex", "hm", "hr", "wm", "wr")])
### FIMD Section 7.3.5 Application
bmi <- function(h, w) {
return(w / (h / 100)^2)
}
init <- mice(selfreport, maxit = 0)
meth <- init$meth
meth["bm"] <- "~bmi(hm,wm)"
pred <- init$pred
pred[, c("src", "id", "web", "bm", "br")] <- 0
imp <- mice(selfreport, pred = pred, meth = meth, seed = 66573, maxit = 2, m = 1)
## imp <- mice(selfreport, pred=pred, meth=meth, seed=66573, maxit=20, m=10)
### Like FIMD Figure 7.6
cd <- complete(imp, 1)
xy <- xy.coords(cd$bm, cd$br - cd$bm)
plot(xy,
col = mdc(2), xlab = "Measured BMI", ylab = "Reported - Measured BMI",
xlim = c(17, 45), ylim = c(-5, 5), type = "n", lwd = 0.7
)
polygon(x = c(30, 20, 30), y = c(0, 10, 10), col = "grey95", border = NA)
polygon(x = c(30, 40, 30), y = c(0, -10, -10), col = "grey95", border = NA)
abline(0, 0, lty = 2, lwd = 0.7)
idx <- cd$src == "krul"
xyc <- xy
xyc$x <- xy$x[idx]
xyc$y <- xy$y[idx]
xys <- xy
xys$x <- xy$x[!idx]
xys$y <- xy$y[!idx]
points(xyc, col = mdc(1), cex = 0.7)
points(xys, col = mdc(2), cex = 0.7)
lines(lowess(xyc), col = mdc(4), lwd = 2)
lines(lowess(xys), col = mdc(5), lwd = 2)
text(1:4, x = c(40, 28, 20, 32), y = c(4, 4, -4, -4), cex = 3)
box(lwd = 1)
Run the code above in your browser using DataLab