# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#                                                                             #
#   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")
