Missing Data Across Waves: What Will Reviewers Accept?

# Synthetic two-wave study: baseline fully observed, follow-up incomplete.
# Requires mice; no automatic installation. Multilevel/multiwave extensions
# need an imputation model that respects their dependence structure.
if (!requireNamespace('mice',quietly=TRUE)) stop('Install mice before running this example.')
set.seed(8008)
n <- 600
baseline <- rnorm(n)
arm <- rep(0:1,each=n/2)
followup <- .5*arm + .8*baseline + rnorm(n)
missing <- rbinom(n,1,plogis(-.5+.9*baseline+.6*arm))==1
full <- data.frame(arm,baseline,followup)
d <- full; d$followup[missing] <- NA
cc <- lm(followup~arm+baseline,d)
methods <- c(arm='',baseline='',followup='norm')
imp <- mice::mice(d,m=30,maxit=5,method=methods,seed=8009,printFlag=FALSE)
# Rubin pooling for a scalar arm coefficient, with normal-approximation CIs.
# The mice::pool result also supplies finite-df inference for the MAR fit.
pooled <- summary(mice::pool(with(imp,lm(followup~arm+baseline))),conf.int=TRUE)
print(pooled)
scalar <- function(delta) {
 vals <- vapply(seq_len(imp$m),function(i) {
  complete <- mice::complete(imp,i)
  # Controlled shift only for missing outcomes in the intervention arm.
  complete$followup[missing & arm==1] <- complete$followup[missing & arm==1]+delta
  f <- lm(followup~arm+baseline,complete)
  c(q=coef(f)['arm'],u=vcov(f)['arm','arm'])
 },numeric(2))
 q <- mean(vals[1,]); se <- sqrt(mean(vals[2,])+(1+1/imp$m)*var(vals[1,]))
 c(delta=delta,estimate=q,lower=q-qnorm(.975)*se,upper=q+qnorm(.975)*se)
}
cat(sprintf('Enrolled=%d; follow-up observed=%d; missing control=%d; missing intervention=%d\n',n,sum(!missing),sum(missing & arm==0),sum(missing & arm==1)))
cat(sprintf('Complete-data arm=%.3f; complete-case adjusted arm=%.3f\n',coef(lm(followup~arm+baseline,full))['arm'],coef(cc)['arm']))
print(t(vapply(c(0,-.5,-1),scalar,numeric(4))))
cat('mice version:',as.character(packageVersion('mice')),'\n')
stopifnot(imp$m==30,all(is.finite(pooled$estimate)))

# Maintainer-only export; run from repository root. Synthetic data only.
if (Sys.getenv("CLINIC_EXPORT_DATA") == "1") {
write.csv(d,'analysis-clinic/examples/case-008-shared-data.csv',row.names=FALSE,quote=FALSE,na='')
}

# DASS Analysis Clinic Case 008; Python dependencies: numpy, scipy.
from pathlib import Path
import numpy as np
from scipy import stats, optimize
HERE = Path(__file__).resolve().parent

def ols(X, y):
    b = np.linalg.lstsq(X, y, rcond=None)[0]
    residual = y - X @ b
    df = len(y) - X.shape[1]
    cov = (residual @ residual / df) * np.linalg.inv(X.T @ X)
    se = np.sqrt(np.diag(cov))
    ci = np.column_stack((b-stats.t.ppf(.975, df)*se,b+stats.t.ppf(.975,df)*se))
    return b, cov, ci

def power(n, d):
    # Matches R power.t.test(strict=FALSE): rejection tail in effect direction.
    return stats.nct.sf(stats.t.ppf(.975,2*n-2),2*n-2,d*np.sqrt(n/2))

# Bayesian normal-regression MI for one incomplete continuous outcome.
# Fully observed arm/baseline, flat coefficient prior and p(sigma^2) proportional
# to 1/sigma^2. Not a general MICE or multilevel imputation implementation.
d=np.genfromtxt(HERE/'case-008-shared-data.csv',delimiter=',',names=True)
arm=d['arm']; baseline=d['baseline']; y=d['followup']; missing=np.isnan(y)
X=np.column_stack((np.ones(len(y)),arm,baseline)); observed=~missing
Xo=X[observed]; yo=y[observed]; bhat=np.linalg.lstsq(Xo,yo,rcond=None)[0]
df=len(yo)-X.shape[1]; sse=np.sum((yo-Xo@bhat)**2); inv=np.linalg.inv(Xo.T@Xo)
rng=np.random.default_rng(8009); completed=[]
for _ in range(30):
 variance=sse/rng.chisquare(df)
 beta=rng.multivariate_normal(bhat,variance*inv)
 out=y.copy(); out[missing]=X[missing]@beta+rng.normal(0,np.sqrt(variance),missing.sum())
 completed.append(out)
def pooled(delta):
 estimates=[]; variances=[]
 for out in completed:
  shifted=out.copy(); shifted[missing & (arm==1)]+=delta
  b,cov,_=ols(X,shifted); estimates.append(b[1]);variances.append(cov[1,1])
 q=np.mean(estimates); se=np.sqrt(np.mean(variances)+(1+1/30)*np.var(estimates,ddof=1))
 return q,q-stats.norm.ppf(.975)*se,q+stats.norm.ppf(.975)*se
print('Observed:',observed.sum(),'missing by arm:',[np.sum(missing & (arm==g)) for g in [0,1]])
print('Complete-case adjusted effect:',bhat[1])
for delta in [0,-.5,-1]:print('delta / pooled effect / approximate CI:',delta,pooled(delta))
print('Python posterior draws differ from mice; the page table reports the R run.')
assert observed.sum()==321 and np.isclose(bhat[1],.474,atol=.0005)
assert pooled(0)[0]>pooled(-.5)[0]>pooled(-1)[0]