# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 11-ipd-component-nma-binary.R # # JAGS code for AD, IPD and combined IPD-AD CNMA models (binary) # # Binomial likelihood. Efthimiou et al., 2022, Statistics in Medicine. # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 0. Packages -------------------------------------------------------------- library(HDInterval) library(igraph) library(MCMCvis) library(netmeta) library(readxl) library(rjags) library(runjags) library(stringr) library(tidyverse) library(openxlsx) source("_utils/utils.R") # 1. Load datasets and network plot ---------------------------------------- data = openxlsx::read.xlsx("_data/_cnma/bin/test-data.xlsx", sheet=1) data.ipd = openxlsx::read.xlsx("_data/_cnma/bin/test-data-ipd.xlsx", sheet="ipd") data.ad = openxlsx::read.xlsx("_data/_cnma/bin/test-data-ipd.xlsx", sheet="ad") # Make sure all datasets are sorted alphabetically by author name data = data %>% arrange(study, arm) data.ipd = data.ipd %>% arrange(study, arm) data.ad = data.ad %>% arrange(study, arm) # Treatment-level netgraph pn = data %>% dplyr::reframe(N = sum(N), event = sum(event), .by = c(study, treat)) %>% pairwise(studlab = study, treat = treat, event = event, n = N, data = .) cairo_pdf("results/bin-netgraph.pdf", family = "Tahoma") netmeta(pn) %>% netgraph( seq = "optimal", col = "#1C1C1C", plastic = FALSE, points = TRUE, pch = 21, cex.points = 5, col.points = "#1C1C1C", bg.points = "#647980", thickness = "se.random", offset = 0.035, multiarm = FALSE, number.of.studies = TRUE) dev.off() # Component-level netgraph ig.elem = data %>% pairwise(studlab = study, treat = comb2, event = event, n = N, data = .) %>% {graph.data.frame(.[,c("treat1", "treat2")], directed = F)} cairo_pdf("results/bin-netgraph-components.pdf", height = 15, width = 15) set.seed(1234) plot(ig.elem, layout = layout_with_fr(ig.elem), vertex.label.cex = 0.8, vertex.label.family = "Tahoma", vertex.label.color = "black", vertex.label.dist = .5, vertex.label.degree = -pi/4, vertex.size = 5, vertex.color = "#647980", edge.color = "gray50") dev.off() # 2. AD-Additive CNMA (Model III) ---------------------------------------------- ## 2.1 Preprate data ----------------------------------------------------------- # Number of events per arm events = data %>% dplyr::select(study, event, arm) %>% pivot_wider(names_from = arm, values_from = event) %>% {.[,-1]} %>% as.matrix() # Sample size per arm n = data %>% dplyr::select(study, N, arm) %>% pivot_wider(names_from = arm, values_from = N) %>% {.[,-1]} %>% as.matrix() # Number of distinct studies K = data %>% distinct(study) %>% nrow() # Number of arms per study Narms = data %>% summarise(max(arm), .by = study) %>% dplyr::pull(2) # Names and number of components comp.names = paste(data$comb2, collapse = " + ") %>% strsplit(" \\+ ") %>% {unique(.[[1]])} Ncomps = length(comp.names) # Component matrices, numbered c1, c2, ... comp.matrices = map(comp.names, ~ make.components(data, comp.name = .x)) %>% {names(.) = paste0("c", 1:length(comp.names)); .} # Compile everything datlist = list(r = events, n = n, K = K, Narms = Narms, Ncomps = Ncomps) %>% append(comp.matrices) ## 2.2 Model formula ----------------------------------------------------------- M = " model { ## Loop through all studies for (i in 1:K) { w[i,1] <- 0 theta[i,1] <- 0 ## Binomial link for (k in 1:Narms[i]) { r[i,k] ~ dbin(p[i,k], n[i,k]) } ## Trial intercept (baseline arm) logit(p[i,1]) <- u[i] ## Add random treatment effect (arms 2+) for (k in 2:Narms[i]) { logit(p[i,k]) <- u[i] + theta[i,k] ## Distribution of random effects theta[i,k] ~ dnorm(md[i,k], precd[i,k]) ## Accounting for correlation in multi-arm trials ## This can either be handled using multivariate normals ## (Lu & Ades, 2004, SiM, tau/2) or conditional univariate distributions ## given all arms from 2 to k-1 (see eq. 14 in https://bit.ly/3GhgrG2) ## The condition approach is chosen here. md[i,k] <- mean[i,k] + sw[i,k] w[i,k] <- theta[i,k] - mean[i,k] sw[i,k] <- sum(w[i,1:(k-1)]) / (k - 1) precd[i,k] <- prec * 2 * (k - 1) / k ## Consistency equation mean[i,k] <- A1[i,k] - B1[i] ## Adapt number of components according (12 here) A1[i,k] <- d[1] * (1 - equals(c1[i,k], 0)) + d[2] * (1 - equals(c2[i,k], 0)) + d[3] * (1 - equals(c3[i,k], 0)) + d[4] * (1 - equals(c4[i,k], 0)) + d[5] * (1 - equals(c5[i,k], 0)) + d[6] * (1 - equals(c6[i,k], 0)) + d[7] * (1 - equals(c7[i,k], 0)) + d[8] * (1 - equals(c8[i,k], 0)) + d[9] * (1 - equals(c9[i,k], 0)) + d[10] * (1 - equals(c10[i,k], 0)) + d[11] * (1 - equals(c11[i,k], 0)) + d[12] * (1 - equals(c12[i,k], 0)) } ## Adapt number of components according (12 here) B1[i] <- d[1] * (1 - equals(c1[i,1], 0)) + d[2] * (1 - equals(c2[i,1], 0)) + d[3] * (1 - equals(c3[i,1], 0)) + d[4] * (1 - equals(c4[i,1], 0)) + d[5] * (1 - equals(c5[i,1], 0)) + d[6] * (1 - equals(c6[i,1], 0)) + d[7] * (1 - equals(c7[i,1], 0)) + d[8] * (1 - equals(c8[i,1], 0)) + d[9] * (1 - equals(c9[i,1], 0)) + d[10] * (1 - equals(c10[i,1], 0)) + d[11] * (1 - equals(c11[i,1], 0)) + d[12] * (1 - equals(c12[i,1], 0)) } ## Non-informative prior distribution baseline arm log-odds for (i in 1:K) { u[i] ~ dnorm(0, 0.01) } ## Informative prior for heterogeneity ## Log-normal derived from Turner, 2014, SiM (Table iv) ## Alternative: use LN(−2.34,1.72^2) for non-pharma vs. control prec <- 1 / tau tau ~ dlnorm(-1.67, inv.sd) inv.sd <- 1 / pow(1.472,2) ## Example comparisons (add more if needed) or.example1 <- exp(d[1] + d[3] - d[1]) # (pl + pe) VS (pl) ## prior distribution for basic parameters for (k in 1:Ncomps) { d[k] ~ dnorm(0, 0.01) } for (k in 1:Ncomps) { ORd[k] <- exp(d[k]) } #monitor# d, ORd, tau, or.example1 } " ## 2.3 Run model --------------------------------------------------------------- res = run.jags(M, data = datlist, summarise = FALSE, n.chains = 4, burnin = 1000, sample = 20000) pdf("results/bin-diagnostics-model-iii.pdf", width = 10) plot(res); dev.off() # Get samples mcmc = coda::as.mcmc(res) # 3. AD-Spike-Slab CNMA (Model V) ---------------------------------------------- ## 3.1 Prepare data ------------------------------------------------------------ # Additional data objects required # Number of interactions (total) datlist$Ninter = calc.interaction.number(Ncomps) datlist$interactions = make.interactions(Ns = K, Nc = Ncomps, na = Narms, component_matrices = comp.matrices, Ninter = datlist$Ninter) ## 3.2 Model formula ----------------------------------------------------------- M = "model { for (i in 1:K) { w[i,1] <- 0 theta[i,1] <- 0 ## Likelihood for (k in 1:Narms[i]) { r[i,k] ~ dbin(p[i,k], n[i,k]) } logit(p[i,1]) <- u[i] for (k in 2:Narms[i]) { logit(p[i,k]) <- u[i] + theta[i,k] ## Distribution of random effects theta[i,k] ~ dnorm(md[i,k], precd[i,k]) ## Accounting for correlation in multi-arm trials md[i,k] <- mean[i,k] + sw[i,k] w[i,k] <- (theta[i,k] - mean[i,k]) sw[i,k] <- sum(w[i,1:(k-1)]) / (k - 1) precd[i,k] <- prec * 2 * (k - 1) / k ## Consistency equations mean[i,k] <- A1[i,k] - B1[i] A1[i,k] <- d[1]*(1 - equals(c1[i,k], 0)) + d[2]*(1 - equals(c2[i,k], 0)) + d[3]*(1 - equals(c3[i,k], 0)) + d[4]*(1 - equals(c4[i,k], 0)) + d[5]*(1 - equals(c5[i,k], 0)) + d[6]*(1 - equals(c6[i,k], 0)) + d[7]*(1 - equals(c7[i,k], 0)) + d[8]*(1 - equals(c8[i,k], 0)) + d[9]*(1 - equals(c9[i,k], 0)) + d[10]*(1 - equals(c10[i,k], 0)) + d[11]*(1 - equals(c11[i,k], 0)) + d[12]*(1 - equals(c12[i,k], 0)) + inprod(gamma[], interactions[i,k,]) } B1[i] <- d[1]*(1 - equals(c1[i,1], 0)) + d[2]*(1 - equals(c2[i,1], 0)) + d[3]*(1 - equals(c3[i,1], 0)) + d[4]*(1 - equals(c4[i,1], 0)) + d[5]*(1 - equals(c5[i,1], 0)) + d[6]*(1 - equals(c6[i,1], 0)) + d[7]*(1 - equals(c7[i,1], 0)) + d[8]*(1 - equals(c8[i,1], 0)) + d[9]*(1 - equals(c9[i,1], 0)) + d[10]*(1 - equals(c10[i,1], 0)) + d[11]*(1 - equals(c11[i,1], 0)) + d[12]*(1 - equals(c12[i,1], 0)) + inprod(gamma[], interactions[i,1,]) } ## Priors for baseline log-odds for (i in 1:K) { u[i] ~ dnorm(0, 0.01) } ## Informative prior for heterogeneity ## Log-normal derived from Turner, 2014, SiM (Table iv) ## Alternative: use LN(−2.34,1.72^2) for non-pharma vs. control prec <- 1 / tau tau ~ dlnorm(-1.67, inv.sd) inv.sd <- 1 / (1.472 * 1.472) ## Stochastic Search Variable Selection (SSVS) ## We use a spike-and-slab prior here instead of Bayesian LASSO (Laplace prior) for (k in 1:Ninter) { IndA[k] ~ dcat(Pind[]) Ind[k] <- IndA[k] - 1 gamma[k] ~ dnorm(0, tauCov[IndA[k]]) } zeta <- pow(eta, -2) eta ~ dnorm(0, 1000) I(0,) tauCov[1] <- zeta tauCov[2] <- zeta * 0.01 # g = 100 Pind[1] <- 0.5 # P(I_j = 1) = 0.5 Pind[2] <- 0.5 ## Priors for basic parameters for (k in 1:Ncomps) { d[k] ~ dnorm(0, 0.01) ORd[k] <- exp(d[k]) } ## Example comparisons ## (pl + ftf + pe + ps + ive) VS (wl) ## Run all.interactions to retrieve all required gammas ## e.g. all.interactions(c(1:4,6), Nc = 12, components = comp.names) or.example1 <- exp( d[1] + d[2] + d[3] + d[4] + d[6] + gamma[1] + gamma[2] + gamma[3] + gamma[5] + gamma[12] + gamma[13] + gamma[15] + gamma[22] + gamma[24] + gamma[32] - d[10] ) ## Ind gives the selection percentage (indicator variable) in mean column #monitor# tau, d, ORd, gamma, Ind, or.example1, eta }" ## 3.3 Run model --------------------------------------------------------------- res = run.jags(M, data = datlist, summarise = FALSE, n.chains = 4, burnin = 1000, sample = 20000) pdf("results/bin-diagnostics-model-v.pdf", width = 10) plot(res); dev.off() # Get samples mcmc = coda::as.mcmc(res) # Plot coefficient against selection frequency X = mcmc %>% as.data.frame() %>% dplyr::select(`gamma[1]`:`gamma[66]`) %>% apply(2, median) Y = mcmc %>% as.data.frame() %>% dplyr::select(`Ind[1]`:`Ind[66]`) %>% apply(2, mean) pdf("results/bin-components-ssvs-m5.pdf") plot(abs(X),Y, pch = 19, xlab = "Coefficient (abs)", ylab = "Selection Frequency") text(abs(X),Y, all.interactions(1:12, 12, comp.names)$label, pos = 2) dev.off() # 4. AD-Spike-Slab CNMA (Model V, include expert opinion) ---------------------- ## 4.1 Prepare data ------------------------------------------------------------ # (Not needed, see previous sections) ## 4.2 Model formula ----------------------------------------------------------- M = " model { for (i in 1:K) { w[i,1] <- 0 theta[i,1] <- 0 for (k in 1:Narms[i]) { r[i,k] ~ dbin(p[i,k], n[i,k]) } logit(p[i,1]) <- u[i] for (k in 2:Narms[i]) { logit(p[i,k]) <- u[i] + theta[i,k] # Distribution of random effects theta[i,k] ~ dnorm(md[i,k], precd[i,k]) md[i,k] <- mean[i,k] + sw[i,k] w[i,k] <- theta[i,k] - mean[i,k] sw[i,k] <- sum(w[i,1:(k-1)]) / (k - 1) precd[i,k] <- prec * 2 * (k - 1) / k # Consistency equations mean[i,k] <- A1[i,k] - B1[i] A1[i,k] <- d[1] * (1 - equals(c1[i,k], 0)) + d[2] * (1 - equals(c2[i,k], 0)) + d[3] * (1 - equals(c3[i,k], 0)) + d[4] * (1 - equals(c4[i,k], 0)) + d[5] * (1 - equals(c5[i,k], 0)) + d[6] * (1 - equals(c6[i,k], 0)) + d[7] * (1 - equals(c7[i,k], 0)) + d[8] * (1 - equals(c8[i,k], 0)) + d[9] * (1 - equals(c9[i,k], 0)) + d[10]* (1 - equals(c10[i,k], 0)) + d[11]* (1 - equals(c11[i,k], 0)) + d[12]* (1 - equals(c12[i,k], 0)) + inprod(gamma[], interactions[i,k,]) } B1[i] <- d[1] * (1 - equals(c1[i,1], 0)) + d[2] * (1 - equals(c2[i,1], 0)) + d[3] * (1 - equals(c3[i,1], 0)) + d[4] * (1 - equals(c4[i,1], 0)) + d[5] * (1 - equals(c5[i,1], 0)) + d[6] * (1 - equals(c6[i,1], 0)) + d[7] * (1 - equals(c7[i,1], 0)) + d[8] * (1 - equals(c8[i,1], 0)) + d[9] * (1 - equals(c9[i,1], 0)) + d[10]* (1 - equals(c10[i,1], 0)) + d[11]* (1 - equals(c11[i,1], 0)) + d[12]* (1 - equals(c12[i,1], 0)) + inprod(gamma[], interactions[i,1,]) } ## Prior distribution for log-odds in baseline arm for (i in 1:K) { u[i] ~ dnorm(0, 0.01) } ## Prior distribution for heterogeneity prec <- 1 / tau tau ~ dlnorm(-1.67, inv.sd) inv.sd <- 1 / (1.472^2) ## SSVS: interaction terms without prior info for (k in c(1:26, 28, 30:34, 36:46, 49:58, 60:66)) { IndA[k] ~ dcat(Pind[]) Ind[k] <- IndA[k] - 1 gamma[k] ~ dnorm(0, tauCov[IndA[k]]) } zeta <- pow(eta, -2) eta ~ dnorm(0, 1000) I(0,) tauCov[1] <- zeta tauCov[2] <- zeta * 0.01 # g = 100 Pind[1] <- 0.5 # P(I_j = 1) = 0.5 Pind[2] <- 0.5 ## SSVS: interaction terms with prior info for (k in c(27, 29, 35, 47, 48, 59)) { IndA[k] ~ dcat(Pind1[]) Ind[k] <- IndA[k] - 1 gamma[k] ~ dnorm(0, tauCov[IndA[k]]) } Pind1[1] <- 0.2 Pind1[2] <- 0.8 ## Prior distribution for basic parameters for (k in 1:Ncomps) { d[k] ~ dnorm(0, 0.01) ORd[k] <- exp(d[k]) } ## Example comparisons ## (pl + ftf + pe + ps + ive) VS (wl) ## Run all.interactions to retrieve all required gammas ## e.g. all.interactions(c(1:4,6), Nc = 12, components = comp.names) or.example1 <- exp( d[1] + d[2] + d[3] + d[4] + d[6] + gamma[1] + gamma[2] + gamma[3] + gamma[5] + gamma[12] + gamma[13] + gamma[15] + gamma[22] + gamma[24] + gamma[32] - d[10] ) ## Ind gives the selection percentage (indicator variable) in mean column #monitor# tau, d, ORd, gamma, Ind, or.example1, eta } " ## 4.3 Run model --------------------------------------------------------------- res = run.jags(M, data = datlist, summarise = FALSE, n.chains = 4, burnin = 1000, sample = 20000) pdf("results/bin-diagnostics-model-v-expert.pdf", width = 10) plot(res); dev.off() # Get samples mcmc = coda::as.mcmc(res) # Plot coefficient against selection frequency X = mcmc %>% as.data.frame() %>% dplyr::select(`gamma[1]`:`gamma[66]`) %>% apply(2, median) Y = mcmc %>% as.data.frame() %>% dplyr::select(`Ind[1]`:`Ind[66]`) %>% apply(2, mean) pdf("results/bin-components-ssvs-m5-expert.pdf") plot(abs(X),Y, pch = 19, xlab = "Coefficient (abs)", ylab = "Selection Frequency") text(abs(X),Y, all.interactions(1:12, 12, comp.names)$label, pos = 2) dev.off() # 5. IPD-AD Additive CNMA (Model VIII) ----------------------------------------- ## 5.1 Data Preparation -------------------------------------------------------- ### 5.1.1 IPD Elements --------------------------------------------------------- datlist.ipd = data.ipd %>% { list(Npatients = nrow(.), studyIPD = .$study %>% factor() %>% as.integer(), armIPD = .$arm, y.ipd = .$event, baseline = .$baseline, age = .$age, gender = .$gender, naIPD = summarise(., max(arm), .by = study) %>% pull(2), Ns.IPD = .$study %>% unique() %>% length) } data.ipd.agg = data.ipd %>% distinct(study, arm, .keep_all = TRUE) comp.names = data.ipd.agg %>% {paste(.$comb2, collapse = " + ")} %>% strsplit(" \\+ ") %>% {unique(.[[1]])} Ncomps = length(comp.names) comp.matrices = map(comp.names, ~ make.components(data.ipd.agg, comp.name = .x)) %>% {names(.) = paste0("cIPD_", 1:length(comp.names)); .} datlist.ipd = datlist.ipd %>% append(comp.matrices) ### 5.1.2 AD Elements ---------------------------------------------------------- # Check if AD components are all included in IPD paste(data.ad$comb2, collapse = " + ") %>% strsplit(" \\+ ") %>% {all(unique(.[[1]]) %in% comp.names)} # Check which AD components are included (numbers) mask = paste(data.ad$comb2, collapse = " + ") %>% strsplit(" \\+ ") %>% {comp.names %in% unique(.[[1]])} components.to.include = which(mask) datlist.ad = data.ad %>% { list(Ns = .$study %>% unique() %>% length(), r = dplyr::select(., study, event, arm) %>% pivot_wider(names_from = arm, values_from = event) %>% {.[,-1]} %>% as.matrix(), n = dplyr::select(., study, N, arm) %>% pivot_wider(names_from = arm, values_from = N) %>% {.[,-1]} %>% as.matrix(), na = summarise(., max(arm), .by = study) %>% dplyr::pull(2)) } comp.matrices = map(comp.names, ~ make.components(data.ad, comp.name = .x)) %>% {names(.) = paste0("c", 1:length(comp.names)); .} datlist.ad = datlist.ad %>% append(comp.matrices) ### 5.1.3 Combine -------------------------------------------------------------- datlist = c(datlist.ipd, datlist.ad) datlist$Nc = length(comp.names) ## 5.2 Model formula ----------------------------------------------------------- M = " model { ##### IPD studies ##### for (i in 1:Npatients) { ## Outcome model for binary data y.ipd[i] ~ dbern(p[i]) logit(p[i]) <- muIPD[i] ## phi combines intercept and random Tx effect (alpha/u + delta) muIPD[i] <- phi.IPD[studyIPD[i], armIPD[i]] + ## Prognostic coeffients beta.baseline * baseline[i] + beta.age * age[i] + beta.gender * gender[i] + ## Component-covariate interactions (baseline, age, gender) gamma.base[1] * cIPD_1[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[3] * cIPD_3[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[4] * cIPD_4[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[5] * cIPD_5[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[6] * cIPD_6[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[7] * cIPD_7[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[8] * cIPD_8[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[9] * cIPD_9[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[10] * cIPD_10[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[11] * cIPD_11[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[12] * cIPD_12[studyIPD[i], armIPD[i]] * baseline[i] ## Additional component-covariate interaction can be added if ## needed, but increases sampling time. ## gamma.age[1] * cIPD_1[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[3] * cIPD_3[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[4] * cIPD_4[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[5] * cIPD_5[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[6] * cIPD_6[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[7] * cIPD_7[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[8] * cIPD_8[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[9] * cIPD_9[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[10] * cIPD_10[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[11] * cIPD_11[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[12] * cIPD_12[studyIPD[i], armIPD[i]] * age[i] + ## gamma.gender[1] * cIPD_1[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[3] * cIPD_3[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[4] * cIPD_4[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[5] * cIPD_5[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[6] * cIPD_6[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[7] * cIPD_7[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[8] * cIPD_8[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[9] * cIPD_9[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[10] * cIPD_10[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[11] * cIPD_11[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[12] * cIPD_12[studyIPD[i], armIPD[i]] * gender[i] } for (i in 1:Ns.IPD) { w1[i,1] <- 0 delta1[i,1] <- 0 ## Define phi as intercept + random Tx effect for (k in 1:naIPD[i]) { phi.IPD[i,k] <- uIPD[i] + delta1[i,k] } ## Define random Tx effect (allow for multi-arm trials) ## Accounting for correlation in multi-arm trials ## This can either be handled using multivariate normals ## (Lu & Ades, 2004, SiM, tau/2) or conditional univariate distributions ## given all arms from 2 to k-1 (see eq. 14 in https://bit.ly/3GhgrG2) ## The conditional approach is chosen here. for (k in 2:naIPD[i]) { delta1[i,k] ~ dnorm(md1[i,k], precd1[i,k]) md1[i,k] <- mean1[i,k] + sw1[i,k] precd1[i,k] <- 2 * prect * (k - 1) / k w1[i,k] <- (delta1[i,k] - mean1[i,k]) sw1[i,k] <- sum(w1[i,1:(k-1)]) / (k - 1) ## consistency equations mean1[i,k] <- AIPD[i,k] - BIPD[i] AIPD[i,k] <- d[1]*(1 - equals(cIPD_1[i,k],0)) + d[2]*(1 - equals(cIPD_2[i,1],0)) + d[3]*(1 - equals(cIPD_3[i,k],0)) + d[4]*(1 - equals(cIPD_4[i,k],0)) + d[5]*(1 - equals(cIPD_5[i,k],0)) + d[6]*(1 - equals(cIPD_6[i,k],0)) + d[7]*(1 - equals(cIPD_7[i,k],0)) + d[8]*(1 - equals(cIPD_8[i,k],0)) + d[9]*(1 - equals(cIPD_9[i,k],0)) + d[10]*(1 - equals(cIPD_10[i,k],0)) + d[11]*(1 - equals(cIPD_11[i,k],0)) + d[12]*(1 - equals(cIPD_12[i,k],0)) } BIPD[i] <- d[1]*(1 - equals(cIPD_1[i,1],0)) + d[2]*(1 - equals(cIPD_2[i,1],0)) + d[3]*(1 - equals(cIPD_3[i,1],0)) + d[4]*(1 - equals(cIPD_4[i,1],0)) + d[5]*(1 - equals(cIPD_5[i,1],0)) + d[6]*(1 - equals(cIPD_6[i,1],0)) + d[7]*(1 - equals(cIPD_7[i,1],0)) + d[8]*(1 - equals(cIPD_8[i,1],0)) + d[9]*(1 - equals(cIPD_9[i,1],0)) + d[10]*(1 - equals(cIPD_10[i,1],0)) + d[11]*(1 - equals(cIPD_11[i,1],0)) + d[12]*(1 - equals(cIPD_12[i,1],0)) } #### AD studies #### for (i in 1:Ns) { w[i,1] <- 0 delta[i,1] <- 0 for (k in 1:na[i]) { ## Binomial likelihood: number of events in arm k of study i r[i,k] ~ dbin(p.ad[i,k], n[i,k]) logit(p.ad[i,k]) <- phi[i,k] phi[i,k] <- u[i] + delta[i,k] } for (k in 2:na[i]) { delta[i,k] ~ dnorm(md[i,k], precd[i,k]) md[i,k] <- mean[i,k] + sw[i,k] precd[i,k] <- 2 * prect * (k - 1) / k w[i,k] <- delta[i,k] - mean[i,k] sw[i,k] <- sum(w[i,1:(k-1)]) / (k - 1) ## consistency equations mean[i,k] <- A1[i,k] - B1[i] A1[i,k] <- d[1] * (1 - equals(c1[i,k], 0)) + d[2] * (1 - equals(c2[i,k], 0)) + d[3] * (1 - equals(c3[i,k], 0)) + d[4] * (1 - equals(c4[i,k], 0)) + d[5] * (1 - equals(c5[i,k], 0)) + d[6] * (1 - equals(c6[i,k], 0)) + d[7] * (1 - equals(c7[i,k], 0)) + d[8] * (1 - equals(c8[i,k], 0)) + d[9] * (1 - equals(c9[i,k], 0)) + d[10] * (1 - equals(c10[i,k], 0)) + d[11] * (1 - equals(c11[i,k], 0)) + d[12] * (1 - equals(c12[i,k], 0)) } B1[i] <- d[1] * (1 - equals(c1[i,1], 0)) + d[2] * (1 - equals(c2[i,1], 0)) + d[3] * (1 - equals(c3[i,1], 0)) + d[4] * (1 - equals(c4[i,1], 0)) + d[5] * (1 - equals(c5[i,1], 0)) + d[6] * (1 - equals(c6[i,1], 0)) + d[7] * (1 - equals(c7[i,1], 0)) + d[8] * (1 - equals(c8[i,1], 0)) + d[9] * (1 - equals(c9[i,1], 0)) + d[10] * (1 - equals(c10[i,1], 0)) + d[11] * (1 - equals(c11[i,1], 0)) + d[12] * (1 - equals(c12[i,1], 0)) } ## prior distribution for baseline arm of study i for (i in 1:Ns) { u[i] ~ dnorm(0, 0.001) } for (i in 1:Ns.IPD) { uIPD[i] ~ dnorm(0, 0.001) } ## prior distribution for heterogeneity tau ~ dnorm(0, 0.1) I(0, ) prect <- 1 / tau.sq tau.sq <- pow(tau, 2) ## prior distribution for basic parameters for (k in 1:Nc) { d[k] ~ dnorm(0, 0.001) ORd[k] <- exp(d[k]) } ## prior distribution for penalization of gamma ## (can be removed) lambda ~ dgamma(1, 0.1) for (k in 1:Nc) { gamma.base[k] ~ ddexp(0, lambda) ## gamma.age[k] ~ ddexp(0, lambda) ## gamma.gender[k] ~ ddexp(0, lambda) } beta.baseline ~ dnorm(0, 1) beta.age ~ dnorm(0, 1) beta.gender ~ dnorm(0, 1) #monitor# d, ORd, tau, gamma.base } " ## 5.3 Run model --------------------------------------------------------------- res = run.jags(M, data = datlist, summarise = FALSE, n.chains = 4, burnin = 500, sample = 1000) # Get samples mcmc = coda::as.mcmc(res) # Get summary sum.res = summary(res) # 6. IPD-AD Additive CNMA (Model VIII, bias parameter) ------------------------- ## 6.1 Data Preparation -------------------------------------------------------- # (not necessary, see 5.1) ## 6.2 Model formula ----------------------------------------------------------- M = " model { ##### IPD studies ##### for (i in 1:Npatients) { ## Outcome model for binary data y.ipd[i] ~ dbern(p[i]) logit(p[i]) <- muIPD[i] ## phi combines intercept and random Tx effect (alpha/u + delta) muIPD[i] <- phi.IPD[studyIPD[i], armIPD[i]] + ## Prognostic coeffients beta.baseline * baseline[i] + beta.age * age[i] + beta.gender * gender[i] + ## Component-covariate interactions (baseline, age, gender) gamma.base[1] * cIPD_1[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[3] * cIPD_3[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[4] * cIPD_4[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[5] * cIPD_5[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[6] * cIPD_6[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[7] * cIPD_7[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[8] * cIPD_8[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[9] * cIPD_9[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[10] * cIPD_10[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[11] * cIPD_11[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[12] * cIPD_12[studyIPD[i], armIPD[i]] * baseline[i] ## Additional component-covariate interaction can be added if ## needed, but increases sampling time. ## gamma.age[1] * cIPD_1[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[3] * cIPD_3[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[4] * cIPD_4[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[5] * cIPD_5[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[6] * cIPD_6[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[7] * cIPD_7[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[8] * cIPD_8[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[9] * cIPD_9[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[10] * cIPD_10[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[11] * cIPD_11[studyIPD[i], armIPD[i]] * age[i] + ## gamma.age[12] * cIPD_12[studyIPD[i], armIPD[i]] * age[i] + ## gamma.gender[1] * cIPD_1[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[3] * cIPD_3[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[4] * cIPD_4[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[5] * cIPD_5[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[6] * cIPD_6[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[7] * cIPD_7[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[8] * cIPD_8[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[9] * cIPD_9[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[10] * cIPD_10[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[11] * cIPD_11[studyIPD[i], armIPD[i]] * gender[i] + ## gamma.gender[12] * cIPD_12[studyIPD[i], armIPD[i]] * gender[i] } for (i in 1:Ns.IPD) { w1[i,1] <- 0 delta1[i,1] <- 0 ## Define phi as intercept + random Tx effect for (k in 1:naIPD[i]) { phi.IPD[i,k] <- uIPD[i] + delta1[i,k] } ## Define random Tx effect (allow for multi-arm trials) ## Accounting for correlation in multi-arm trials ## This can either be handled using multivariate normals ## (Lu & Ades, 2004, SiM, tau/2) or conditional univariate distributions ## given all arms from 2 to k-1 (see eq. 14 in https://bit.ly/3GhgrG2) ## The conditional approach is chosen here. for (k in 2:naIPD[i]) { delta1[i,k] ~ dnorm(md1[i,k], precd1[i,k]) md1[i,k] <- mean1[i,k] + sw1[i,k] precd1[i,k] <- 2 * prect * (k - 1) / k w1[i,k] <- (delta1[i,k] - mean1[i,k]) sw1[i,k] <- sum(w1[i,1:(k-1)]) / (k - 1) ## consistency equations mean1[i,k] <- AIPD[i,k] - BIPD[i] AIPD[i,k] <- d[1]*(1 - equals(cIPD_1[i,k],0)) + d[2]*(1 - equals(cIPD_2[i,1],0)) + d[3]*(1 - equals(cIPD_3[i,k],0)) + d[4]*(1 - equals(cIPD_4[i,k],0)) + d[5]*(1 - equals(cIPD_5[i,k],0)) + d[6]*(1 - equals(cIPD_6[i,k],0)) + d[7]*(1 - equals(cIPD_7[i,k],0)) + d[8]*(1 - equals(cIPD_8[i,k],0)) + d[9]*(1 - equals(cIPD_9[i,k],0)) + d[10]*(1 - equals(cIPD_10[i,k],0)) + d[11]*(1 - equals(cIPD_11[i,k],0)) + d[12]*(1 - equals(cIPD_12[i,k],0)) } BIPD[i] <- d[1]*(1 - equals(cIPD_1[i,1],0)) + d[2]*(1 - equals(cIPD_2[i,1],0)) + d[3]*(1 - equals(cIPD_3[i,1],0)) + d[4]*(1 - equals(cIPD_4[i,1],0)) + d[5]*(1 - equals(cIPD_5[i,1],0)) + d[6]*(1 - equals(cIPD_6[i,1],0)) + d[7]*(1 - equals(cIPD_7[i,1],0)) + d[8]*(1 - equals(cIPD_8[i,1],0)) + d[9]*(1 - equals(cIPD_9[i,1],0)) + d[10]*(1 - equals(cIPD_10[i,1],0)) + d[11]*(1 - equals(cIPD_11[i,1],0)) + d[12]*(1 - equals(cIPD_12[i,1],0)) } #### AD studies #### for (i in 1:Ns) { w[i,1] <- 0 delta[i,1] <- 0 for (k in 1:na[i]) { ## Binomial likelihood: number of events in arm k of study i r[i,k] ~ dbin(p.ad[i,k], n[i,k]) logit(p.ad[i,k]) <- phi[i,k] phi[i,k] <- u[i] + delta[i,k] + bias[i] } for (k in 2:na[i]) { delta[i,k] ~ dnorm(md[i,k], precd[i,k]) md[i,k] <- mean[i,k] + sw[i,k] precd[i,k] <- 2 * prect * (k - 1) / k w[i,k] <- delta[i,k] - mean[i,k] sw[i,k] <- sum(w[i,1:(k-1)]) / (k - 1) ## consistency equations mean[i,k] <- A1[i,k] - B1[i] A1[i,k] <- d[1] * (1 - equals(c1[i,k], 0)) + d[2] * (1 - equals(c2[i,k], 0)) + d[3] * (1 - equals(c3[i,k], 0)) + d[4] * (1 - equals(c4[i,k], 0)) + d[5] * (1 - equals(c5[i,k], 0)) + d[6] * (1 - equals(c6[i,k], 0)) + d[7] * (1 - equals(c7[i,k], 0)) + d[8] * (1 - equals(c8[i,k], 0)) + d[9] * (1 - equals(c9[i,k], 0)) + d[10] * (1 - equals(c10[i,k], 0)) + d[11] * (1 - equals(c11[i,k], 0)) + d[12] * (1 - equals(c12[i,k], 0)) } B1[i] <- d[1] * (1 - equals(c1[i,1], 0)) + d[2] * (1 - equals(c2[i,1], 0)) + d[3] * (1 - equals(c3[i,1], 0)) + d[4] * (1 - equals(c4[i,1], 0)) + d[5] * (1 - equals(c5[i,1], 0)) + d[6] * (1 - equals(c6[i,1], 0)) + d[7] * (1 - equals(c7[i,1], 0)) + d[8] * (1 - equals(c8[i,1], 0)) + d[9] * (1 - equals(c9[i,1], 0)) + d[10] * (1 - equals(c10[i,1], 0)) + d[11] * (1 - equals(c11[i,1], 0)) + d[12] * (1 - equals(c12[i,1], 0)) } ## prior distribution for baseline arm of study i for (i in 1:Ns) { u[i] ~ dnorm(0, 0.001) } for (i in 1:Ns.IPD) { uIPD[i] ~ dnorm(0, 0.001) } ## Bias term for unadjusted AD studies for (i in 1:Ns) { bias[i] ~ dnorm(mu.bias, prec.bias) } mu.bias ~ dnorm(0, 0.001) prec.bias <- pow(sd.bias, -2) sd.bias ~ dunif(0, 5) ## prior distribution for heterogeneity tau ~ dnorm(0, 0.1) I(0, ) prect <- 1 / tau.sq tau.sq <- pow(tau, 2) ## prior distribution for basic parameters for (k in 1:Nc) { d[k] ~ dnorm(0, 0.001) ORd[k] <- exp(d[k]) } ## prior distribution for penalization of gamma ## (can be removed) lambda ~ dgamma(1, 0.1) for (k in 1:Nc) { gamma.base[k] ~ ddexp(0, lambda) ## gamma.age[k] ~ ddexp(0, lambda) ## gamma.gender[k] ~ ddexp(0, lambda) } beta.baseline ~ dnorm(0, 1) beta.age ~ dnorm(0, 1) beta.gender ~ dnorm(0, 1) #monitor# d, ORd, tau, bias, gamma.base } " ## 6.3 Run model --------------------------------------------------------------- res = run.jags(M, data = datlist, summarise = FALSE, n.chains = 4, burnin = 500, sample = 1000) # Get samples mcmc = coda::as.mcmc(res) # Get summary sum.res = summary(res) # 7. IPD Additive CNMA (Model VIII) -------------------------------------------- ## 7.1 Data Preparation -------------------------------------------------------- # (not necessary, see 5.1) ## 7.2 Model formula ----------------------------------------------------------- M = " model { ##### IPD studies ##### for (i in 1:Npatients) { ## Outcome model for binary data y.ipd[i] ~ dbern(p[i]) logit(p[i]) <- muIPD[i] ## phi combines intercept and random Tx effect (alpha/u + delta) muIPD[i] <- phi.IPD[studyIPD[i], armIPD[i]] + ## Prognostic coefficients beta.baseline * baseline[i] + beta.age * age[i] + beta.gender * gender[i] + ## Component-covariate interactions (baseline) gamma.base[1] * cIPD_1[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[3] * cIPD_3[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[4] * cIPD_4[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[5] * cIPD_5[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[6] * cIPD_6[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[7] * cIPD_7[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[8] * cIPD_8[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[9] * cIPD_9[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[10] * cIPD_10[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[11] * cIPD_11[studyIPD[i], armIPD[i]] * baseline[i] + gamma.base[12] * cIPD_12[studyIPD[i], armIPD[i]] * baseline[i] } for (i in 1:Ns.IPD) { w1[i,1] <- 0 delta1[i,1] <- 0 ## Define phi as intercept + random Tx effect for (k in 1:naIPD[i]) { phi.IPD[i,k] <- uIPD[i] + delta1[i,k] } ## Define random Tx effect (allow for multi-arm trials) for (k in 2:naIPD[i]) { delta1[i,k] ~ dnorm(md1[i,k], precd1[i,k]) md1[i,k] <- mean1[i,k] + sw1[i,k] precd1[i,k] <- 2 * prect * (k - 1) / k w1[i,k] <- (delta1[i,k] - mean1[i,k]) sw1[i,k] <- sum(w1[i,1:(k-1)]) / (k - 1) ## consistency equations mean1[i,k] <- AIPD[i,k] - BIPD[i] AIPD[i,k] <- d[1]*(1 - equals(cIPD_1[i,k],0)) + d[2]*(1 - equals(cIPD_2[i,k],0)) + d[3]*(1 - equals(cIPD_3[i,k],0)) + d[4]*(1 - equals(cIPD_4[i,k],0)) + d[5]*(1 - equals(cIPD_5[i,k],0)) + d[6]*(1 - equals(cIPD_6[i,k],0)) + d[7]*(1 - equals(cIPD_7[i,k],0)) + d[8]*(1 - equals(cIPD_8[i,k],0)) + d[9]*(1 - equals(cIPD_9[i,k],0)) + d[10]*(1 - equals(cIPD_10[i,k],0)) + d[11]*(1 - equals(cIPD_11[i,k],0)) + d[12]*(1 - equals(cIPD_12[i,k],0)) } BIPD[i] <- d[1]*(1 - equals(cIPD_1[i,1],0)) + d[2]*(1 - equals(cIPD_2[i,1],0)) + d[3]*(1 - equals(cIPD_3[i,1],0)) + d[4]*(1 - equals(cIPD_4[i,1],0)) + d[5]*(1 - equals(cIPD_5[i,1],0)) + d[6]*(1 - equals(cIPD_6[i,1],0)) + d[7]*(1 - equals(cIPD_7[i,1],0)) + d[8]*(1 - equals(cIPD_8[i,1],0)) + d[9]*(1 - equals(cIPD_9[i,1],0)) + d[10]*(1 - equals(cIPD_10[i,1],0)) + d[11]*(1 - equals(cIPD_11[i,1],0)) + d[12]*(1 - equals(cIPD_12[i,1],0)) } ## Prior distribution for baseline arm of study i for (i in 1:Ns.IPD) { uIPD[i] ~ dnorm(0, 0.001) } ## Prior distribution for heterogeneity tau ~ dnorm(0, 0.1) I(0, ) prect <- 1 / tau.sq tau.sq <- pow(tau, 2) ## Prior distribution for basic parameters for (k in 1:Nc) { d[k] ~ dnorm(0, 0.001) ORd[k] <- exp(d[k]) } ## Prior distribution for penalization of gamma lambda ~ dgamma(1, 0.1) for (k in 1:Nc) { gamma.base[k] ~ ddexp(0, lambda) } beta.baseline ~ dnorm(0, 1) beta.age ~ dnorm(0, 1) beta.gender ~ dnorm(0, 1) #monitor# d, ORd, tau, gamma.base } " ## 7.3 Run model --------------------------------------------------------------- res = run.jags(M, data = datlist, summarise = FALSE, n.chains = 4, burnin = 500, sample = 1000) # Get samples mcmc = coda::as.mcmc(res) # Get summary sum.res = summary(res)