# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#                                                                             #
#   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/bin/_cnma/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)













