chapter_16.qmd

---
title: "Prediction of individual treatment effect using data from multiple studies"
authors:   
  - name: Orestis Efthimiou
    orcid: 0000-0002-0955-7572
    affiliations:
      - ref: ubern
affiliations:
  - id: smartdas
    name: Smart Data Analysis and Statistics B.V.
    city: Utrecht
  - id: ubern
    name: Institute of Social and Preventive Medicine (ISPM)
    city: Bern, Switzerland
format:
  html:
    toc: true
    number-sections: true
execute:
  cache: true
bibliography: 'https://api.citedrive.com/bib/0d25b38b-db8f-43c4-b934-f4e2f3bd655a/references.bib?x=eyJpZCI6ICIwZDI1YjM4Yi1kYjhmLTQzYzQtYjkzNC1mNGUyZjNiZDY1NWEiLCAidXNlciI6ICIyNTA2IiwgInNpZ25hdHVyZSI6ICI0MGFkYjZhMzYyYWE5Y2U0MjQ2NWE2ZTQzNjlhMWY3NTk5MzhhNzUxZDNjYWIxNDlmYjM4NDgwOTYzMzY5YzFlIn0=/bibliography.bib'
---

``` {r}
#| message: false
#| warning: false
#| echo: false

# List of required packages
required_packages <- c("dplyr", "table1", "kableExtra", "bipd")

# Install required packages
for (pkg in required_packages) {
    if (!requireNamespace(pkg, quietly = TRUE)) {
        install.packages(pkg)
    }
}
```

In this chapter, we discuss statistical methods for developing models to predict patient-level treatment effects using data from multiple randomized and non-randomized studies. We will first present prediction models that assume a constant treatment effect and discuss how to address heterogeneity in baseline risk. Subsequently, we will discuss approaches that allow for treatment effect modification by adopting two different approaches in an IPD-MA context, namely the risk modelling and the effect modelling approach. For both approaches, we will first discuss how to combine IPD from RCTs comparing the same two treatments. We will then discuss how these methods can be extended to include randomized data from multiple treatments, real-world data, and published aggregate data. We will discuss statistical software to implement these approaches and provide example code as supporting information. Real examples will be used throughout to illustrate the main methods.

## Estimating heterogeneous treatment effects in pairwise meta-analysis

```{r}
#| echo: false
#| message: false
library(table1)
library(dplyr)
library(kableExtra)
```

We hereby provide code for estimating patient-level treatment effects for the case when we have patient-level data from multiple randomized trials.

### Example of a continuous outcome

#### Setup
We  start by simulating an artificial dataset using the R package **bipd**: 

```{r ds}
library(bipd)
ds <- generate_ipdma_example(type = "continuous")
```

Let us have a look at the dataset:
```{r ds2}
head(ds)
```

The simulated dataset contains information on the following variables:

- the trial indicator `studyid`
- the treatment indicator `treat`, which takes the values 0 for control and 1 for active treatment
- two prognostic variables `z1` and `z2`
- the continuous outcome `y`

```{r}
#| label: tbl-summary-continuous_outcome-data
#| tbl-cap: The simulated dataset with a continuous outcome
#| echo: false
#| warning: false
dstbl <- ds %>% mutate(treat = as.factor(treat),
                       studyid = as.factor(studyid))
table1(~ z1+z2+studyid | treat, data = dstbl)
```

#### Model fitting
We synthesize the evidence using a Bayesian random effects meta-analysis model. The model is given in Equation 16.7 of the book. First we need set up the data and create the model:
```{r MAcont1}
ipd <- with(ds, ipdma.model.onestage(y = y, study = studyid, treat = treat,
                                     X = cbind(z1, z2), 
                                     response = "normal", 
                                     shrinkage = "none"), 
                                     type = "random")
```

The JAGS model can be accessed as follows:

```{r MAcont2}
ipd$model.JAGS
```
We can fit the treatment effect model as follows:

```{r MAcont3}
#| results: hide
#| message: false
samples <- ipd.run(ipd, n.chains = 2, n.iter = 20,
                   pars.save = c("alpha", "beta", "delta", "sd", "gamma"))
trtbenefit <- round(treatment.effect(ipd, samples, newpatient = c(z1 = 1, z2 = 0.5)), 2)
```

Here are the estimated model parameters:

```{r MAcont4}
summary(samples)
```

#### Prection
We can now predict the individualized treatment effect for a new patient with covariate values `z1 = 1` and `z2 = 0.5`.

```{r MAcont5}
round(treatment.effect(ipd, samples, newpatient = c(z1 = 1, z2 = 0.5)), 2)
```

This means that the predicted outcome for patient with covariate values `z1 = 1` and `z2 = 0.5` will differ by `r trtbenefit["0.5"]` units when receiving the active treatment (`treat = 1`) as compared to the control treatment (`treat = 0`).

We can also predict treatment benefit for all patients in the sample, and look at the distribution of predicted benefit.

```{r benMA}
#| warning: false
#| label: fig-predben_continuous_outcome
#| fig-cap: "Distribution of predicted treatment benefit in each trial. Dashed lines represent the trial mean."
#| fig-width: 8
#| fig-height: 7
library(dplyr)
library(ggplot2)

ds <- ds %>% mutate(benefit = NA,
                    study = paste("Trial", studyid)) 

for (i in seq(nrow(ds))) {
  newpat <- as.matrix(ds[i, c("z1", "z2")])
  ds$benefit[i] <- treatment.effect(ipd, samples, newpatient = newpat)["0.5"]
}

summbenefit <- ds %>% group_by(study) %>% 
  summarize(mediabenefit = median(benefit), meanbenefit = mean(benefit))

ggplot(ds, aes(x = benefit)) + 
  geom_histogram(aes(y = after_stat(density)), alpha = 0.3) + 
  geom_density() +
  geom_vline(data = summbenefit, aes(xintercept = meanbenefit), 
             linewidth = 0.5, lty = 2) + 
  facet_wrap(~study) + 
  ylab("Density") +
  xlab("Predicted treatment benefit")  + theme_bw()
```

#### Penalization
Let us repeat the analysis, but this time while penalizing the treatment-covariate coefficients using a Bayesian LASSO prior. 

```{r MAcont6}
ipd <- with(ds, ipdma.model.onestage(y = y, study = studyid, 
                                     treat = treat,
                                     X = cbind(z1, z2), 
                                     response = "normal", 
                                     shrinkage = "laplace"), 
            type = "random")

samples <- ipd.run(ipd, n.chains = 2, n.iter = 20, 
                   pars.save = c("alpha", "beta", "delta", "sd", "gamma"))

round(treatment.effect(ipd, samples, newpatient = c(1,0.5)), 2)
```

### Example of a binary outcome


#### Setup

We now present the case of a binary outcome. We first generate a dataset as before, using the **bipd** package.

```{r MAbin1}
ds2 <- generate_ipdma_example(type = "binary")
head(ds2)
```

The simulated dataset contains information on the following variables:

- the trial indicator `studyid`
- the treatment indicator `treat`, which takes the values 0 for control and 1 for active treatment
- two prognostic variables `w1` and `w2`
- the binary outcome `y`

```{r}
#| label: tbl-summary-binary_outcome-data
#| tbl-cap: The simulated dataset with a binary outcome
#| echo: false
#| warning: false
dstbl <- ds2 %>% mutate(treat = as.factor(treat),
                        studyid = as.factor(studyid))
table1(~ w1+w2+studyid | treat, data = dstbl)
```

#### Model fitting
We use a Bayesian random effects model with binomial likelihood. This is similar to the model 16.7 of the book, but with a Binomial likelihood, i.e. 

$$ 
y_{ij}\sim \text{Binomial}(\pi_{ij}) \\
$$ 
$$ 
\text{logit}(\pi_{ij})=a_j+\delta_j t_{ij}+ \sum_{l=1}^{L}\beta_l x_{ij}+ \sum_{l=1}^{L}\gamma_l x_{ij} t_{ij}
$$
The remaining of the model is as in the book. 
We can  penalize the estimated parameters for effect modification ($\gamma$'s), using a Bayesian LASSO. We can do this using again the *bipd* package:
```{r MAbin2}
ipd2 <- with(ds2, ipdma.model.onestage(y = y, study = studyid, treat = treat,
                                       X = cbind(w1, w2), 
                                       response = "binomial", 
                                       shrinkage = "laplace"), 
             type = "random", hy.prior = list("dunif", 0, 1))

ipd2$model.JAGS
```

```{r}
#| message: false
#| echo: false
#| results: hide
samples <- ipd.run(ipd2, n.chains = 2, n.iter = 20, 
                   pars.save = c("alpha", "beta", "delta", "sd", "gamma"))
trtpred <- round(treatment.effect(ipd2, samples, newpatient = c(w1 = 1.6, w2 = 1.3)), 2)
```

```{r}
#| message: false
#| eval: false
samples <- ipd.run(ipd2, n.chains = 2, n.iter = 20, 
                   pars.save = c("alpha", "beta", "delta", "sd", "gamma"))
summary(samples)
```

```{r}
#| message: false
#| echo: false
summary(samples)
```

The predicted treatment benefit for a new patient with covariates `w1 = 1.6` and `w2 = 1.3` is given as:

```{r}
#| message: false
round(treatment.effect(ipd2, samples, newpatient = c(w1 = 1.6, w2 = 1.3)), 2)
```

In other words, the aforementioned patient `r trtpred["0.5"]` (95\% Credibility Interval:  `r trtpred["0.025"]` to  `r trtpred["0.975"]`)


## Estimating heterogeous treatment effects in network meta-analysis
###  Example of a continuous outcome
#### Setup
We use again the bipd package to simulate a dataset:
```{r NMA1}
ds3 <- generate_ipdnma_example(type = "continuous")
head(ds3)
```

Let us look into the data a bit in more detail:

```{r NMA2}
#| label: tbl-nma-summary-continuous_outcome-data
#| tbl-cap: The simulated dataset with a continuous outcome
#| echo: false
#| warning: false
dstbl <- ds3 %>% mutate(treat = as.factor(treat),
                        studyid = as.factor(studyid))
table1(~ z1+z2+studyid | treat, data=dstbl)
```

#### Model fitting
We will use the model shown in Equation 16.8 in the book. In addition, we will use Bayesian LASSO to penalize the treatment-covariate interactions.

```{r NMA3}
ipd3 <- with(ds3, ipdnma.model.onestage(y = y, study = studyid, treat = treat, 
                                        X = cbind(z1, z2), 
                                        response = "normal", 
                                        shrinkage = "laplace", 
                                        type = "random"))
ipd3$model.JAGS


samples <- ipd.run(ipd3, n.chains = 2, n.iter = 20, 
                   pars.save = c("alpha", "beta", "delta", "sd", "gamma"))
summary(samples)
```

As before, we can use the `treatment.effect()` function of *bipd* to estimate relative effects for new patients. 
```{r NMA4}
treatment.effect(ipd3, samples, newpatient= c(1,2))
```
This gives us the relative effects for all treatments versus the reference. To obtain relative effects between active treatments we need some more coding:

```{r NMA5}
samples.all=data.frame(rbind(samples[[1]], samples[[2]]))
newpatient= c(1,2)
newpatient <- (newpatient - ipd3$scale_mean)/ipd3$scale_sd

median(
  samples.all$delta.2.+samples.all$gamma.2.1.*
    newpatient[1]+samples.all$gamma.2.2.*newpatient[2]
-
  (samples.all$delta.3.+samples.all$gamma.3.1.*newpatient[1]+
     samples.all$gamma.3.2.*newpatient[2])
)

quantile(samples.all$delta.2.+samples.all$gamma.2.1.*
           newpatient[1]+samples.all$gamma.2.2.*newpatient[2]
         -(samples.all$delta.3.+samples.all$gamma.3.1.*newpatient[1]+
             samples.all$gamma.3.2.*newpatient[2])
         , probs = 0.025)

quantile(samples.all$delta.2.+samples.all$gamma.2.1.*
           newpatient[1]+samples.all$gamma.2.2.*newpatient[2]
         -(samples.all$delta.3.+samples.all$gamma.3.1.*newpatient[1]+
             samples.all$gamma.3.2.*newpatient[2])
         , probs = 0.975)
```

### Modeling patient-level relative effects using randomized and observational evidence for a network of treatments

We will now follow Chapter 16.3.5 from the book. In this analysis we will not use penalization, and we will assume fixed effects. For an example with penalization and random effects, see part 2 of this vignettte.

#### Setup
We generate a very simple dataset of three studies comparing three treatments. We will assume 2 RCTs and 1 non-randomized trial:

```{r NMA6}
ds4 <- generate_ipdnma_example(type = "continuous")
ds4 <- ds4 %>% filter(studyid %in% c(1,4,10)) %>%
  mutate(studyid = factor(studyid) %>%
           recode_factor(
             "1" = "1",
             "4" = "2",
             "10" = "3"),
         design = ifelse(studyid == "3", "nrs", "rct"))
```
The sample size is as follows:

```{r}
#| echo: false
t1=table(ds4$treat, ds4$studyid)
row.names(t1)=c("treat A:", "treat B:", "treat C:")
colnames(t1)=paste("s",1:3, sep="")
t1
```

#### Model fitting
We will use the design-adjusted model, equation 16.9 in the book. We will fit a two-stage fixed effects meta-analysis and we will use a variance inflation factor. The code below is used to specify the analysis of each individual study. Briefly, in each study we adjust the treatment effect for the prognostic factors `z1` and `z2`, as well as their interaction with `treat`.

```{r NMA7}
library(rjags)
first.stage <- "
model{

for (i in 1:N){
	y[i] ~ dnorm(mu[i], tau)  
	mu[i] <- a + inprod(b[], X[i,]) + inprod(c[,treat[i]], X[i,]) + d[treat[i]] 
}
sigma ~ dunif(0, 5)
tau <- pow(sigma, -2)

a ~ dnorm(0, 0.001)

for(k in 1:Ncovariate){
	b[k] ~ dnorm(0,0.001)
}

for(k in 1:Ncovariate){
	c[k,1] <- 0
}

tauGamma <- pow(sdGamma,-1)
sdGamma ~ dunif(0, 5)

for(k in 1:Ncovariate){
	for(t in 2:Ntreat){
		c[k,t] ~ ddexp(0, tauGamma)
	}
}

d[1] <- 0
for(t in 2:Ntreat){
	d[t] ~ dnorm(0, 0.001)
}
}"
```

Subsequently, we estimate the relative treatment effects in the first (randomized) study comparing treatments A and B:

```{r}
#| echo: false
# Keep track of study results
results <- data.frame(study	 = character(),
                      BA_est = numeric(),
                      BA_se = numeric(),
                      CA_est = numeric(),
                      CA_se = numeric())
```

```{r}
model1.spec <- textConnection(first.stage) 
data1 <- with(ds4 %>% filter(studyid == 1), 
              list(y = y,
                   N = length(y), 
                   X = cbind(z1,z2),  
                   treat = treat,
                   Ncovariate = 2, 
                   Ntreat = 2))
jags.m <- jags.model(model1.spec, data = data1, n.chains = 2, n.adapt = 500,
                     quiet =  TRUE)
params <- c("d", "c") 
samps4.1 <- coda.samples(jags.m, params, n.iter = 50)
samps.all.s1 <- data.frame(as.matrix(samps4.1))

samps.all.s1 <- samps.all.s1[, c("c.1.2.", "c.2.2.", "d.2.")]
delta.1 <- colMeans(samps.all.s1)
cov.1 <- var(samps.all.s1)
```

```{r}
#| echo: false
results <- results %>% add_row(data.frame(study = "study 1",
                                          BA_est = delta.1["d.2."],
                                          BA_se = sqrt(diag(cov.1)["d.2."])))  
```

We repeat the analysis for the second (randomized) study comparing treatments A and C:
```{r}
model1.spec <- textConnection(first.stage) 
data2 <- with(ds4 %>% filter(studyid == 2), 
              list(y = y,
                   N = length(y), 
                   X = cbind(z1,z2),  
                   treat = ifelse(treat == 3, 2, treat),
                   Ncovariate = 2, 
                   Ntreat = 2))
jags.m <- jags.model(model1.spec, data = data2, n.chains = 2, n.adapt = 100,
                     quiet =  TRUE)
params <- c("d", "c") 
samps4.2 <- coda.samples(jags.m, params, n.iter = 50)
samps.all.s2 <- data.frame(as.matrix(samps4.2))
samps.all.s2 <- samps.all.s2[, c("c.1.2.", "c.2.2.", "d.2.")]
delta.2 <- colMeans(samps.all.s2)
cov.2 <- var(samps.all.s2)
```

```{r}
#| echo: false
results <- results %>% add_row(data.frame(study = "study 2",
                                          CA_est = delta.2["d.2."],
                                          CA_se = sqrt(diag(cov.2)["d.2."])))   
```

Finally, we analyze the third (non-randomized) study comparing treatments A, B, and C:
```{r}
model1.spec <- textConnection(first.stage) 
data3 <- with(ds4 %>% filter(studyid == 3), 
              list(y = y,
                   N = length(y), 
                   X = cbind(z1,z2),  
                   treat = treat,
                   Ncovariate = 2, 
                   Ntreat = 3))
jags.m <- jags.model(model1.spec, data = data3, n.chains = 2, n.adapt = 100,
                     quiet = TRUE)
params <- c("d", "c") 
samps4.3 <- coda.samples(jags.m, params, n.iter = 50)
samps.all.s3 <- data.frame(as.matrix(samps4.3))

samps.all.s3 <- samps.all.s3[, c("c.1.2.", "c.2.2.", "d.2.", "c.1.3.", 
                                 "c.2.3.", "d.3.")]
delta.3 <- colMeans(samps.all.s3)
cov.3 <- var(samps.all.s3)
```

```{r}
#| echo: false
# Save the results from the third analysis
results <- results %>% add_row(data.frame(study = "study 3",
                                          BA_est = delta.3["d.2."],
                                          BA_se = sqrt(diag(cov.3)["d.2."]),
                                          CA_est = delta.3["d.3."],
                                          CA_se = sqrt(diag(cov.3)["d.3."])))  
```

The corresponding treatment effect estimates are depicted below:

```{r}
#| echo: false
#| label: tbl-results_nma_stage1
#| tbl-cap: "Treatment effect estimates."

results <- results %>% mutate(BA = ifelse(is.na(BA_est), "", paste(sprintf("%.3f", BA_est), "(SE = ", sprintf("%.3f", BA_se), ")")),
                              CA = ifelse(is.na(CA_est), "", paste(sprintf("%.3f", CA_est), "(SE = ", sprintf("%.3f", CA_se), ")"))) %>% 
  select(study, BA, CA)

kable(results,
      booktabs = TRUE, digits = 2,
      col.names = c("study", "B versus A", "C versus A")) %>%
  kable_styling(latex_options = "hold_position")
```

We can now fit the second stage of the network meta-analysis. The corresponding JAGS model is specified below:
```{r}
second.stage <-
"model{
  
  #likelihood
  y1 ~ dmnorm(Mu1, Omega1)
  y2 ~ dmnorm(Mu2, Omega2)
  y3 ~ dmnorm(Mu3, Omega3*W)

  
  Omega1 <- inverse(cov.1)
  Omega2 <- inverse(cov.2)
  Omega3 <- inverse(cov.3)

  Mu1 <- c(gamma[,1], delta[2])
  Mu2 <- c(gamma[,2], delta[3])  
  Mu3 <- c(gamma[,1], delta[2],gamma[,2], delta[3])
  
  #parameters
  for(i in 1:2){
    gamma[i,1] ~ dnorm(0, 0.001)
    gamma[i,2] ~ dnorm(0, 0.001)
  }
  
  delta[1] <- 0
  delta[2] ~ dnorm(0, 0.001)
  delta[3] ~ dnorm(0, 0.001)
  
}
"
```

We can fit as follows:
```{r}
model1.spec <- textConnection(second.stage) 
data3 <- list(y1 = delta.1, y2 = delta.2, y3 = delta.3, 
              cov.1 = cov.1, cov.2 = cov.2, cov.3 = cov.3, W = 0.5)

jags.m <- jags.model(model1.spec, data = data3, n.chains = 2, n.adapt = 50,
                     quiet = TRUE)
params <- c("delta", "gamma") 
samps4.3 <- coda.samples(jags.m, params, n.iter = 50)
```

```{r}
summary(samps4.3)

# calculate  treatment effects
samples.all = data.frame(rbind(samps4.3[[1]], samps4.3[[2]]))
newpatient = c(1,2)

median(
  samples.all$delta.2. + samples.all$gamma.1.1.*newpatient[1] +
    samples.all$gamma.2.1.*newpatient[2]
)

quantile(samples.all$delta.2.+samples.all$gamma.1.1.*newpatient[1]+
           samples.all$gamma.2.1.*newpatient[2]
         , probs = 0.025)

quantile(samples.all$delta.2.+samples.all$gamma.1.1.*newpatient[1]+
           samples.all$gamma.2.1.*newpatient[2]
         , probs = 0.975)

```


## Version info {.unnumbered}
This chapter was rendered using the following version of R and its packages:

```{r}
#| echo: false
#| message: false
#| warning: false
sessionInfo()
```

## References {.unnumbered}