Basic usage
The primary function is misl(). Supply a dataset and
specify:
m– the number of multiply imputed datasets to createmaxit– the number of Gibbs sampling iterations per datasetcon_method,bin_method,cat_method– learners for each outcome type
Use list_learners() to see available options:
| learner | description | continuous | binomial | categorical | package | installed |
|---|---|---|---|---|---|---|
| glm | Linear / logistic regression | TRUE | TRUE | FALSE | stats | TRUE |
| mars | Multivariate adaptive regression splines | TRUE | TRUE | FALSE | earth | TRUE |
| multinom_reg | Multinomial regression | FALSE | FALSE | TRUE | nnet | TRUE |
| rand_forest | Random forest | TRUE | TRUE | TRUE | ranger | TRUE |
| boost_tree | Gradient boosted trees | TRUE | TRUE | TRUE | xgboost | TRUE |
Running misl() with a single learner per outcome type is
the fastest option and is well suited for exploratory work:
misl_imp <- misl(
sim_data,
m = 5,
maxit = 5,
con_method = "glm",
bin_method = "glm",
cat_method = "multinom_reg",
seed = 42,
quiet = TRUE
)misl() returns a list of length m. Each
element contains:
$datasets– the fully imputed tibble$trace– a long-format tibble of convergence statistics
# Number of imputed datasets
length(misl_imp)
#> [1] 5
# Confirm no missing values remain
anyNA(misl_imp[[1]]$datasets)
#> [1] FALSE
# Compare observed vs imputed distribution for sbp
obs_mean <- mean(sim_data$sbp, na.rm = TRUE)
imp_mean <- mean(misl_imp[[1]]$datasets$sbp)
cat("Observed mean sbp (non-missing):", round(obs_mean, 2), "\n")
#> Observed mean sbp (non-missing): 147.31
cat("Imputed mean sbp (all values): ", round(imp_mean, 2), "\n")
#> Imputed mean sbp (all values): 147.65
# Check binary imputations are valid
table(misl_imp[[1]]$datasets$smoker)
#>
#> 0 1
#> 212 88
# Check categorical imputations stay within observed levels
levels(misl_imp[[1]]$datasets$group)
#> [1] "A" "B" "C"Multiple learners and stacking
When multiple learners are supplied, misl uses
cross-validation to learn optimal ensemble weights via the
stacks package. Use cv_folds to reduce the
number of folds and speed up computation:
Analysing the imputed datasets
After imputation, fit your analysis model to each of the
m datasets and pool the results using Rubin’s rules. Here
we implement pooling manually using base R:
# Fit a linear model to each imputed dataset
models <- lapply(misl_imp, function(imp) {
lm(sbp ~ age + bmi + smoker + group, data = imp$datasets)
})
# Pool point estimates and standard errors using Rubin's rules
m <- length(models)
ests <- sapply(models, function(fit) coef(fit))
vars <- sapply(models, function(fit) diag(vcov(fit)))
Q_bar <- rowMeans(ests) # pooled estimate
U_bar <- rowMeans(vars) # within-imputation variance
B <- apply(ests, 1, var) # between-imputation variance
T_total <- U_bar + (1 + 1 / m) * B # total variance
pooled <- data.frame(
term = names(Q_bar),
estimate = round(Q_bar, 4),
std.error = round(sqrt(T_total), 4),
conf.low = round(Q_bar - 1.96 * sqrt(T_total), 4),
conf.high = round(Q_bar + 1.96 * sqrt(T_total), 4)
)
print(pooled)
#> term estimate std.error conf.low conf.high
#> (Intercept) (Intercept) 154.9631 5.1777 144.8148 165.1114
#> age age -0.0421 0.0577 -0.1552 0.0711
#> bmi bmi -0.1890 0.1699 -0.5219 0.1439
#> smoker smoker -0.7820 1.5082 -3.7381 2.1740
#> groupB groupB -0.1393 1.7862 -3.6402 3.3617
#> groupC groupC -0.9780 1.6066 -4.1269 2.1709For a full implementation of Rubin’s rules including degrees of
freedom and p-values, the mice package provides
pool() and can be used directly with misl
output:
Convergence: trace plots
The $trace element records the mean and standard
deviation of imputed values at each iteration, which can be used to
assess convergence:
# Extract trace data across all m datasets
trace <- do.call(rbind, lapply(misl_imp, function(r) r$trace))
# Plot mean of imputed sbp values across iterations for each dataset
sbp_trace <- subset(trace, variable == "sbp" & statistic == "mean")
plot(
sbp_trace$iteration,
sbp_trace$value,
col = sbp_trace$m,
pch = 16,
xlab = "Iteration",
ylab = "Mean imputed sbp",
main = "Trace plot: mean of imputed sbp values",
xaxt = "n"
)
axis(1, at = 1:5)
legend("topright", legend = paste("m =", 1:5),
col = 1:5, pch = 16, cex = 0.8)
Stable traces that mix well across datasets indicate the algorithm has converged.