# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 00-multilevel-imputation.R # # MULTILEVEL MULTIPLE IMPUTATION BY CHAINED EQUATIONS (ML-MICE) # # IPD Meta-Analysis template with system-missing variable support # # # # Data: 2400 patients, 20 studies # # Outcome: y (continuous, 0–30; ~25% missing) # # Covariates: cov1–7 (continuous, 0–10; ~20% missing each) # # New var: comorbidity (binary 0/1; sporadic missings in all studies, # # completely absent in 2 studies) # # Complete: treat, baseline, age, sex # # Cluster: study → id.study (integer) # # # # All methods are single-stage multilevel FCS: # # - y, cov1–7: 2l.pmm (miceadds) — multilevel PMM via blmer # # - comorbidity: 2l.jomo (micemd) — joint multilevel MCMC model # # # # 2l.2stage.* methods are architecturally incompatible with system-missing # # data: they misclassify missingness patterns when any cluster has entirely # # missing data, crashing in sig2 initialisation regardless of workarounds. # # # # References: # # - van Buuren (2018) – Flexible Imputation of Missing Data # # - Resche-Rigon & White (2016) – 2-stage MI for clustered data # # - Quartagno & Carpenter (2016) – jomo multilevel joint imputation # # - BDA3 §18.4: hierarchical imputation when variables are "not asked" # # in some surveys — conceptual basis for cross-study borrowing # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 0. Packages -------------------------------------------------------------- library(readxl) library(dplyr) library(tidyr) library(purrr) library(mice) library(micemd) # provides 2l.jomo library(miceadds) library(mitml) library(foreach) library(doParallel) library(ggplot2) library(blme) # Bayesian regularised mixed models via blmer # 1. Import and stack all sheets ------------------------------------------- data.path = "_data/data.xlsx" sheet.names = excel_sheets(data.path) dat.raw = sheet.names %>% map(~ read_xlsx(data.path, sheet = .x)) %>% bind_rows() # 2. Add comorbidity variable ---------------------------------------------- # comorbidity: binary (0/1), added at patient level. # Missingness structure: # - 2 studies have NO comorbidity data (system-missing / "not assessed") # - All remaining studies have ~15% sporadic (item-level) missings # The two system-missing studies are the first two alphabetically. # Adjust `system.missing.studies` to match your real data situation. set.seed(2024) system.missing.studies = sort(unique(dat.raw$study))[1:2] dat.raw = dat.raw %>% mutate( comorbidity = case_when( # System-missing studies: always NA study %in% system.missing.studies ~ NA_real_, # Other studies: 0/1 with ~15% sporadic item-missing TRUE ~ { n = n() vals = rbinom(n, size = 1, prob = 0.35) # ~35% prevalence vals[runif(n) < 0.15] = NA_real_ # 15% item-missing vals } ) ) # Verify missingness structure dat.raw %>% group_by(study) %>% summarise( n = n(), n.missing = sum(is.na(comorbidity)), pct.missing = round(100 * mean(is.na(comorbidity)), 1), .groups = "drop" ) %>% print(n = Inf) # 3. Encode cluster and coerce types --------------------------------------- # id.study must be a consecutive integer (micemd requirement). # Encode before adding study-level vars so the key is stable. dat.raw = dat.raw %>% mutate(id.study = as.integer(factor(study))) # Create a dictionary for study names and IDs dictionary = table(dat.raw$study, dat.raw$id.study) %>% {rownames(.) -> x; names(x) = colnames(.); x} # 4. Add study-level predictors -------------------------------------------- # Two predictors that are constant within each study (e.g. publication year, # risk-of-bias score). Replace NA_real_ with your actual values. # In the predictor matrix these are set to type 1 (fixed effect only), # not type 2, because they have no within-cluster variance. study.level.vars = tibble( study = sheet.names, year = NA_real_, # ← replace with your actual values, e.g. publication year (mean-centred) rob = NA_real_ # ← replace with your actual values, e.g. risk-of-bias score ) # Example: simulated L2 values (remove and replace with real data) set.seed(42) study.level.vars = study.level.vars %>% mutate( year = as.numeric(scale(seq_along(sheet.names))), # mean-centred publication year rob = sample(0:1, size = n(), replace = TRUE) # risk-of-bias: 0 = low, 1 = high ) dat = dat.raw %>% left_join(study.level.vars, by = "study") %>% select(-study) %>% # keep id.study (integer) as key mutate( treat = factor(treat), # 5-level, complete sex = factor(sex), # binary, complete comorbidity = factor(comorbidity), # binary, system+item missing y = as.numeric(y), # continuous outcome across(cov1:cov7, as.numeric) # continuous covariates # year: numeric, already set in study.level.vars # rob: binary integer (0/1), already set in study.level.vars ) dat %>% summarise(across(everything(), ~ sum(is.na(.)))) %>% print() # 5. System-missing indicator for comorbidity ------------------------------- # Best practice (van Buuren 2018 §2.7; BDA3 §18.4): create a binary # indicator R.comorbidity = 1 if the study assessed comorbidity, 0 if not. # Include as a fixed predictor so that imputed values for system-missing # studies are informed by the known fact that the variable was not collected. # Crucially, R.comorbidity is NEVER itself imputed (method = ""). study.assessed.comorbidity = dat %>% group_by(id.study) %>% summarise(assessed = any(!is.na(comorbidity)), .groups = "drop") dat = dat %>% left_join(study.assessed.comorbidity, by = "id.study") %>% mutate(R.comorbidity = as.integer(assessed)) %>% select(-assessed) %>% select(id.study, everything()) # 6. Predictor matrix ------------------------------------------------------ # Encoding convention (mice / micemd): # -2 = cluster indicator (id.study) # 0 = not used as predictor # 1 = fixed effect predictor # 2 = fixed + random effect predictor (within-cluster random intercept) # comorbidity is excluded as a predictor for the 2l.2stage rows (y, cov1–7). # In the 2 system-missing studies, comorbidity is entirely NA within the # cluster. Including it as a predictor in stage 1 causes lmer to fail and # return NULL for sig2 — this is the source of the sig2 crash, and it # occurs regardless of whether the entry is type 1 or type 2. Removing it # from these rows means it is not used as a within-study predictor; it # remains a predictor in its own row via 2l.2stage.bin. init = mice(dat, maxit = 0, print = FALSE) outlist = init$loggedEvents$out pmatrix = init$predictorMatrix # Cluster variable pmatrix[, "id.study"] = -2 # All active predictors default to type 1 (fixed effect only) pmatrix[pmatrix == 1] = 1 # Only baseline gets a random intercept (type 2) across all imputed rows. # Using type 2 for many predictors risks an unidentifiable between-study # covariance matrix with only 20 clusters; restricting to baseline keeps # the model parsimonious while still accounting for the most important # source of within-study clustering. pmatrix[, "baseline"] = ifelse(pmatrix[, "baseline"] != 0, 2, 0) # Study-level predictors: fixed effect only (no within-cluster variance) pmatrix[, "year"] = ifelse(pmatrix[, "year"] != 0, 1, 0) pmatrix[, "rob"] = ifelse(pmatrix[, "rob"] != 0, 1, 0) # System-missing indicator: fixed effect only, never predicted pmatrix[, "R.comorbidity"] = ifelse(pmatrix[, "R.comorbidity"] != 0, 1, 0) # Exclude comorbidity as predictor for 2l.2stage rows rows.2stage = c("y", paste0("cov", 1:7)) pmatrix[rows.2stage, "comorbidity"] = 0 # Remove variables flagged as problematic (logged events) if (!is.null(outlist) && length(outlist) > 0) { pmatrix[, outlist] = 0 } # Never predict a variable from itself diag(pmatrix) = 0 print(pmatrix[rowSums(pmatrix != 0) > 0, colSums(pmatrix != 0) > 0]) # 7. Imputation methods ---------------------------------------------------- # Method assignment follows the nature of missingness per variable: # 2l.pmm (y, cov1–cov7) # Single-stage multilevel predictive mean matching (Grund et al. 2018). # Fits a mixed model via lmer and draws imputed values from observed # donors close to the predicted mean, keeping imputations within the # observed range by construction. # 2l.2stage.pmm was tested but is architecturally incompatible with # system-missing data. Internally, 2l.2stage.* classifies missingness # per cluster as either sporadic or systematic before fitting stage 1. # When any cluster has entirely missing data for any variable (here: # comorbidity in 2 studies), the method misclassifies the missingness # pattern in other variables within those clusters, crashing in # 2l.2stage.norm.intern regardless of pre-filled initial values. # This cannot be worked around — 2l.2stage.* requires all variables # to have at least some observed data in every cluster. # 2l.jomo (comorbidity) # Single-stage joint multilevel imputation via MCMC (Quartagno & # Carpenter 2016). Handles both sporadic item-missing and the 2 # system-missing studies by fitting a joint model across all clusters # simultaneously, borrowing information from studies that assessed # comorbidity to impute the studies that did not (BDA3 §18.4). # Requires comorbidity to be a factor so jomo routes it into its # categorical (latent normal) component rather than treating it as # continuous. method = init$method method["y"] = "2l.pmm" method[paste0("cov", 1:7)] = "2l.pmm" method["comorbidity"] = "2l.jomo" method[c("treat", "sex", "baseline", "age", "id.study", "year", "rob", "R.comorbidity")] = "" print(method[method != ""]) # 8. Post-processing constraints ------------------------------------------ # For 2l.pmm, squeeze() is redundant — PMM donors are always observed # values. Kept as a safety net for unexpected edge cases. post = make.post(dat) squeeze.range = function(var) { lo = floor(min(dat[[var]], na.rm = TRUE)) hi = ceiling(max(dat[[var]], na.rm = TRUE)) sprintf("imp[[j]][, i] <- squeeze(imp[[j]][, i], c(%d, %d))", lo, hi) } post["y"] = squeeze.range("y") for (v in paste0("cov", 1:7)) post[v] = squeeze.range(v) # comorbidity is factor; squeeze not needed — 2l.jomo routes factors into # jomo's categorical matrix and returns valid 0/1 levels only # 9. Run imputation (parallelised batches) ---------------------------------- # m = 10 batches x 5 imputations = 50 total datasets. # 50 imputations is appropriate for ~25% missingness (van Buuren 2018 §2.8). # blme_use = TRUE replaces lmer with blmer in 2l.pmm, adding weakly # informative priors on variance components to avoid singular fits. # Batch-wise imputation with try() allows recovery from individual failures. n.batches = 10 m.per.batch = 5 maxit = 50 cores = min(detectCores() - 1, n.batches) cl = makeCluster(cores) registerDoParallel(cl) imp.list = foreach( i = seq_len(n.batches), .verbose = TRUE, .errorhandling = "pass" ) %dopar% { require(mice) require(micemd) require(miceadds) require(blme) try( mice( dat, predictorMatrix = pmatrix, method = method, m = m.per.batch, maxit = maxit, post = post, seed = i * 100, blme_use = TRUE, # use blmer instead of lmer in 2l.pmm ) ) } stopCluster(cl) # Report failures failed = sapply(imp.list, inherits, "try-error") if (any(failed)) { warning(sum(failed), " batch(es) failed: ", paste(which(failed), collapse = ", ")) imp.list = imp.list[!failed] } cat("\nSuccessful batches:", length(imp.list), "| Total imputed datasets:", length(imp.list) * m.per.batch, "\n") # 10. Convergence diagnostics ---------------------------------------------- # Trace plots: imputed means and SDs across MICE iterations should show # free mixing with no trends (van Buuren 2018 §6.5). # Density plots: imputed distribution should overlap with observed. # Load imp.list (imputed data element) load("_data/_imputation/imp.list.rda") # Pick one imputed set imp.diag = imp.list[[1]] imputed.vars = c("y", "comorbidity", paste0("cov", 1:7)) plot(imp.diag, y = imputed.vars, layout = c(5, 4)) densityplot(imp.diag, ~ y + cov1 + cov2 + cov3 + cov4, ylab = "Density", scales = list(x = "free")) # Inspect comorbidity imputation in system-missing studies densityplot(imp.diag, ~ comorbidity | treat, ylab = "Density", scales = list(x = "free")) # 11. Combine into mitml list ---------------------------------------------- implist = imp.list %>% map(~ mids2mitml.list(.x)) %>% do.call(c, .) %>% as.mitml.list() # Add study names back & create binary treatment indicator map(implist, function(x) { x$study = recode(x$id.study, !!!dictionary) select(x, study, id.study, everything()) -> x split(x, ~id.study) %>% map_dfr(function(y) { y$treat.bin = as.numeric(as.numeric(factor(y$treat))==1); y }) }) -> implist # Save the implist as a .rda file save(implist, file="_data/_imputation/implist.rda")