mi-predictions.md

Can we get away with single imputation when the goal is to obtain well-calibrated prediction intervals

Of course not! 😁

library(foreach)

set.seed(123)

nsim <- 200
cov1 <- covm <- numeric(nsim)

predfunc <- \(x, fit, alpha = 0.05) {
  tval <- qt(1 - alpha/2, fit$df.residual)
  summfit <- summary(fit)
  sigma <- summfit$sigma
  xmat <- model.matrix(~x)
  se <- sqrt(1 + rowSums((xmat %*% summfit$cov.unscaled) * xmat))
  xfit <- xmat %*% fit$coefficients
  return(list(fit = cbind(xfit, xfit, xfit) + (tval * se * sigma) %x% t(c(0, -1, 1)),
              var = se^2 * sigma^2))
}

N <- 30
P <- 5
V <- 0.5 + diag(0.5, P)

cl <- parallel::makeCluster(parallel::detectCores() - 4)
doParallel::registerDoParallel(cl)

x <- foreach (i = 1:nsim, .combine = 'c') %dopar% {
cat("\rIteration: ", i)

Xtrain <- matrix(rnorm(N * P), N, P) %*% chol(V)
b <- sqrt(0.5 / c(1,2,3,4,5) %*% V %*% c(1,2,3,4,5))
ytrain <- Xtrain %*% (c(1,2,3,4,5)*c(b)) + rnorm(N, 0, sqrt(1-0.5))

train <- data.frame(y = ytrain, Xtrain[,-5])
fit <- lm(y ~ ., data = train)

#preds <- predict(fit, interval = "pred", se = TRUE)

#preds$fit |> head()

#predfunc(Xtrain[,-5], fit)$fit |> head()
#preds$fit |> head()


Ntest <- 10000
Xtest <- matrix(rnorm(Ntest * P), Ntest, P) %*% chol(V)
ytest <- Xtest %*% (c(1,2,3,4,5)*c(b)) + rnorm(Ntest, 0, sqrt(1-0.5))


#predint <- predfunc(Xtest[,-5], fit)
#mean(ytest > predint$fit[,2] & ytest < predint$fit[,3])


m <- 10
Xamp <- mice::ampute(train, prop = 0.5)
Ximp1 <- mice::mice(Xamp$amp, m = 1, method = "norm.predict", print = FALSE)
Ximpm <- mice::mice(Xamp$amp, m = m, method = "norm", print = FALSE)

fitimp1 <- with(Ximp1, lm(y ~ X1 + X2 + X3 + X4))$analyses
fitimpm <- with(Ximpm, lm(y ~ X1 + X2 + X3 + X4))$analyses

impvals1 <- purrr::map(fitimp1, \(x) {
  predint <- predfunc(Xtest[,-5], x)
  data.frame(
    q = predint$fit[,1],
    u = predint$var
  )
}) |>
  unlist() |>
  array(dim = c(Ntest, 2, m)) |>
  apply(1, \(x) {
    pooled <- mice::pool.scalar(x[1,], x[2,])
    if (is.na(pooled$b)) {
      pooled$t <- pooled$ubar
      pooled$df <- fit$df.residual
    }
    pooled$qbar + c(-1, 1) * sqrt(pooled$t) * qt(0.975, pooled$df)
  }) |>
  t()

impvals1m <- purrr::map(fitimpm, \(x) {
  predint <- predfunc(Xtest[,-5], x)
  data.frame(
    q = predint$fit[,1],
    u = predint$var
  )
}) |>
  unlist() |>
  array(dim = c(Ntest, 2, m)) |>
  apply(1, \(x) {
    pooled <- mice::pool.scalar(x[1,], x[2,])
    if (is.na(pooled$b)) {
      pooled$t <- pooled$ubar
      pooled$df <- fit$df.residual
    }
    pooled$qbar + c(-1, 1) * sqrt(pooled$t) * qt(0.975, pooled$df)
  }) |>
  t()

c(imp1 = mean(ytest > impvals1[,1] & ytest < impvals1[,2]),
  impm = mean(ytest > impvals1m[,1] & ytest < impvals1m[,2]))
}

# Single imputation
mean(x[2*0:(nsim-1)+1])
#> [1] 0.8638355
# Multiple imputation
mean(x[2*1:nsim])
#> [1] 0.9328245

^{Created on 2024-10-15 with reprex v2.0.2}