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

# 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(Nt))
prog.labels = prog.cols


# 13. Printed posterior summaries --------------------------------------------

cat("\n--- Treatment effects vs reference (d[t]) ---\n")
d.params = paste0("d[", seq_len(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(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))

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",
                       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(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, 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))




