# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 08-ipd-multivariate-network-ma.R # # Two-stage Bayesian IPD network meta-analysis (IPD-NMA), continuous # # Design follows the multivariate NMA framework of: # # Simmonds & Higgins (2007), Riley et al. (2010), Debray et al. (2015) # # # # STAGE 1 — Per-study individual-level models (runjags, method = "rjags") # # For each study i and each imputation m, fit a JAGS linear model: # # y[j] ~ Normal(mu[j], sigma^2) # # mu[j] = intercept + trt.eff[trt[j]] # # + inprod(beta[1:Nprog], prog.mat[j,]) # prognostic # # + inprod(tmod[trt[j], 1:Nprog], prog.mat[j,]) # moderators # # trt.eff[1] = 0 (within-study reference arm = locally lowest treatment) # # tmod[1, p] = 0 for all p (interactions reference-anchored) # # Pool Stage 1 MCMC draws across all M imputations by concatenation. # # Compute per-study mean vector theta.hat[i] and sampling covariance Sigma[i] # # # # STAGE 2 — Multivariate NMA (runjags, method = "parallel") # # For each study i with na[i] arms, Stage 1 parameter vector has dimension # # D = (na-1) * (1 + Nprog): positions 1..(na-1) = trt.eff[2..na]; # # positions (na-1)*p+1 .. (na-1)*(p+1) (p=1..Nprog) = tmod[2..na, p]. # # Likelihood: y.hat[i] ~ dmnorm(Mu[i], inv(Sigma[i])) where Mu[i] contains: # # - treatment effect dimensions: random effects with heterogeneity tau, # # multi-arm precision (Lu & Ades 2006) # # - interaction dimensions: fixed pooled parameters tmod.global[t, p] # # Priors (Gelman 2006; BUGS Book §5.6; BDA3 §5.7): # # d[t] ~ Normal(0, 1), tmod.global[t,p] ~ Normal(0, 0.1), tau ~ Unif(0,2) # # # # Covariates detected from data automatically; Nprog never hardcoded. # # Within-study centring (Simmonds & Higgins 2007) separates individual-level # # from between-study ecological effects. # # Multiple imputation: implist = list of M fully-imputed long-format # # data frames (study, treat, y, baseline, age, sex, cov1..covK). # # Stage 1 pools imputations by concatenating MCMC draws before mean/ # # covariance — propagates MI uncertainty into Stage 2 (Rubin 1987). # # Assumptions: treatment codes 1..Nt; treatment 1 = global reference; # # all covariates fully imputed; no structural missingness in y. # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 0. Packages -------------------------------------------------------------- library(runjags) # run.jags() for both stages library(coda) library(dplyr) library(ggplot2) library(purrr) library(matrixStats) library(HDInterval) path = "_data/_imputation/implist.rda" load(path) # Select first 25 sets to speed up computations implist = implist[1:25] # 1. Global settings -------------------------------------------------------- set.seed(2026) ref.trt = 1L # global reference treatment # Stage 1 MCMC (per-study; method = "rjags" avoids process-spawn overhead) s1.adapt = 1000 s1.burn = 3000 s1.iter = 5000 s1.thin = 5 s1.chains = 3 Nprog = 10 # number of predictors to be used arm.counts = c(2,3) # number of arm counts in data # Stage 2 MCMC s2.adapt = 3000 s2.burn = 10000 s2.iter = 20000 s2.thin = 10 s2.chains = 4 # 2. Covariate detection and within-study centring helpers ------------------- # Detect covariate columns from a data frame and return their names in a # canonical order: baseline, age, sex, then all cov* columns (sorted). detect.prog.cols = function(df) { fixed = c("baseline", "age", "sex") extras = sort(grep("^cov[0-9]+$", names(df), value = TRUE)) c(fixed, extras) } # Centre covariates within each study (subtract study mean for continuous; # grand-mean centre for sex since it is binary) and scale by pooled # within-study SD. Returns the augmented data frame plus scaling info. centre.and.scale = function(df, sds = NULL) { prog.cols = detect.prog.cols(df) study.levs = sort(unique(df$study)) df = df %>% group_by(study) %>% # All factors have to be converted to numeric first mutate(sex = sex %>% as.character %>% as.numeric) %>% # Center and scale all other variables mutate( baseline.c = baseline - mean(baseline, na.rm = TRUE), age.c = age - mean(age, na.rm = TRUE), across(all_of(setdiff(prog.cols, c("baseline","age","sex"))), ~ .x - mean(.x, na.rm = TRUE), .names = "{.col}.c") ) %>% ungroup() %>% # Binary variable(s), in this case sex, are grand mean centered mutate(sex.c = sex - mean(sex, na.rm = TRUE)) # Compute or re-use pooled within-study SDs cent.cols = c("baseline.c", "age.c", paste0(setdiff(prog.cols, c("baseline","age","sex")), ".c")) if (is.null(sds)) { sds = sapply(cent.cols, function(v) sd(df[[v]], na.rm = TRUE)) } for (col in cent.cols) { sc = sub("\\.c$", ".cs", col) s = sds[col] df[[sc]] = if (s > 0) df[[col]] / s else df[[col]] } list(df = df, sds = sds, prog.cols = prog.cols) } # Build the [N x Nprog] prognostic matrix from a centred-and-scaled data frame. make.prog.mat = function(df.cs, prog.cols) { extra.cols = setdiff(prog.cols, c("baseline", "age", "sex")) mat.cols = c("baseline.cs", "age.cs", "sex.c", paste0(extra.cols, ".cs")) mat = as.matrix(df.cs[, mat.cols, drop = FALSE]) if (any(is.na(mat))) { warning("NAs in prognostic matrix — replacing with 0.") mat[is.na(mat)] = 0 } mat } # 3. Dry run on first imputation to fix Nprog, Nt, scaling SDs -------------- # Make sure that treat is an integer map(implist, ~mutate(.,treat = treat %>% as.character %>% as.integer)) -> implist dry.df = implist[[1]] %>% filter(!is.na(y)) dry.cs = centre.and.scale(dry.df) prog.cols = dry.cs$prog.cols Nprog = length(prog.cols) # total moderators; derived from data sds.ref = dry.cs$sds # pooled within-study SDs (used in all imputations) Nt = max(dry.df$treat, na.rm = TRUE) M = length(implist) # Grand means for each covariate (used as centering reference in prediction) grand.means = sapply(prog.cols, function(v) mean(as.numeric(as.character(dry.df[[v]])), na.rm = TRUE)) cat(sprintf("Data: %d studies | %d treatments | %d prognostic factors\n", length(unique(dry.df$study)), Nt, Nprog)) cat("Prognostic factors:", paste(prog.cols, collapse = ", "), "\n") # 4. Stage 1 JAGS model (single generic model for any arm count) ------------- # na and Nprog are passed as JAGS data; trt.eff[1]=0, tmod[1,p]=0 are constraints. model.stage1 = " model { for (j in 1:N) { y[j] ~ dnorm(mu[j], prec.y) mu[j] <- intercept + trt.eff[trt[j]] + inprod(beta[1:Nprog], prog.mat[j, 1:Nprog]) + inprod(tmod[trt[j], 1:Nprog], prog.mat[j, 1:Nprog]) } prec.y <- 1 / (sigma * sigma) sigma ~ dunif(0, 5) intercept ~ dnorm(0, 0.01) # Treatment effects; arm 1 (within-study reference) fixed to 0 trt.eff[1] <- 0 for (k in 2:na) { trt.eff[k] ~ dnorm(0, 0.01) } # Prognostic main effects and moderator interactions; arm 1 fixed to 0 for (p in 1:Nprog) { beta[p] ~ dnorm(0, 0.01) tmod[1, p] <- 0 for (k in 2:na) { tmod[k, p] ~ dnorm(0, 0.01) } } } " # 5. Stage 1 fitting: per study, pool across imputations --------------------- # For each study, concatenate MCMC draws from all M imputations → theta.hat, inv.Sigma. # Parameter vector: positions 1..(na-1) = trt.eff[2..na]; then tmod[2..na, p] by p. # Helper: build column names for the Stage 1 parameter vector in the correct order stage1.param.names = function(na, Nprog) { nm1 = na - 1 trt.cols = paste0("trt.eff[", 2:na, "]") mod.cols = as.vector(outer(2:na, seq_len(Nprog), function(k, p) paste0("tmod[", k, ",", p, "]"))) # outer(rows=k, cols=p) flattened column-first → for p=1: k=2..na; for p=2: k=2..na; ... c(trt.cols, mod.cols) } # All unique studies and their arm structures (same across imputations) study.info = implist[[1]] %>% group_by(study) %>% summarise( na = n_distinct(treat), t.global = list(sort(unique(treat))), # global treatment codes, sorted .groups = "drop" ) %>% arrange(study) ns = nrow(study.info) study.ids = study.info$study cat(sprintf("\nFitting Stage 1 to %d studies x %d imputations.\n", ns, M)) # Storage for Stage 2 inputs (grouped by arm count later) s1.results = vector("list", ns) for (si in seq_len(ns)) { study.name = study.ids[si] na = study.info$na[si] D = (na - 1) * (1 + Nprog) t.global = study.info$t.global[[si]] # length na, globally coded cat(sprintf(" Study %-20s na = %d D = %d\n", study.name, na, D)) # Extract and fit each imputation, collect MCMC samples samps.imp = vector("list", M) for (m in seq_len(M)) { df.study = implist[[m]] %>% filter(study == study.name, !is.na(y)) # Within-study centre and scale using shared reference SDs cs = centre.and.scale(df.study, sds = sds.ref) prog.mat = make.prog.mat(cs$df, prog.cols) # Local treatment coding: sorted global codes → positions 1..na trt.local = match(df.study$treat, t.global) jags.data = list( N = nrow(df.study), na = na, Nprog = Nprog, y = df.study$y, trt = trt.local, prog.mat = prog.mat ) fit = run.jags(model = model.stage1, data = jags.data, monitor = c("trt.eff", "tmod"), n.chains = s1.chains, adapt = s1.adapt, burnin = s1.burn, sample = s1.iter, thin = s1.thin, summarise = FALSE, plots = FALSE) samps.imp[[m]] = as.matrix(as.mcmc.list(fit)) } # Pool imputations by concatenating draws samps.all = do.call(rbind, samps.imp) # Extract parameters in Stage 2 vector order col.order = stage1.param.names(na, Nprog) samps.vec = samps.all[, col.order, drop = FALSE] theta.hat = matrixStats::colMedians(samps.vec) Sigma.hat = cov(samps.vec) # Regularise covariance matrix before inversion (add small ridge) Sigma.reg = Sigma.hat + diag(1e-6, D) inv.Sigma = solve(Sigma.reg) s1.results[[si]] = list( study = study.name, na = na, D = D, t.global = t.global, # global treatment codes for arms 1..na theta.hat = theta.hat, Sigma.hat = Sigma.hat, inv.Sigma = inv.Sigma ) } # 6. Assemble Stage 2 inputs (group studies by arm count) -------------------- # For each arm count present in the data, build: # y_na [N_na x D_na] — Stage 1 mean vectors # invV_na [N_na x D_na x D_na] — Stage 1 inverse covariance matrices # t_na [N_na x na] — global treatment codes per arm position arm.counts = sort(unique(sapply(s1.results, `[[`, "na"))) # e.g. c(2, 3) s2.inputs = lapply(arm.counts, function(na) { idx = which(sapply(s1.results, `[[`, "na") == na) N = length(idx) D = (na - 1) * (1 + Nprog) y.mat = matrix(NA_real_, nrow = N, ncol = D) invV.arr = array(NA_real_, dim = c(N, D, D)) t.mat = matrix(NA_integer_, nrow = N, ncol = na) for (ii in seq_len(N)) { r = s1.results[[idx[ii]]] y.mat[ii, ] = r$theta.hat invV.arr[ii,,] = r$inv.Sigma t.mat[ii, ] = r$t.global } list(na = na, N = N, D = D, y = y.mat, invV = invV.arr, t = t.mat) }) names(s2.inputs) = as.character(arm.counts) cat(sprintf("\nStage 2 groups: %s\n", paste(sapply(s2.inputs, function(g) sprintf("%d-arm: %d studies", g$na, g$N)), collapse = " | "))) # 7. Stage 2 JAGS model builder ---------------------------------------------- # Generates a JAGS model string for any set of arm-count groups. # Parameter vector convention (for a na-arm study, nm1 = na-1): # positions 1..nm1 : treatment effects # positions nm1*p + kk # (p=1..Nprog, kk=1..nm1) : interaction for covariate p, arm (kk+1) # Multi-arm treatment-effect precision matrix (Lu & Ades 2006): # Cor(delta_j, delta_l) = 0.5 for all j != l # => Sigma_RE = tau^2 * (0.5*J + 0.5*I) # => Prec_RE = (1/tau^2) * inv(0.5*J + 0.5*I) # Diagonal : 2*(na-1)/na * prec (= 4/3 * prec for na=3) # Off-diagonal: -2/na * prec (= -2/3 * prec for na=3) build.stage2.model = function(arm.counts, Nprog) { L = function(...) paste0(..., "\n") # line builder out = L("model {") out = paste0(out, L("")) # Global priors out = paste0(out, L("## Pooled treatment effects (d[1] = 0 by definition)"), L(" d[1] <- 0"), L(" for (t in 2:Nt) {"), L(" d[t] ~ dnorm(0, 1)"), L(" }"), L(""), L("## Pooled moderator effects (tmod.global[1,p] = 0 by definition)"), L(" for (p in 1:Nprog) {"), L(" tmod.global[1, p] <- 0"), L(" }"), L(" for (t in 2:Nt) {"), L(" for (p in 1:Nprog) {"), L(" tmod.global[t, p] ~ dnorm(0, 0.1)"), L(" }"), L(" }"), L(""), L("## Between-study heterogeneity in treatment effects"), L(" tau ~ dunif(0, 2)"), L(" prec <- 1 / (tau * tau)"), L("") ) for (na.str in as.character(sort(arm.counts))) { na = as.integer(na.str) nm1 = na - 1 D = nm1 * (1 + Nprog) # Variable name suffixes for this arm count suf = na.str # "2", "3", etc. Nv = paste0("N", suf) yv = paste0("y", suf) iVv = paste0("invV", suf) tv = paste0("t", suf) Mv = paste0("Mu", suf) Dv = paste0("D", suf) out = paste0(out, L(sprintf("## %d-arm studies ", na))) out = paste0(out, L(sprintf(" for (i in 1:%s) {", Nv))) out = paste0(out, L(sprintf(" %s[i, 1:%s] ~ dmnorm(%s[i, 1:%s], %s[i, 1:%s, 1:%s])", yv, Dv, Mv, Dv, iVv, Dv, Dv))) if (na == 2) { # Scalar random effect for the single treatment contrast out = paste0(out, L(sprintf(" %s[i, 1] ~ dnorm(d[%s[i,2]] - d[%s[i,1]], prec)", Mv, tv, tv))) } else { # Multivariate random effect for the nm1 treatment contrasts thv = paste0("theta", suf) prcv = paste0("prec.", suf) # Precision matrix entries (defined once outside study loop, see below) out = paste0(out, L(sprintf(" %s[i, 1:%d] ~ dmnorm(%s[i, 1:%d], %s[1:%d, 1:%d])", Mv, nm1, thv, nm1, prcv, nm1, nm1)), L(sprintf(" for (kk in 1:%d) {", nm1)), L(sprintf(" %s[i, kk] <- d[%s[i, kk+1]] - d[%s[i, 1]]", thv, tv, tv)), L(" }") ) } # Interaction terms: position nm1*p + kk (p=1..Nprog, kk=1..nm1) out = paste0(out, L(" for (p in 1:Nprog) {"), L(sprintf(" for (kk in 1:%d) {", nm1)), L(sprintf(" %s[i, %d*p + kk] <- tmod.global[%s[i, kk+1], p] - tmod.global[%s[i, 1], p]", Mv, nm1, tv, tv)), L(" }"), L(" }") ) out = paste0(out, L(" }"), L("")) } # Precision matrices for multi-arm groups (defined once in global scope) for (na.str in as.character(sort(arm.counts[arm.counts >= 3]))) { na = as.integer(na.str) nm1 = na - 1 prcv = paste0("prec.", na.str) diag.v = 2 * nm1 / na # e.g. 4/3 for na=3 offdiag.v = -2 / na # e.g. -2/3 for na=3 out = paste0(out, L(sprintf("## Precision matrix for %d-arm RE (Cor=0.5, diag=%.4g*prec, off=%.4g*prec)", na, diag.v, offdiag.v))) for (j in seq_len(nm1)) { for (l in seq_len(nm1)) { val = if (j == l) diag.v else offdiag.v out = paste0(out, L(sprintf(" %s[%d,%d] <- %.10g * prec", prcv, j, l, val))) } } out = paste0(out, L("")) } paste0(out, "}\n") } model.stage2 = build.stage2.model(arm.counts, Nprog) cat("\n--- Stage 2 JAGS model ---\n") cat(model.stage2) # 8. Assemble Stage 2 JAGS data list ---------------------------------------- # Collect all relevant data s2.data = list(s2.inputs = s2.inputs, Nt = Nt, Nprog = Nprog) # To "skip" L1 for the example, load this data list path = "_data/_mvnma/s2.data.rda" load(path) jags.data.s2 = list(Nt = s2.data$Nt, Nprog = s2.data$Nprog) for (na.str in names(s2.data$s2.inputs)) { g = s2.data$s2.inputs[[na.str]] na = g$na D = g$D suf = na.str jags.data.s2[[paste0("N", suf)]] = g$N jags.data.s2[[paste0("y", suf)]] = g$y jags.data.s2[[paste0("invV", suf)]] = g$invV jags.data.s2[[paste0("t", suf)]] = g$t jags.data.s2[[paste0("D", suf)]] = D } # 9. Stage 2 initial values -------------------------------------------------- make.s2.inits = function(Nt, Nprog, n.chains, seed = 1L) { set.seed(seed) lapply(seq_len(n.chains), function(ch) { d_init = rnorm(Nt, 0, 0.2) d_init[ref.trt] = NA # skip fixed reference node tmod_init = matrix(0, nrow = Nt, ncol = Nprog) tmod_init[ref.trt, ] <- NA # skip fixed reference row list( d = d_init, tmod.global = tmod_init, tau = runif(1, 0.1, 0.5) )}) } # 10. Stage 2 fitting (runjags, parallel chains) ---------------------------- s2.params = c("d", "tmod.global", "tau") fit.s2 = run.jags(model = model.stage2, data = jags.data.s2, monitor = s2.params, inits = make.s2.inits(s2.data$Nt, s2.data$Nprog, s2.chains), n.chains = s2.chains, adapt = s2.adapt, burnin = s2.burn, sample = s2.iter, thin = s2.thin, summarise = FALSE) samp.s2 = as.mcmc.list(fit.s2) post.mat = as.matrix(samp.s2) # 11. Convergence diagnostics ------------------------------------------------ cat("\n--- Gelman-Rubin R-hat (target < 1.1) ---\n") rhat = gelman.diag(samp.s2, multivariate = FALSE)$psrf[, 1] %>% {.[!is.na(.)]} top.rhat = sort(rhat, decreasing = TRUE)[seq_len(min(10, length(rhat)))] print(round(top.rhat, 3)) cat(sprintf("%d / %d parameters with R-hat < 1.1\n", sum(rhat < 1.1, na.rm = TRUE), length(rhat))) ess = effectiveSize(samp.s2) %>% {.[!.==0]} cat(sprintf("\nESS — Min: %d Median: %d\n", round(min(ess)), round(median(ess)))) # 12. Posterior summary helper ----------------------------------------------- posterior.summary = function(mat, params = NULL, prob = 0.95) { if (!is.null(params)) { params = intersect(params, colnames(mat)) mat = mat[, params, drop = FALSE] } as.data.frame(t(apply(mat, 2, function(x) { h = HDInterval::hdi(x, credMass = prob) c( mean = mean(x), sd = sd(x), hdi.lo = unname(h["lower"]), median = unname(quantile(x, 0.500)), hdi.hi = unname(h["upper"]), P.pos = mean(x > 0) ) }))) } trt.labels = paste0("Trt", seq_len(s2.data$Nt)) prog.labels = c("baseline", "age", "sex", "cov1", "cov2", "cov3", "cov4", "cov5", "cov6", "cov7") # 13. Printed posterior summaries -------------------------------------------- { cat("\n--- Treatment effects vs reference (d[t]) ---\n") d.params = paste0("d[", seq_len(s2.data$Nt), "]") ps.d = posterior.summary(post.mat, d.params) rownames(ps.d) = trt.labels print(round(ps.d, 3)) cat("\n--- Pooled moderator effects (tmod.global[t, p]) ---\n") for (t in seq_len(s2.data$Nt)) { if (t == ref.trt) next gp = paste0("tmod.global[", t, ",", seq_len(Nprog), "]") gp = intersect(gp, colnames(post.mat)) gs = posterior.summary(post.mat, gp) rownames(gs) = prog.labels[seq_along(gp)] cat(sprintf("\n Trt%d vs reference:\n", t)) print(round(gs, 3)) } cat("\n--- Heterogeneity ---\n") print(round(posterior.summary(post.mat, "tau"), 3)) } # 14. Forest plot: treatment effects ------------------------------------------ ps.trt = posterior.summary(post.mat, d.params) ps.trt$treatment = factor(trt.labels, levels = rev(trt.labels)) ns = Reduce("+", lapply( Filter(function(k) grepl("^[0-9]+$", k), names(s2.data$s2.inputs)), function(k) s2.data$s2.inputs[[k]]$N )) p.forest = ggplot( subset(ps.trt, treatment != trt.labels[ref.trt]), aes(x = mean, y = treatment, xmin = hdi.lo, xmax = hdi.hi) ) + geom_vline(xintercept = 0, linetype = "dashed", colour = "grey50") + geom_errorbarh(height = 0.25, linewidth = 0.8) + geom_point(size = 3, colour = "steelblue4") + labs( title = "IPD-NMA Stage 2: Pooled treatment effects vs reference", subtitle = sprintf("Multivariate RE model | %d studies | tau median = %.2f", s2.data$ns, median(post.mat[, "tau"])), x = "Mean difference (95% CrI)", y = NULL ) + theme_bw(base_size = 12) + theme(panel.grid.minor = element_blank()) print(p.forest) # 15. Tile plot: pooled moderator effects tmod.global[t, p] ------------------- gamma.grid = expand.grid(t = seq_len(s2.data$Nt), p = seq_len(Nprog), stringsAsFactors = FALSE) gamma.grid$param = paste0("tmod.global[", gamma.grid$t, ",", gamma.grid$p, "]") gamma.grid = subset(gamma.grid, param %in% colnames(post.mat) & t != ref.trt) gamma.grid$mean = colMeans(post.mat[, gamma.grid$param, drop = FALSE]) gamma.grid$q2.5 = apply(post.mat[, gamma.grid$param, drop = FALSE], 2, quantile, 0.025) gamma.grid$q97.5 = apply(post.mat[, gamma.grid$param, drop = FALSE], 2, quantile, 0.975) gamma.grid$P.pos = colMeans(post.mat[, gamma.grid$param, drop = FALSE] > 0) gamma.grid$treatment = factor(trt.labels[gamma.grid$t], levels = trt.labels) gamma.grid$covariate = factor(prog.labels[gamma.grid$p], levels = rev(prog.labels)) p.tile = ggplot(gamma.grid, aes(x = treatment, y = covariate, fill = mean)) + geom_tile(colour = "white", linewidth = 0.5) + geom_text(aes(label = sprintf("%.2f", mean)), size = 3.2) + scale_fill_gradient2(low = "steelblue4", mid = "white", high = "firebrick3", midpoint = 0, name = "Post.\nmean") + labs(title = "Stage 2: Pooled moderator effects tmod.global[t, p]", subtitle = "Differential covariate slope vs reference treatment", x = NULL, y = NULL) + theme_bw(base_size = 12) + theme(panel.grid = element_blank()) print(p.tile) # 16. Trace plots ----------------------------------------------------------- trace.candidates = c( paste0("d[", seq_len(min(4, s2.data$Nt)), "]"), paste0("tmod.global[2,", seq_len(min(2, Nprog)), "]"), "tau") trace.params = intersect(trace.candidates, colnames(post.mat)) for (p in trace.params) traceplot(samp.s2[, p], main = p, ask = FALSE) # 17. Patient-level prediction function -------------------------------------- # Computes the posterior distribution of the personalised treatment effect # for a new patient vs the reference treatment, using the Stage 2 posterior. # Patient covariates are centred using grand means from the NMA data and # scaled using the same pooled within-study SDs as Stage 1. # Arguments: # patient.covs : named list/vector with raw covariate values; # names must match prog.cols # post.mat : Stage 2 posterior matrix (from combined.samp or samp.s2) # grand.means : named vector of grand means (from dry run) # sds.ref : named vector of pooled within-study SDs (from dry run) # prog.cols : covariate names in the canonical order # Nt, ref.trt : as defined globally # Returns: data frame with one row per non-reference treatment, columns: # treatment, mean, sd, hdi.lo, median, hdi.hi, P.pos predict.patient = function(patient.covs, post.mat, grand.means, sds.ref, prog.cols, Nt, ref.trt, binary.vars = NULL, prob = 0.95) { # Centre and scale patient covariates to match Stage 1 scale. x.cs = sapply(prog.cols, function(v) { cent = patient.covs[[v]] - grand.means[[v]] s.key = paste0(v, ".c") s = sds.ref[s.key] if (v %in% binary.vars || is.na(s) || s == 0) cent else cent / s }) # Personalised draws for the reference treatment ref.d.col = paste0("d[", ref.trt, "]") ref.mod.cols = intersect( paste0("tmod.global[", ref.trt, ",", seq_along(prog.cols), "]"), colnames(post.mat) ) ref.draws = post.mat[, ref.d.col] + as.numeric(post.mat[, ref.mod.cols, drop = FALSE] %*% x.cs) # For each non-reference treatment, compute personalised effect vs ref.trt results = lapply(setdiff(seq_len(Nt), ref.trt), function(t) { d.col = paste0("d[", t, "]") mod.cols = intersect( paste0("tmod.global[", t, ",", seq_along(prog.cols), "]"), colnames(post.mat) ) pers.draws = (post.mat[, d.col] + as.numeric(post.mat[, mod.cols, drop = FALSE] %*% x.cs)) - ref.draws h = HDInterval::hdi(pers.draws, credMass = prob) data.frame( treatment = paste0("Trt", t), mean = mean(pers.draws), sd = sd(pers.draws), hdi.lo = unname(h["lower"]), median = unname(quantile(pers.draws, 0.5)), hdi.hi = unname(h["upper"]), P.pos = mean(pers.draws > 0), row.names = NULL ) }); do.call(rbind, results) } # 18. Example prediction ---------------------------------------------------- # Fill in raw covariate values for the target patient; names must match prog.cols. # Unspecified cov* variables default to their grand mean (i.e. zero centred effect). example.patient = as.list(grand.means) # start from grand means (= zero after centring) example.patient[["baseline"]] = 20 # override with patient-specific values example.patient[["age"]] = 40 example.patient[["sex"]] = 1 (pred = predict.patient( patient.covs = example.patient, post.mat = post.mat, grand.means = grand.means, sds.ref = sds.ref, prog.cols = prog.cols, Nt = Nt, ref.trt = 1, binary.vars = "sex" # names of vars to centre-only (no scaling) )) cat("\n--- Personalised treatment effects vs reference ---\n") print(as_tibble(pred))