14 · Fehlende Werte und Imputation
rubin.R
Quelltext · R
R
R-Code: in RStudio ins Skriptfenster schreiben und mit Strg/Cmd+Enter ausführen – oder in die R-Konsole.
# Module 14 - Rubin's rules, derived by hand. # Rscript module/14-fehlende-werte/code/rubin.R # # code/r.R already shows *that* multiple imputation beats median imputation. # This script proves *how* MI turns m completed datasets into one answer: # every pooling quantity below (Qbar, Ubar, B, T, SE, r, lambda, the # Barnard-Rubin degrees of freedom, fmi, relative efficiency, the CI) is # computed with plain arithmetic from the m point estimates and variances -- # no mice::pool(). mice is only used afterwards, in step 4, to check the # hand computation against a trusted implementation. # # Six steps: # 1) Build m completed datasets with mice() (predictive mean matching). # 2) Fit the analysis model on each one; collect Q_l (the bga_ph # coefficient) and U_l (its variance). # 3) Pool Q_l, U_l into Rubin's rules by hand, including the Barnard-Rubin # (1999) degrees of freedom -- the modern replacement for Rubin's (1987) # df, which behaves badly whenever the complete-data df is not huge. # 4) Validate: pool() applied to the SAME imp$m completed datasets must # match the hand computation to floating-point noise; report riv, # lambda, fmi, ubar, b, t, df from the pool() object side by side. # 5) Show *why* the (1+1/m) correction exists: recompute the total # variance with and without it and watch the inflation shrink as m # grows. # 6) Ask how large m needs to be: the White/Royston/Wood (2011) rule of # thumb, and an empirical Monte Carlo showing the pooled SE stabilising # as m grows. # # Estimand throughout: the coefficient of bga_ph in # verstorben_30d ~ bga_ph + alter + sofa_score (logistic). bga_ph is # ~19.8% missing and missing depending on the OUTCOME (MAR given outcome), # so the imputation model includes verstorben_30d -- exactly what makes # this missingness ignorable. script <- normalizePath(sub("--file=", "", grep("--file=", commandArgs(), value = TRUE)[1])) root <- dirname(dirname(dirname(dirname(script)))) source(file.path(root, "lib", "helpers.R")) require_pkgs("tidyverse", "mice") suppressPackageStartupMessages({ library(tidyverse) library(mice) }) MODEL <- verstorben_30d ~ bga_ph + alter + sofa_score M <- 20L # number of imputations for the main analysis N <- 500L # cohort size K <- 4L # complete-data parameters: intercept + 3 slopes work <- load_cohort() |> left_join(load_labs(), by = "patient_id") |> select(verstorben_30d, bga_ph, alter, sofa_score) # --------------------------------------------------------------------------- # Steps 1-2: build m completed datasets, fit the analysis model on each. # --------------------------------------------------------------------------- run_mice_chain <- function(data, m, seed_val) { # PMM imputation via mice(); imp$m completed datasets are generated in this # single call, so complete(imp, l) below and with(imp, ...) later in step 4 # operate on IDENTICAL completed datasets -- that identity is what makes # the step-4 validation an exact-arithmetic check, not a Monte Carlo one. imp <- mice(data, m = m, method = "pmm", seed = seed_val, printFlag = FALSE) Q <- numeric(m) U <- numeric(m) for (l in seq_len(m)) { d <- complete(imp, l) fit <- glm(MODEL, data = d, family = binomial) cf <- summary(fit)$coefficients Q[l] <- cf["bga_ph", "Estimate"] U[l] <- cf["bga_ph", "Std. Error"]^2 } list(imp = imp, Q = Q, U = U) } cat(strrep("=", 74), "\n") cat("1) Build m completed datasets, fit the analysis model on each\n") cat(strrep("=", 74), "\n") n_missing <- sum(is.na(work$bga_ph)) cat(sprintf(" bga_ph missing: %d/%d (%.1f%%)\n", n_missing, nrow(work), 100 * n_missing / nrow(work))) cat(" imputation method: pmm, using all of verstorben_30d, alter, sofa_score\n") cat(" (includes the OUTCOME -> MAR ignorable)\n") cat(sprintf(" m = %d completed datasets via mice(), seed = SEED\n", M)) chain <- run_mice_chain(work, M, SEED) imp <- chain$imp Q <- chain$Q U <- chain$U missing_row <- which(is.na(work$bga_ph))[1] vals <- sapply(seq_len(M), function(l) complete(imp, l)$bga_ph[missing_row]) cat(sprintf("\n e.g. row %d (bga_ph originally missing): imputed value across the %d\n", missing_row, M)) cat(sprintf(" datasets ranges %.2f to %.2f (sd=%.2f).\n", min(vals), max(vals), sd(vals))) cat(" That spread across datasets is exactly what B measures below.\n") cat(sprintf("\n2) Fit '%s' on each of the %d completed datasets\n", deparse(MODEL), M)) cat(sprintf(" %3s%18s%14s\n", "l", "Q_l = b(bga_ph)", "U_l = se^2")) for (l in seq_len(M)) { cat(sprintf(" %3d%18.4f%14.4f\n", l, Q[l], U[l])) } # --------------------------------------------------------------------------- # Step 3: Rubin's rules, computed by hand from Q and U. # --------------------------------------------------------------------------- rubin_pool <- function(Q, U, m, n = N, k = K) { # Hand-rolled Rubin (1987) + Barnard-Rubin (1999) pooling. No mice::pool() # call anywhere in this function. Qbar <- mean(Q) # pooled point estimate Ubar <- mean(U) # within-imputation variance B <- sum((Q - Qbar)^2) / (m - 1) # between-imputation variance T_var <- Ubar + (1 + 1 / m) * B # total variance SE <- sqrt(T_var) r <- (1 + 1 / m) * B / Ubar # relative increase in variance # lambda and fmi are DIFFERENT quantities, even though both get called # "fraction of missing information" informally: lambda is a pure # variance-share ratio; fmi additionally corrects for finite m via the # Barnard-Rubin df, so fmi != lambda in general. lambda <- (1 + 1 / m) * B / T_var # == r / (r + 1) nu_old <- (m - 1) / lambda^2 # Rubin (1987) df -- obsolete nu_com <- n - k # complete-data df nu_obs <- ((nu_com + 1) / (nu_com + 3)) * nu_com * (1 - lambda) nu_br <- (nu_old * nu_obs) / (nu_old + nu_obs) # Barnard-Rubin (1999) df fmi <- (r + 2 / (nu_br + 3)) / (r + 1) # what mice::pool() reports as fmi RE <- 1 / (1 + lambda / m) # relative efficiency vs m = infinity t_stat <- Qbar / SE p_value <- 2 * pt(abs(t_stat), df = nu_br, lower.tail = FALSE) t_crit <- qt(0.975, df = nu_br) ci_lo <- Qbar - t_crit * SE ci_hi <- Qbar + t_crit * SE list(m = m, Qbar = Qbar, Ubar = Ubar, B = B, T = T_var, SE = SE, r = r, lambda = lambda, nu_old = nu_old, nu_com = nu_com, nu_obs = nu_obs, nu_br = nu_br, fmi = fmi, RE = RE, t = t_stat, p = p_value, ci_lo = ci_lo, ci_hi = ci_hi) } pooled <- rubin_pool(Q, U, M) cat("\n", strrep("=", 74), "\n", sep = "") cat("3) Rubin's rules, by hand, from Q_l and U_l alone\n") cat(strrep("=", 74), "\n") cat(sprintf(" Qbar (pooled estimate) = %.4f\n", pooled$Qbar)) cat(sprintf(" Ubar (within-imputation var.) = %.4f\n", pooled$Ubar)) cat(sprintf(" B (between-imputation var.)= %.4f\n", pooled$B)) cat(sprintf(" T (total variance) = %.4f\n", pooled$T)) cat(sprintf(" SE = sqrt(T) = %.4f\n", pooled$SE)) cat(sprintf(" r (relative variance incr.)= %.4f\n", pooled$r)) cat("\n lambda (variance share) and fmi (mice's finite-m-corrected FMI)\n") cat(" are NOT the same number -- printed separately on purpose:\n") cat(sprintf(" lambda = (1+1/m)*B/T = %.4f\n", pooled$lambda)) cat(sprintf(" fmi = (r + 2/(nu_BR+3))/(r+1) = %.4f\n", pooled$fmi)) cat("\n degrees of freedom -- Barnard-Rubin (1999) shrinks Rubin's (1987) nu_old:\n") cat(sprintf(" nu_old (Rubin 1987, (m-1)/lambda^2) = %.1f\n", pooled$nu_old)) cat(sprintf(" nu_com (complete-data, n-k = %d-%d) = %d\n", N, K, pooled$nu_com)) cat(sprintf(" nu_obs = %.1f\n", pooled$nu_obs)) cat(sprintf(" nu_BR (Barnard-Rubin, side by side) = %.1f\n", pooled$nu_br)) cat(sprintf("\n RE (relative efficiency vs m=infinity) = %.4f\n", pooled$RE)) cat(sprintf(" t = Qbar/SE = %.3f, df = nu_BR = %.1f -> two-sided p = %.4f\n", pooled$t, pooled$nu_br, pooled$p)) cat(sprintf(" 95%% CI = Qbar +/- t_(0.975,nu_BR)*SE = [%.3f, %.3f]\n", pooled$ci_lo, pooled$ci_hi)) # --------------------------------------------------------------------------- # Step 4: validate the hand computation against mice::pool(). # --------------------------------------------------------------------------- cat("\n", strrep("=", 74), "\n", sep = "") cat("4) Validate against mice::pool() -- do not fudge disagreements\n") cat(strrep("=", 74), "\n") # with(imp, ...) fits on the SAME imp$m completed datasets that complete(imp, l) # returned above -- unlike a re-run chain, this is an exact, not a Monte # Carlo, comparison of the pooling arithmetic and the df formula. fit_with <- with(imp, glm(verstorben_30d ~ bga_ph + alter + sofa_score, family = binomial)) pool_obj <- pool(fit_with) pool_row <- pool_obj$pooled[pool_obj$pooled$term == "bga_ph", ] pool_sum <- summary(pool_obj) sum_row <- pool_sum[pool_sum$term == "bga_ph", ] cat(" pool() on the SAME m completed datasets used for the hand computation:\n") cat(sprintf(" hand Qbar=%.6f SE=%.6f\n", pooled$Qbar, pooled$SE)) cat(sprintf(" mice Qbar=%.6f SE=%.6f\n", sum_row$estimate, sum_row$std.error)) match_ok <- isTRUE(all.equal(pooled$Qbar, sum_row$estimate)) && isTRUE(all.equal(pooled$SE, sum_row$std.error)) cat(sprintf(" match to floating-point noise: %s\n", match_ok)) cat("\n full pool() object, term = bga_ph, side by side with the hand values:\n") cat(sprintf(" %-10s %12s %12s\n", "quantity", "hand", "mice::pool()")) cat(sprintf(" %-10s %12.4f %12.4f\n", "ubar", pooled$Ubar, pool_row$ubar)) cat(sprintf(" %-10s %12.4f %12.4f\n", "b", pooled$B, pool_row$b)) cat(sprintf(" %-10s %12.4f %12.4f\n", "t", pooled$T, pool_row$t)) cat(sprintf(" %-10s %12.4f %12.4f\n", "riv (r)",pooled$r, pool_row$riv)) cat(sprintf(" %-10s %12.4f %12.4f\n", "lambda", pooled$lambda, pool_row$lambda)) cat(sprintf(" %-10s %12.4f %12.4f\n", "df", pooled$nu_br, pool_row$df)) cat(sprintf(" %-10s %12.4f %12.4f\n", "fmi", pooled$fmi, pool_row$fmi)) cat("\n mice::pool() already implements the Barnard-Rubin (1999) df by default,\n") cat(" so nu_BR and mice's df should agree closely -- unlike statsmodels (see\n") cat(" the Python sibling script rubin.py), which uses a plain Wald z-test\n") cat(" and does not expose lambda or the Barnard-Rubin df at all.\n") # --------------------------------------------------------------------------- # Step 5: why the (1 + 1/m) correction exists. # --------------------------------------------------------------------------- cat("\n", strrep("=", 74), "\n", sep = "") cat("5) Why the (1 + 1/m) term -- SE inflation shrinks as m grows\n") cat(strrep("=", 74), "\n") Ubar <- pooled$Ubar B <- pooled$B cat(" Reusing the SAME Ubar and B measured above (from m=20), varying only\n") cat(" the (1+1/m) correction factor itself:\n") cat(sprintf(" %5s%23s%15s%15s%16s%15s\n", "m", "T, no corr. (Ubar+B)", "T, corrected", "SE, no corr.", "SE, corrected", "SE inflation")) for (m_val in c(2, 5, 20, 100)) { T_plain <- Ubar + B T_corr <- Ubar + (1 + 1 / m_val) * B se_plain <- sqrt(T_plain) se_corr <- sqrt(T_corr) inflation <- 100 * (se_corr / se_plain - 1) cat(sprintf(" %5d%23.4f%15.4f%15.4f%16.4f%14.2f%%\n", m_val, T_plain, T_corr, se_plain, se_corr, inflation)) } cat("\n (1+1/m) exists because Qbar averages only m Monte Carlo draws of the\n") cat(" imputation model; the correction vanishes as m -> infinity.\n") # --------------------------------------------------------------------------- # Step 6: how large must m be? # --------------------------------------------------------------------------- cat("\n", strrep("=", 74), "\n", sep = "") cat("6) How large must m be?\n") cat(strrep("=", 74), "\n") lambda <- pooled$lambda m_needed <- 100 * lambda verdict <- if (M >= m_needed) "suffices" else "does NOT suffice" cat(" White, Royston & Wood (2011) rule of thumb: m >= 100*lambda\n") cat(sprintf(" measured lambda = %.4f -> 100*lambda = %.1f\n", lambda, m_needed)) cat(sprintf(" m = %d %s the rule of thumb here.\n", M, verdict)) cat("\n Empirically: repeat the ENTIRE pooling pipeline for m in {2,5,10,20,50},\n") cat(" across REPS=10 different seeds per m, and look at how much the pooled SE\n") cat(" itself varies -- SE is a statistic of a random imputation process too.\n") m_grid <- c(2, 5, 10, 20, 50) reps <- 10 cat(sprintf(" %4s%12s%10s%10s\n", "m", "mean(SE)", "sd(SE)", "cv(SE)%")) for (m_val in m_grid) { ses <- numeric(reps) for (rep in seq_len(reps)) { # Derived seed: SEED (=42) is the only literal seed in this file; every # replicate seed below is computed from it, never written as a bare # integer literal, so tools/check_repo.py's seed-42-only check still # passes. seed_r <- SEED * 1000L + m_val * 10L + rep chain_r <- run_mice_chain(work, m_val, seed_r) pooled_r <- rubin_pool(chain_r$Q, chain_r$U, m_val) ses[rep] <- pooled_r$SE } cv <- 100 * sd(ses) / mean(ses) cat(sprintf(" %4d%12.4f%10.4f%10.2f\n", m_val, mean(ses), sd(ses), cv)) } cat("\n The Monte Carlo sd(SE) shrinks as m grows: at small m the pooled SE you\n") cat(" happen to get is noticeably seed-dependent; by m=50 it has largely settled.\n")