```
library(ggplot2)
library(dplyr)
library(tidyr)
library(brms)
library(survival)
library(knitr)
library(posterior)
library(bayesplot)
library(simsurv)
library(RBesT)
library(survminer)
library(gt)
library(here)
# instruct brms to use cmdstanr as backend and cache all Stan binaries
options(brms.backend="cmdstanr", cmdstanr_write_stan_file_dir=here("_brms-cache"))
# create cache directory if not yet available
dir.create(here("_brms-cache"), FALSE)
# allow for wider printouts without line-breaking
options(width=120)
```

# 14 Time-to-event data

In this case study the design of a time-to-event trial is presented. The trial evaluates a combination treatment of an active drug given in addition to standard of care (SoC) therapy. Two SoC therapies are in clinical use, which varies geographically. The SoC therapies alone are used as control and compared to the combination of the active drug given in addition to the respective SoC. This results in four arms of the trial and special emphasis is put on the parametrization of the problem. That is, the treatment effect (and each standard of care) is anticipated to be similar to one another such that the mean treatment effect is of particular interest. In addition, a historical data set is incorporated in the analysis of the trial.

This case study demonstrates:

- modeling censored time-to-event data in
`brms`

(use of`cens`

formula modifier) with borrowing from historical control data, - the use of custom contrasts allowing for custom definition of the model parametrization, which in turn allows to specify priors on relevant bits of the problem (a priori we expect the treatment effect not to differ substantially between standard of cares),
- an additional alternative approaches of using custom coded contrasts to include with greater flexibility historical data.

To run the R code of this section please ensure to load these libraries and options first:

## 14.1 Background

The outcome of interest in clinical trials is often the occurence of an event. Such events might be negative for patients such as a worsening of a disease, hospitalization or death. They might also be a positive such as clearing the blood of malaria parasites, or being released from hospital.

We use time-to-event (or survival) analysis when it is important to patients whether the event of interest occurs earlier or later. Additionally, time-to-event methods let us deal with the fact that not all patients will have the event of interest during a trial, which leads to censored observations.

In this case study we focus on an Oncology late phase trial evaluating a test treatment in combination with standard of care (SoC) using progression free survival as endpoint in the indication of gastric cancer. Due to geographic differences in the SoC, two different SoC treatments are considered as control treatments. One active treatment is tested in combination with each SoC and compared to the control treatment using only the SoC, which leads to four trial arms in total. From clinical experience with the two SoCs on expects a similar PFS for both therapies. In addition, it is expected that the treatment effect of the active drug is similar when combined with each SoC respectivley. It is therefore desirable to express this prior knowledge when setting up the model.

By way of aligning the parametrization of the model we can encode the prior knowledge on expected similarities (no difference between the two SoCs and consistent treatment effect). For example, defining the average treatment effect and the difference between the two treatment effects as parameters, one can place a prior with limited width on the difference parameter. This expresses the a-priori expectation that differences in the treatment effect per SoC are not too large. This results in efficency gains for the estimation of the average treatment effect while only partially pooling the data.

As furthermore historical data is available on each SoC treatment, its use may allow for un-equal randomization between active and control arms. However, the available literature data only reports outcomes for a treatment arm which lumps together data from the two SoC treatments. Since about 50% of these patients were treated with either SoC, the reported data can be considered to report the average effect of both SoCs. How this historical data can be included in the analysis is outlined below.

## 14.2 Modeling time-to-event data with `brms`

Modeling time to event data of clinical trials is oftentimes run with the semi-parametric Cox proportional hazards model. This model avoids the need to specify a hazard function \(h(t)\), which is the rate of events provided no event as happened yet and it crucially defines the probability density of the survival times. Using the assumption of proportional hazards allows to quantify effects of explanatory variables on the hazard and thereby one can quantify the difference between two treatment groups in a trial, for example. The hazard for a patient \(i\) with covariates \(x_i\) is often modelled as \(h_i(t) =
\exp\{x_i \, \beta\} \, h_0(t)\). The baseline hazard function \(h_0(t)\) is not required in the estimation of the effects \(\beta\) and the Cox model parameters are estimated using a *partial* likelihood approach. However, as the Bayesian framework is based on the (full) likelihood principle, the semi-parametric Cox proportional hazards model is not applicale. Nonetheless, `brms`

does offer a variant of the Cox proportional hazards model, which is designed to result in numerically very similar results if used with non-informative priors when compared to the respective Cox proportional hazards model. Instead of marginalizing out the baseline hazard function \(h_0(t)\), the baseline hazard function is modelled using an almost non-parametric approach based on a parametric spline approximation to \(h_0(t)\). The estimated regression coefficients \(\beta\) for the Cox `brms`

family correspond to logarithmic hazard ratios.

In contrast to the semi-parametric Cox model, parametric modeling of time to event data assumes a functional form for the hazard function \(h(t)\) or equivalently for the probability density function \(f(t)\) of the event times at \(t=T\). Both functions are related to one another by the basic relationship \(h(t) = \frac{f(t)}{S(t)}\), where \(S(t) = P(T>t) =
\int_t^\infty f(u)\, du\) is the survivor function (probability that an event occurs past time \(t\)). In `brms`

the convention is to define via the `family`

argument to `brm`

the probability density of the event times. A description of the parametrization for these densities are found in the vignette “Parametrization of Response Distributions in `brms`

”.

In this case study the Weibull distribution is used as literature data suggested its appropiateness and a parametric modeling approach may lead to greater statistical efficiency. The Weibull probability density is parametrized in `brms`

in terms of shape \(\alpha\) and scale parameter \(s\). Instead of directly modeling the scale \(s\) (as done in many other time to event programs), `brms`

models the mean of the Weibull distribution \(\mu\). Therefore, the scale \(s\) is set equal to \(s=\frac{\mu}{\Gamma(1 + \frac{1}{\alpha})}\) in the Weibull probability density

\[ f_{\mbox{Weibull}}(t) = \frac{\alpha}{s} \left( \frac{t}{s} \right)^{\alpha-1} \, \exp\left(-\left(\frac{t}{s}\right)^\alpha\right).\]

In this form, the model fulfills the property of an *accelerated failure time* (AFT) model whenever explanatory variables are introduced. This follows from considering the survivor function

\[S_{\mbox{Weibull}}(t) = \exp\left(-\left(\frac{t}{s}\right)^\alpha\right)\]

and modeling the mean \(\mu_i\) linearly on the \(\log\) scale as a function of covariates \(x_i\) for subject \(i\), \(\log(\mu_i) = \beta_0 + x_i' \, \beta\). Defining as the reference covariate level \(x_i = 0\) motivates the definition of a reference survivor function \(S_{\mbox{Weibull},0}(t)\) for which the linear predictor \(\log(\mu_0)\) is equal to the intercept \(\beta_0\). As a consequence, the survivor function of any patient \(i\) is related to the reference survivor function by

\[ S_{i}(t) = S_{\mbox{Weibull},0}\left(\frac{t}{\exp(x_i' \, \beta)} \right),\]

which is the defining property of AFT models. The regression coefficients \(\beta\) are then interpretable as relative speedup/slowdown of the process progression. That is, an increased time scale (slowdown) leads to a delay of an event. Given that modeling the scale of the Weibull distribution is a common convention (instead of the mean \(\mu\) as in `brms`

), it is useful to recall that this merely means that we need to offset the intercept as estimated by `brms`

with \(-\log\Gamma\left(1 + \frac{1}{\alpha} \right)\) in order to obtain the scale on the log scale, \(\log(s)\) (which is the what the `survival`

R package would report with its `survreg`

routine). The `survival`

R package does futhermore use a log-linear representation of the statistical problem leading to an estimation of the inverse shape parameter \(\sigma=\frac{1}{\alpha}\).

The Weibull AFT model as estimated by `brms`

can be converted into the respective proportional hazard model by transforming the scale \(s\) with the relation \(s = \lambda^{-\frac{1}{\alpha}}\) to an alternative scale parameter \(\lambda\). Doing so casts the hazard function into

\[h_{\mbox{Weibull}}(t) = \lambda \, \alpha \, t^{\alpha-1}.\]

When we now model \(\lambda\) as a function of covariates \(x_i\) linearly on the \(\log\) scale one arrives at

\[ h_i(t) = \exp(x_i' \, \xi) \, h_{\mbox{Weibull},0}(t),\]

where \(h_{\mbox{Weibull},0}(t)\) is defined by the reference covariate level \(x_i=0\) such that \(\lambda = \exp(\xi_0)\) for \(h_{\mbox{Weibull},0}(t)\). To now convert the AFT regression coefficients \(\beta\) as estimated by `brms`

to their respective proportional hazard coefficients \(\xi\) one may just apply the transformation resulting in the relationship

\[\xi = -\alpha \, \beta,\]

which converts from logarithmic speedup/slowdown \(\beta\) (AFT model) into logarithmic hazard ratios \(\xi\) (proportional hazard model).

## 14.3 Data

We demonstrate here the trial analysis using a fake data simulation of the trial design. As key modeling choices for the parametric modeling approach are motivated from the historical data we start with the presentation of the historical data set and then describe the details of the trial simulation.

### 14.3.1 Historical data

The CheckMate 649 trial (Janjigian et al. 2021) included 789 gastric cancer patients which were treated with either `chemoA`

or `chemoB`

and their progression free survival (PFS) was reported. Both chemotherapies were used in a roughly 1:1 ratio such that we will consider the reported data as “average” between both chemo therapies (despite lack of randomization occured wrt to chemotherapy assignment).

The data of trial has here been reconstructed from the published survival curves:

To now evaluate if the Weibull probability density function is an appropiate choice for this data we consider if the non-parametric estimate of the survivor function is compatible with properties of Weibull distributed event times. Specifically, if we transform the survivor function \(S_{\mbox{Weibull}}(t)\) of the Weibull probability density function with the complementary log-log transformation, we obtain that

\[ \mbox{cloglog}(S_{\mbox{Weibull}}(t)) = \log( - \log(S_{\mbox{Weibull}}(t))) = -\alpha \, \log(s) + \alpha \, \log(t)\]

holds. Therefore, we can visualize an estimate of the survivor function on a transformed scale as a function of \(\log(t)\) and we should observe a straight line with slope \(\alpha\) and intercept \(-\alpha \, \log(s)\):

```
<- survfit(Surv(time, status) ~ 1, data=hdata2)
km ggsurvplot(km, data=hdata2, fun = "cloglog")
```

While the line is not perfectly straight over the entire follow-up, the assumption of a straight line within the time period including the bulk of all events falls in the range of a good approximation with a line (as can be seen from the numbers at risk of the plot above).

As an additional evaluation of the Weibull distributional form we model the CheckMate 649 data via `brms`

and perform a posterior predictive check. The priors are specified as illustrated in the later section on the used priors. The result looks reasonable.

```
<- bf(time | cens(1-status) ~ 1, family=weibull())
model_weibull_hist
<- prior(normal(log(6), log(4)/1.64), class=Intercept) +
prior_weibull_hist prior(gamma(3, 2.7), class=shape)
<- brm(model_weibull_hist,
fit_hist_checkmate data=hdata2,
prior=prior_weibull_hist,
seed=234235,
refresh=0)
```

```
Running MCMC with 4 sequential chains...
Chain 1 finished in 1.3 seconds.
Chain 2 finished in 1.3 seconds.
Chain 3 finished in 1.3 seconds.
Chain 4 finished in 1.2 seconds.
All 4 chains finished successfully.
Mean chain execution time: 1.3 seconds.
Total execution time: 5.5 seconds.
```

` fit_hist_checkmate`

```
Family: weibull
Links: mu = log; shape = identity
Formula: time | cens(1 - status) ~ 1
Data: hdata2 (Number of observations: 789)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.44 0.04 2.37 2.51 1.00 2794 2455
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 1.22 0.04 1.14 1.30 1.00 2865 2574
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
```

With the posterior we can now perform a posterior predictive check. We do so by drawing 100 samples of the model parameters from the posterior and draw for each sample a realization of the event times for all patients. Hence, for each sample one obtains a fake data set simulated according to the Weibull distribution. Finally, each fake data set is summarized by a Kaplan-Meier estimate and drawn in light color while the actually observed data is summarized in the same way and drawn using a thick line. If the fitted model is appropiate, then the observed Kaplan Meier estimate should look just like a random realization among the many fake data generated ones. However, we observe that over the entire follow-up time the observed data is leaving for late follow-up times the bulk of the fake data generated Kaplan Meier curves. This may suggest a deviation from the Weibull distribution or the possibility that time-varying covariates are helpful. Nonetheless, if one focusses on the follow-up time frame, which includes the majoritiy of events (up to 15 month), then the posterior predictive check does look reasonable. This is why we use the Weibull distribution without a time-varying covariate for simplicity.

```
<- pp_check(fit_hist_checkmate, status_y=hdata2$status,
p_full_fup type="km_overlay", ndraws=100) +
scale_y_continuous(breaks=seq(0,1,by=0.1)) +
xlab("Time [month]") + coord_cartesian(xlim=c(0, 35))
<- p_full_fup + coord_cartesian(xlim=c(0, 15))
p_sub_fup ::ggarrange(p_full_fup, p_sub_fup, common.legend = TRUE) ggpubr
```

Note that the above predictive check does not simulate the censoring process. Thus, the censoring process is ignored and non-informative censoring is thereby assumed.

### 14.3.2 Trial simulation

Here we consider a fake data simulation of the trial of interest evaluating a novel treatment of gastric cancer for which an established treatment is available already. The trial is a randomized trial to study the efficacy and safety of adding a drug `CompoundX`

in combination with a monoclonal antibody `mAb`

+ standard of care (SoC). Two chemotherapy backbones, `ChemoA`

and `ChemoB`

were used as SoC. There are in total four treatment arms including two control groups, which are `mAb`

+`ChemoA`

(`controlChemoA`

) and `mAb`

+`ChemoB`

(`controlChemoB`

), and two corresponding active arms with the test treatment administrated in addition, that is `CompoundX`

+`mAb`

+`ChemoA`

(`activeChemoA`

), `CompoundX`

+`mAb`

+`ChemoB`

(`activeChemoB`

).

Synthetic data were generated for the randomised study. We let the sample size of both active and control group be 100, respectively. The sample size has been chosen to ensure adequate trial power of 80%. Efficacy was assessed using progression-free survival (PFS) measured in units of months. For patients not experiencing a progression event at the time of the analysis, their event time will be right censored at the time of their last valid tumor assessment. For patients experiencing an event, their event time will be interval censored, with the upper limit of the interval being the study day on which the progression event is observed, and the lower limit being the day of the last preceding valid tumor assessment at which the patient’s progression-free status was confirmed. While interval censoring is a supported feature of `brms`

, this case study will for simplicity use the actual event times and leave it to the reader to explore interval censoring with `brms`

.

There are two historical datasets available. The first one with 400 patients on a control arm, `mAb`

+`ChemoA`

(`controlChemoA`

), is simulated together with the randomised study. And the PFS data of 789 patients receiving first-line programmed cell death (PD-1) inhibitor+chemotherapy are available from a randomised, open-label, phase 3 trial, CheckMate 649 (Janjigian et al. 2021), and it will be introduced later. These data are highly relevant and we would like to integrate the historical information into our analysis.

Simulated data will be generated using `simsurv`

function via simulating survival times from the Weibull model using the proportional hazard form of the model. We assumed a 20% hazard rate (HR) reduction due to active treatment administration (`CompoundX`

) to the active group in comparison to the control group and a 5% difference in HR for the two chemotherapy SoCs (`ChemoB`

is 5% worse than `ChemoA`

on the hazard scale). We also assumed that both of the historical data have a 2% worse outcome on the hazards scale and we represent the two studies with `hist1`

and `hist2`

respectively (`hist1`

for simulated historidal data and `hist2`

for CheckMate 649).

## Show the code

```
set.seed(46767)
# n per group
<- 100
n_grp <- 400
n_hist # use month as time-unit
<- 1 / 6
rate_1 <- 1 / 10
rate_cens <- c(trt=-0.2, ## roughly 20% HR reduction
beta soc_alt=0.05, ## alternative chemotherapy is 5% worse
hist1=0.02) ## simulated historical data has a 2% worse outcome
# covariates of simulated trial data
<- data.frame(id = seq(1, 2*n_grp), trt = c(0L, 1L),
covs soc_alt=rbinom(2*n_grp, 1L, 0.3), hist1=0L, hist2=0L)
# covariates of historical data
<- data.frame(id = seq(2*n_grp+1, 2*n_grp+1 +n_hist - 1),
hcovs trt = c(0), soc_alt=0, hist1=1L, hist2=0L)
<- function(lambda, gamma, lambda_cens, x, betas) {
simulate_trial ## simulate censoring times, note that we do not simulate end of
## trial with maxt for now. Also note the different parametrization of
## simsurv which models log(shape) with a linear predictor.
<- simsurv(lambdas = lambda_cens, gammas = 1, x=x)
cens <- simsurv(lambdas = lambda, gammas = gamma, x=x, betas=betas)
events names(cens) <- paste0(names(cens), "_cens")
bind_cols(events, select(cens, -id_cens), select(x, -id)) %>%
rename(censtime=eventtime_cens) %>%
mutate(event=1L*(eventtime <= censtime),
y=if_else(event==1, eventtime, censtime),
status=NULL, status_cens=NULL) %>%
relocate(id, y, event) %>%
mutate(soc=factor(soc_alt, c(0, 1), c("ChemoA", "ChemoB")),
trt_ind=trt,
trt=factor(trt_ind, c(0,1), c("control", "active")),
arm=factor(paste0(c("control", "active")[trt_ind + 1], soc),
levels=c("activeChemoA", "controlChemoA",
"activeChemoB", "controlChemoB")))
}
# using shape=1 for simulation (corresponding to exponential survival times)
<- simulate_trial(rate_1, 1, rate_cens, covs, beta)
sim <- simulate_trial(rate_1, 1, rate_cens, hcovs, beta) %>%
hdata1 mutate(id=id+max(sim$id))
```

We can check how the simulated data looks like and its sample size.

`gt_preview(sim) %>% fmt_number(where(is.double), decimals=1)`

id | y | event | eventtime | censtime | trt | soc_alt | hist1 | hist2 | soc | trt_ind | arm | |
---|---|---|---|---|---|---|---|---|---|---|---|---|

1 | 1 | 7.7 | 0 | 12.3 | 7.7 | control | 0 | 0 | 0 | ChemoA | 0 | controlChemoA |

2 | 2 | 0.1 | 0 | 20.4 | 0.1 | active | 0 | 0 | 0 | ChemoA | 1 | activeChemoA |

3 | 3 | 4.7 | 0 | 10.3 | 4.7 | control | 0 | 0 | 0 | ChemoA | 0 | controlChemoA |

4 | 4 | 2.7 | 0 | 19.0 | 2.7 | active | 0 | 0 | 0 | ChemoA | 1 | activeChemoA |

5 | 5 | 3.6 | 1 | 3.6 | 7.5 | control | 0 | 0 | 0 | ChemoA | 0 | controlChemoA |

6..199 | ||||||||||||

200 | 200 | 11.8 | 0 | 17.1 | 11.8 | active | 0 | 0 | 0 | ChemoA | 1 | activeChemoA |

`table(sim$arm)`

```
activeChemoA controlChemoA activeChemoB controlChemoB
72 73 28 27
```

`table(hdata1$arm)`

```
activeChemoA controlChemoA activeChemoB controlChemoB
0 400 0 0
```

## 14.4 Models

### 14.4.1 Model specification

Survival data will be fitted by a Bayesian generalized linear model.

The survival function follows a Weibull distribution with shape parameter \(\alpha\) and scale parameter \(s\),

\[T|\alpha,s \sim f_{\mbox{Weibull}}(T).\]

A common shape parameter \(\alpha\) is assumed for all treatment arms, with the scale parameter allowed to vary by treatment arm. In place of directly modeling the scale parameter \(s\), the mean of the Weibull distribution is modelled on the logarithmic scale such that the scale is set to \(s = \frac{\mu}{\Gamma(1 + \frac{1}{\alpha})}\) and \(\log(\mu)\) is modelled as a linear function of the covariates. We use \(\lambda_{\mbox{activeChemoA}}\) to represent the logarithm of the mean survival time for patients receiving `CompoundX`

+`mAb`

+`ChemoA`

(\(\lambda_{\mbox{activeChemoA}} = \log(\mu_{\mbox{activeChemoA}})\)), and \(\lambda_{\mbox{controlChemoA}}\) for patients receiving `mAb`

+`ChemoA`

. We define \(\lambda_{\mbox{activeChemoB}}\) and \(\lambda_{\mbox{controlChemoB}}\) similarly for patients receiving the `ChemoB`

chemotherapy backbone.

With 4 different patient groups we may define 4 parameters in total. The default parametrization in `R`

for categorical variables is the treatment contrast representation. This representation delcares a reference category to form the overall intercept and defines differences to the overall intercept for each categorical level. As this does automatically define the parametrization of the model, we here define the contrasts in a customized manner. Doing so ensures that we control directly the exact parametrization of the model. This is important in this problem as we can then setup priors in a tailored manner. More details in defining custom contrasts are explained by Ben Bolker here.

In addition to the overall mean (global intercept)

\[\gamma = \frac{1}{4}( \lambda_{\mbox{activeChemoA}} + \lambda_{\mbox{controlChemoA}} + \lambda_{\mbox{activeChemoB}} + \lambda_{\mbox{controlChemoB}} ), \]

we define for the four groups the following contrasts of interest:

- Average difference between the active and control arms in each chemotherapy backbone:

\[\delta_{\mbox{EffectAvg}} = \frac{1}{2} { (\lambda_{\mbox{activeChemoA}} - \lambda_{\mbox{controlChemoA}}) + ( \lambda_{\mbox{activeChemoB}} - \lambda_{\mbox{controlChemoB}} ) } \]

- Half of the difference in treatment effect between the two chemotherapy backbones:

\[\delta_{\mbox{effect}} = \frac{1}{2} { (\lambda_{\mbox{activeChemoA}} - \lambda_{\mbox{controlChemoA}}) - ( \lambda_{\mbox{activeChemoB}} - \lambda_{\mbox{controlChemoB}}) }\]

- Difference between the two control arms:

\[\delta_{\mbox{control}} = - \lambda_{\mbox{controlChemoA}} + \lambda_{\mbox{controlChemoB}} \]

To formalize these contrast definitions, we summarize the arms to parameter mapping as matrix:

Parameter | $\lambda_{\mbox{activeChemoA}}$ | $\lambda_{\mbox{controlChemoA}}$ | $\lambda_{\mbox{activeChemoB}}$ | $\lambda_{\mbox{controlChemoB}}$ |
---|---|---|---|---|

\(\gamma\) |
1/4 |
1/4 |
1/4 |
1/4 |

\(\delta_{\mbox{EffectAvg}}\) |
1/2 |
-1/2 |
1/2 |
-1/2 |

\(\delta_{\mbox{effect}}\)$ |
1/2 |
-1/2 |
-1/2 |
1/2 |

\(\delta_{\mbox{control}}\) |
0 |
-1 |
0 |
1 |

In this form it is straightforward to derive the definition of a parameter and how it relates to the strata of the data based on the parameter definitions. However, in order to define the design matrix, which maps parameters to data strata, we will need the inverse of the above matrix. By convention `R`

requires for the specification of custom contrats, the matrix mapping parameters to the data means such that the above matrix is referred to as \(C^{-1}\) in the following.

The linear model for \(\log(\mu)\) is parameterised in terms of the defined contrasts with additional parameters allowing for heterogeneity between data from the two separate historical data strata. For the 4 different treatment arms the linear predictor is then given by:

\[\begin{align} \begin{bmatrix} \lambda_{\mbox{activeChemoA}} \\ \lambda_{\mbox{controlChemoA}} \\ \lambda_{\mbox{activeChemoB}} \\ \lambda_{\mbox{controlChemoB}} \end{bmatrix} &= C \, \begin{bmatrix} \gamma \\ \delta_{\mbox{EffectAvg}} \\ \delta_{\mbox{effect}} \\ \delta_{\mbox{control}} \end{bmatrix} + I_{\mbox{hist1}} \, \beta_{\mbox{hist1}} + I_{\mbox{hist2}} \, \beta_{\mbox{hist2}} \end{align}\]

The indicators \(I_{\mbox{hist1}}\) and \(I_{\mbox{hist2}}\) represent that a patient was in the simulated historical stratum, or in the historical CheckMate 649 stratum, respectively. The corresponding parameters \(\beta_{\mbox{hist1}}\) and \(\beta_{\mbox{hist2}}\) capture differences in outcome between these two historical strata, and patients in our randomised comparison. In order to borrow strength from the historical data, the prior on these two regression coefficients are set reasonably narrow which in effect emulates a random effects meta-analysis.

As the historical data set reported in CheckMate 649 reported the outcome for the control treatment with a mixed population of 50% being on either SoC, we consider this data to report an average SoC effect. Basic algebra shows that \[\frac{1}{2} \, (\lambda_{\mbox{activeChemoA}} + \lambda_{\mbox{activeChemoB}}) = \gamma - \frac{1}{2} \, \delta_{\mbox{effect}}. \] Thus, we may simply add an additional row to the matrix \(C\) such that the full matrix \(C\) is:

Treatment arm | $\gamma$ | $\delta_{\mbox{EffectAvg}}$ | $\delta_{\mbox{effect}}$ | $\delta_{\mbox{control}}$ |
---|---|---|---|---|

activeChemoA |
1 |
1/2 |
1 |
-1/2 |

controlChemoA |
1 |
-1/2 |
0 |
-1/2 |

activeChemoB |
1 |
1/2 |
-1 |
1/2 |

controlChemoB |
1 |
-1/2 |
0 |
1/2 |

controlAverage |
1 |
-1/2 |
0 |
0 |

### 14.4.2 Priors

Priors for each random variable are described in the following table.

Parameter | Prior |
---|---|

\(\alpha\) |
\(\Gamma(3, 2.7)\) |

\(\gamma\) |
\(\mbox{Normal}( \log(\frac{8}{\log(2)}), \frac{\log(4)}{1.64})\) |

\(\delta_{\mbox{effectAvg}}\) |
\(\mbox{Normal}(0, \frac{\log(2)}{1.64})\) |

\(\delta_{\mbox{effect}}\) |
\(\mbox{Normal}(0, \frac{\log(1.25)}{1.64})\) |

\(\delta_{\mbox{control}}\) |
\(\mbox{Normal}(0, \frac{\log(1.25)}{1.64})\) |

\(\beta_{\mbox{hist1}}\) |
\(\mbox{Normal}(0, \frac{\log(1.8)}{1.64})\) |

\(\beta_{\mbox{hist2}}\) |
\(\mbox{Normal}(0, \frac{\log(1.8)}{1.64})\) |

These priors parametrize the central 90% credible interval for the respective normal priors used. As the parameters are additive on the \(\log\) scale, the standard deviation of \(\log(2)/\Phi^{-1}(95\%) \approx \log(2)/1.64\) defines a 90% credible interval of \([-\log(2), \log(2)]\) and implies on the back-transformed scale a respective interval of \([1/2, 2]\), i.e. a doubling or a halving of the progression speed.

The common shape parameter, \(\alpha\), is set to have median value of unity and it’s 90% credible interval is \([0.3, 2.33]\), which admits substantial deviation from unity, which is the case of exponential survival time distribution.

The prior for \(\gamma\), the intercept, is set with a mean approximately corresponding to expected median PFS on anti-PD-1 with SoC, as estimated from CheckMate 649 (approx. 8 months). The width of the central 90% credible interval for the intercept prior \(\gamma\) is set to admit for a 4 fold increase or decrease. A priori it is expected that the outcome for `controlChemoA`

and `controlChemoB`

will be similar (Janjijian et al, mPFS by chemotherapy not provided, but mOS was 14 months for both anti-PD1+Chemotherapy backbones (considered as similar to `mAb`

+`ChemoA`

and`mAb`

+`ChemoB`

in our simulated data), hence a prior for \(\delta_{\mbox{control}}\) is centered on 0 with a standard deviation corresponding to a 25% increase or decrease is chosen for the 90% credible interval. Similarly, the treatment benefit from adding `CompoundX`

is expected to be similar for both arms, so the prior for \(\delta_{\mbox{effect}}\) is again centered on 0 with the prior admitting a 25% increase or decrease a priori in terms of the 90% central credible interval. Priors for \(\beta_{\mbox{hist1}}\) and \(\beta_{\mbox{hist2}}\) are centered at a mean of 0, corresponding to no difference in outcome for the two historical studies as compared to this trial’s randomised part, but with the same variance allowing substantial differences between the randomised study and historical data. The use of Student-t priors for these parameters also robustifies against possible data discordance between the randomised comparison and historical data.

These prior definitions provide an uninformative prior, that allows for a broad range of shapes for the survival curve.

## 14.5 Implementation/Results

### 14.5.1 Construct customized coded contrasts

First we define customized contrasts which define the parametrization of the model as defined in the *Model sepcification* section.

The mapping of data strata to parameters is:

```
<- matrix(c(1/4, 1/4, 1/4, 1/4,
cc_inv 1/2, -1/2, 1/2, -1/2,
1/2, -1/2, -1/2, 1/2,
0, -1, 0, 1),
nrow=4, ncol=4, byrow=TRUE,
dimnames=list(contrast=c("intercept", "deltaEffectAvg",
"deltaEffect", "deltaControl"),
arm=c("activeChemoA", "controlChemoA",
"activeChemoB", "controlChemoB")))
::fractions(cc_inv) MASS
```

```
arm
contrast activeChemoA controlChemoA activeChemoB controlChemoB
intercept 1/4 1/4 1/4 1/4
deltaEffectAvg 1/2 -1/2 1/2 -1/2
deltaEffect 1/2 -1/2 -1/2 1/2
deltaControl 0 -1 0 1
```

`intercept`

(\(\gamma\)) is the overall mean while `deltaControl`

(\(\delta_{\mbox{control}}\)) is the difference between the two control arms.

By way of defining `deltaEffectAvg`

(\(\delta_{\mbox{EffectAvg}}\)) we obtain the average treatment effect. Setting `deltaEffect`

(\(\delta_{\mbox{effect}}\)) to half of the difference between the two treatment effects one arrives at the following symmetric relations to obtain the chemo therapy specific treatment effect:

\[\delta_{\mbox{EffectAvg}} = \frac{1}{2} \, [ (\lambda_{\mbox{activeChemoA}} - \lambda_{\mbox{controlChemoA}}) + (\lambda_{\mbox{activeChemoB}} - \lambda_{\mbox{controlChemoB}})]\]

\[\delta_{\mbox{effect}} = \frac{1}{2} \, [ (\lambda_{\mbox{activeChemoA}} - \lambda_{\mbox{controlChemoA}}) - (\lambda_{\mbox{activeChemoB}} - \lambda_{\mbox{controlChemoB}})]\]

\[ \Leftrightarrow \delta_{\mbox{effectChemoA}} = \delta_{\mbox{effectAvg}} + \delta_{\mbox{effect}}\]

\[ \Leftrightarrow \delta_{\mbox{effectChemoB}} = \delta_{\mbox{effectAvg}} - \delta_{\mbox{effect}}\]

`R`

requires to define the inverse relationship, i.e. the mapping from parameters to data strata (as needed to define the design matrix). Therefore, the contrast matrix \(C\) is:

```
<- solve(cc_inv)
cc ::fractions(cc) MASS
```

```
intercept deltaEffectAvg deltaEffect deltaControl
activeChemoA 1 1/2 1 -1/2
controlChemoA 1 -1/2 0 -1/2
activeChemoB 1 1/2 -1 1/2
controlChemoB 1 -1/2 0 1/2
```

In order to setup custom contrasts of discrete explanatory variables like treatment arms, `R`

requires to define the explanatory variable as `factor`

, which has been defined accordingly in the simulation routine of this case study:

`levels(sim$arm)`

`[1] "activeChemoA" "controlChemoA" "activeChemoB" "controlChemoB"`

As `R`

sets up the overall intercept separatley, we need to drop the overall intercept definition from \(C\) and define the contrast of the `factor`

variable to be equal to the matrix \(C\):

`contrasts(sim$arm) <- cc[,-1]`

Then to check if the contrasts work as expected, we can use a dummy fit to see if things get used correctly and the results here are consistent with our definitions:

```
<- lm(y ~ 1 + arm, sim)
lfit1 summary(lfit1)
```

```
Call:
lm(formula = y ~ 1 + arm, data = sim)
Residuals:
Min 1Q Median 3Q Max
-5.506 -2.909 -1.461 1.261 17.064
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.5023 0.3349 13.444 <2e-16 ***
armdeltaEffectAvg 0.9641 0.6698 1.439 0.152
armdeltaEffect -0.4516 0.6698 -0.674 0.501
armdeltaControl 0.6574 0.9526 0.690 0.491
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.229 on 196 degrees of freedom
Multiple R-squared: 0.02433, Adjusted R-squared: 0.009398
F-statistic: 1.629 on 3 and 196 DF, p-value: 0.1839
```

We can see that above model specification for `lm`

estimates the parameter coefficient for the three contrasts `armdeltaEffectAvg`

(\(\delta_{\mbox{EffectAvg}}\)), `armdeltaEffect`

(\(\delta_{\mbox{effect}}\)) and `armdeltaControl`

(\(\delta_{\mbox{control}}\)).

### 14.5.2 Model fit with `brms`

Now we start with analyzing the simulated randomized study and specify the `brms`

model using a Weibull distribution. In a first version, we do not include the literature historical data which requires the definition of a custom design matrix as will be explained below.

First we define with a call to `bf`

the model formula and model family. This defines the model’s likelihood and the parametrization.

`<- bf(y | cens(1-event) ~ 1 + arm, family=weibull()) model_weibull `

The left hand side of the formula, `y | cens(1-event) ~ ...`

, denotes with \(y\) the data being modeled - the time variable -, while `cens`

is a response modifier added with a vertical bar`|`

to the left hand side. Here, `cens`

defines which data rows correspond to (right) censored observations. The censoring information is passed into `cens`

by the variable `1-event`

. The right hand side defines the linear regressor formula used for the mean parameter of the model, which is in this case modeled by default with a \(\log\) link function. The formula `1 + arm`

defines an overall intercept and the covariate `arm`

(coded as `factor`

in the data) is used with it’s custom constrasts as explanatory variable. Finally, the `family`

argument specifies that the event and censoring times stored in \(y\) are distributed according to a `weibull`

probability density distribution.

With the model (and data) being defined, we are left to specify the priors for the model. With the help of the `get_prior`

function we can report which default priors and parameters were instantiated by `brms`

and hence need an assignment of a prior.

`gt(get_prior(model_weibull, sim))`

prior | class | coef | group | resp | dpar | nlpar | lb | ub | source |
---|---|---|---|---|---|---|---|---|---|

b | default | ||||||||

b | armdeltaControl | default | |||||||

b | armdeltaEffect | default | |||||||

b | armdeltaEffectAvg | default | |||||||

student_t(3, 1, 2.5) | Intercept | default | |||||||

gamma(0.01, 0.01) | shape | 0 | default |

Here we define the priors with standard deviations stored in variables, which are passed as `stanvars`

to the `brms`

model. Doing so allows to change the prior quickly, since the Stan model as generated by `brms`

remains unchanged whenever the standard deviation of a prior is changed. A priori it is expected that the outcome for the two control groups will be similar, so the prior on `deltaControl`

(\(\delta_{\mbox{control}}\)) is centered on \(0\). The same is a priori assumed for the `deltaEffect`

(\(\delta_{\mbox{effect}}\)) parameter.

```
<- prior(normal(meanInter, log(4)/1.64),
model_weibull_prior class="Intercept") +
prior(normal(0, sdEffectAvg), coef=armdeltaEffectAvg) +
prior(normal(0, sdDeltaEffect), coef=armdeltaEffect) +
prior(normal(0, sdDeltaControl), coef=armdeltaControl) +
prior(gamma(3, 2.7), class=shape)
```

The prior for the intercept, \(\gamma\), is set to a mean approximately corresponding to an expected median PFS on anti-PD-1 + SoC, as estimated from CheckMate 649 (here 8 month)

`<- stanvar(log(8/log(2)), "meanInter") prior_mean `

The standard deviations are setup with reference to a \(90 \%\) credible interval. Hence, \(log(2)/1.64\) for the standard deviation of the prior for `deltaEffectAvg`

(\(\delta_{\mbox{effectAvg}}\)) means that the event process can increase/decrease in progression speed by a factor of \(2\) at most.

```
<- stanvar(log(2)/1.64, "sdEffectAvg") +
prior_sd stanvar(log(1.25)/1.64, "sdDeltaEffect")+
stanvar(log(1.25)/1.64, "sdDeltaControl")
```

Now the model and prior specifications are complete and we are ready to run the model in `brms`

. First fit the model using the randomized study data

```
<- brm(model_weibull, data=sim, prior=model_weibull_prior,
fit_weibull stanvars=prior_sd + prior_mean,
seed=43534, refresh=0)
```

```
Running MCMC with 4 sequential chains...
Chain 1 finished in 0.5 seconds.
Chain 2 finished in 0.5 seconds.
Chain 3 finished in 0.5 seconds.
Chain 4 finished in 0.5 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.5 seconds.
Total execution time: 2.5 seconds.
```

When looking at the model summary, it is preferable to use `robust=TRUE`

to get the median estimate of the posterior rather than the mean (which would change due to transforms)

`summary(fit_weibull, prob=0.9, robust=TRUE)`

```
Family: weibull
Links: mu = log; shape = identity
Formula: y | cens(1 - event) ~ 1 + arm
Data: sim (Number of observations: 200)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.16 0.12 1.98 2.38 1.00 3998 2848
armdeltaEffectAvg 0.26 0.19 -0.05 0.58 1.00 3806 2833
armdeltaEffect 0.01 0.11 -0.17 0.18 1.00 4514 3328
armdeltaControl 0.06 0.12 -0.13 0.26 1.00 4265 3056
Further Distributional Parameters:
Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
shape 0.97 0.08 0.85 1.10 1.00 4014 3098
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
```

Summarize posterior draws based on point estimates, estimation errors and quantiles. And we transform the estimates to the original scale

```
<- posterior_summary(fit_weibull, prob=c(0.05, 0.95), robust=TRUE)
post_sum print(exp(post_sum[,"Estimate"]), digits=2)
```

```
b_Intercept b_armdeltaEffectAvg b_armdeltaEffect b_armdeltaControl shape Intercept
8.6e+00 1.3e+00 1.0e+00 1.1e+00 2.6e+00 8.5e+00
lprior lp__
1.0e+00 5.4e-140
```

These estimates are the AFT regression coefficients such that the effect of treatment is to *slowdown* the time to an event and thus the treatment effect estimate \(\delta_{\mbox{effectAvg}}\) is positive. To convert to hazard ratios, one needs to multiply with the negative of the estimated shape parameter. Doing so for all 3 regression coefficients by using the `rvar`

facility (provided from the `posterior`

R package) yields:

```
<- as_draws_rvars(fit_weibull)
post_rv <- with(post_rv, exp(-1 * shape * c(b_armdeltaEffectAvg, b_armdeltaEffect, b_armdeltaControl)))
post_rv_HR names(post_rv_HR) <- c("deltaEffectAvg", "deltaEffect", "deltaControl")
print(summarize_draws(post_rv_HR), digits=2)
```

```
# A tibble: 3 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 post_rv_HR[deltaEffectAvg] 0.788 0.777 0.146 0.145 0.568 1.05 1.00 3872. 2876.
2 post_rv_HR[deltaEffect] 1.00 0.995 0.104 0.103 0.839 1.18 1.00 4515. 3209.
3 post_rv_HR[deltaControl] 0.950 0.947 0.110 0.107 0.777 1.13 1.00 4282. 2919.
```

Obtain posterior mean estimates (mean rate). Median of the posterior estimates indicated a roughly 20% increase in survival in `activeChemoA`

compared to `controlChemoA`

.

```
<- data.frame(y=0, event=0, arm=levels(sim$arm))
nd <- posterior_epred(fit_weibull, newdata=nd)
post_mean colnames(post_mean) <- nd$arm
quantile(post_mean[,"activeChemoA"] / post_mean[,"controlChemoA"],
probs=c(0.05, 0.5, 0.95))
```

```
5% 50% 95%
0.9425351 1.3043160 1.8179356
```

Now extract posterior coefficients with a 90% credible interval

```
<- fixef(fit_weibull, prob=c(0.05, 0.95), robust=TRUE)
post_est gt(as.data.frame(post_est)) %>% fmt_number(everything(), decimals=2)
```

Estimate | Est.Error | Q5 | Q95 |
---|---|---|---|

2.16 | 0.12 | 1.98 | 2.38 |

0.26 | 0.19 | −0.05 | 0.58 |

0.01 | 0.11 | −0.17 | 0.18 |

0.06 | 0.12 | −0.13 | 0.26 |

Define trial successfull if the `deltaEffectAvg`

exceed 0 (which means the survival time is increased) and the result is not significant here.

`"armdeltaEffectAvg","Q5"] > 0 post_est[`

`[1] FALSE`

### 14.5.3 Incoporate simulated historical data

Since we have the simulated historical data `hist1`

on one control arm, we can incorporate the information by adding the corresponding covariate `hist1`

(indicator variable being \(1\) for the historical study and \(0\) otherwise) when specifying the formula for meta-analytic-combined (MAC) analysis:

`<- bf(y | cens(1-event) ~ 1 + arm + hist1, family=weibull()) model_mac_weibull `

And `get_prior`

function shows that in addition to the priors for the weibull model on randomised data, the co-data model needs one more prior on the historical data.

`gt(get_prior(model_mac_weibull, sim))`

prior | class | coef | group | resp | dpar | nlpar | lb | ub | source |
---|---|---|---|---|---|---|---|---|---|

b | default | ||||||||

b | armdeltaControl | default | |||||||

b | armdeltaEffect | default | |||||||

b | armdeltaEffectAvg | default | |||||||

b | hist1 | default | |||||||

student_t(3, 1, 2.5) | Intercept | default | |||||||

gamma(0.01, 0.01) | shape | 0 | default |

We use a student with 6 degrees of freedom for the regression coefficient for the historical data prior. Using a Student-t distribution robustifies against possible data discordance between the randomised study and the historical data, since very large deviations between the two data sets are admitted by the heavy tails of the prior. The prior is centered at a mean of \(0\), corresponding to no expected difference in the outcome for historical and current study.

```
<- model_weibull_prior +
model_weibull_mac_prior prior(student_t(6, 0, sdHist), coef=hist1)
```

A variance allowing substantial differences, a \(1.8\) fold of increase/decrease, between the randomised study and historical data is set.

`<- stanvar(log(1.8)/1.64, "sdHist") prior_sdHist `

Now we combine the simulated trial data and historical data to fit the co-data MAC model. Don’t forget to set contrasts of arm in the combined data to the contrast matrix.

```
<- bind_rows(sim, hdata1)
comb_data contrasts(comb_data$arm) <- cc[,-1]
<- brm(model_mac_weibull, data=comb_data,
fit_mac_weibull prior=model_weibull_mac_prior,
stanvars=prior_sd + prior_mean + prior_sdHist,
seed=43534, refresh=0)
```

```
Running MCMC with 4 sequential chains...
Chain 1 finished in 1.7 seconds.
Chain 2 finished in 1.7 seconds.
Chain 3 finished in 1.6 seconds.
Chain 4 finished in 1.7 seconds.
All 4 chains finished successfully.
Mean chain execution time: 1.7 seconds.
Total execution time: 7.1 seconds.
```

And we can compare the results from two models. Inclusion of historical data on `controlChemoA`

affect the estimates of coefficients for overall mean and `EffectAverage`

.

`summary(fit_weibull, prob=0.9, robust=TRUE)`

```
Family: weibull
Links: mu = log; shape = identity
Formula: y | cens(1 - event) ~ 1 + arm
Data: sim (Number of observations: 200)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.16 0.12 1.98 2.38 1.00 3998 2848
armdeltaEffectAvg 0.26 0.19 -0.05 0.58 1.00 3806 2833
armdeltaEffect 0.01 0.11 -0.17 0.18 1.00 4514 3328
armdeltaControl 0.06 0.12 -0.13 0.26 1.00 4265 3056
Further Distributional Parameters:
Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
shape 0.97 0.08 0.85 1.10 1.00 4014 3098
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
```

`summary(fit_mac_weibull, prob=0.9, robust=TRUE)`

```
Family: weibull
Links: mu = log; shape = identity
Formula: y | cens(1 - event) ~ 1 + arm + hist1
Data: comb_data (Number of observations: 600)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.12 0.10 1.96 2.30 1.00 4909 3769
armdeltaEffectAvg 0.29 0.18 -0.01 0.59 1.00 4051 2881
armdeltaEffect 0.01 0.11 -0.17 0.18 1.00 4748 3261
armdeltaControl 0.07 0.12 -0.13 0.27 1.00 4343 2769
hist1 -0.20 0.14 -0.43 0.02 1.00 4039 2944
Further Distributional Parameters:
Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
shape 1.00 0.04 0.94 1.07 1.00 4712 3074
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
```

Next, we check the informativeness of the historical data by only including the historical data. Note that we have to request that factor levels for which no realization is present in the data set must be kept, since only the `controlChemoA`

arm is present in the historical data.

```
contrasts(hdata1$arm) <- cc[,-1]
<- brm(model_mac_weibull, data=hdata1,
fit_mac_prior_weibull prior=model_weibull_mac_prior,
drop_unused_levels=FALSE,
stanvars=prior_sd + prior_mean + prior_sdHist,
seed=43534, refresh=0)
```

```
Running MCMC with 4 sequential chains...
Chain 1 finished in 1.1 seconds.
Chain 2 finished in 1.1 seconds.
Chain 3 finished in 1.1 seconds.
Chain 4 finished in 1.1 seconds.
All 4 chains finished successfully.
Mean chain execution time: 1.1 seconds.
Total execution time: 4.8 seconds.
```

`gt(as.data.frame(fixef(fit_mac_prior_weibull, robust=TRUE))) %>% fmt_number(everything(), decimals=2)`

Estimate | Est.Error | Q2.5 | Q97.5 |
---|---|---|---|

1.73 | 0.47 | 0.76 | 2.69 |

0.01 | 0.42 | −0.80 | 0.85 |

0.00 | 0.14 | −0.28 | 0.27 |

0.00 | 0.14 | −0.27 | 0.27 |

0.01 | 0.38 | −0.86 | 0.87 |

We see that we essentially sample the prior for most regression coefficients, but we do substantially reduce the standard error of the intercept estimate which is almost halved in comparison to the prior standard deviation of 0.85.

### 14.5.4 Including historical data of average SoC

So far we have setup a custom parametrization allowing to define priors aligned with prior expectations (small difference of the treatment effect and small difference between SoC). However, the CheckMate 649 historical data reports the effect of an average SoC treatment. Thus, the data of CheckMate 649 is incompatible with the 4 defined categorical factor levels so far as it does not correspond to any of the 4 treatment arms in the current data set. However, as the custom factor contrasts and `factor`

level definitions are translated by `R`

into a design matrix, we can customize the created design matrix directly. This requires to setup the design matrix manually giving us the freedom to define the design matrix appropiatley for the CheckMate 649 data set.

There is an alternative way to use the customized contrast and setup analysis by using respective indicator flags, which is created using `model.matrix`

. The `model.matrix`

function is used internally by `brms`

to create the design matrix for the regression problem

```
<- model.matrix(~ 1 + arm, comb_data)
Xind gt_preview(Xind)
```

(Intercept) | armdeltaEffectAvg | armdeltaEffect | armdeltaControl | |
---|---|---|---|---|

1 | 1 | -0.5 | 0 | -0.5 |

2 | 1 | 0.5 | 1 | -0.5 |

3 | 1 | -0.5 | 0 | -0.5 |

4 | 1 | 0.5 | 1 | -0.5 |

5 | 1 | -0.5 | 0 | -0.5 |

6..599 | ||||

600 | 1 | -0.5 | 0 | -0.5 |

By adding the generated design matrix to the original data set, we can then simply delcare the model formula using these created indicator variables:

`<- cbind(comb_data, Xind[,-1]) comb_data_ind `

Thus we use those indicator variables in the fit as if these were continuous covariates which replaces the use of the categorical `factor`

variable `arm`

used in the first approach.

```
<- bf(y | cens(1-event) ~ 1 + armdeltaEffectAvg +
model_mac_weibull_ind + armdeltaControl + hist1,
armdeltaEffect family=weibull())
```

And query what prior parameters need to be set

`gt(get_prior(model_mac_weibull_ind, comb_data_ind))`

prior | class | coef | group | resp | dpar | nlpar | lb | ub | source |
---|---|---|---|---|---|---|---|---|---|

b | default | ||||||||

b | armdeltaControl | default | |||||||

b | armdeltaEffect | default | |||||||

b | armdeltaEffectAvg | default | |||||||

b | hist1 | default | |||||||

student_t(3, 1, 2.5) | Intercept | default | |||||||

gamma(0.01, 0.01) | shape | 0 | default |

Fit the model using indicator flags with the same priors as before and compare with the previous one fitted with `arm`

variable, the results are the same.

```
<- prior(normal(meanInter, log(4)/1.64),
model_weibull_mac_prior_ind class="Intercept") +
prior(normal(0, sdEffectAvg), coef=armdeltaEffectAvg) +
prior(normal(0, sdDeltaEffect), coef=armdeltaEffect) +
prior(normal(0, sdDeltaControl), coef=armdeltaControl) +
prior(normal(0, sdHist), coef=hist1)
<- brm(model_mac_weibull_ind, data=comb_data_ind,
fit_mac_weibull_ind prior=model_weibull_mac_prior_ind,
stanvars=prior_sd + prior_mean + prior_sdHist,
seed=43534, refresh=0)
```

```
Running MCMC with 4 sequential chains...
Chain 1 finished in 1.7 seconds.
Chain 2 finished in 1.6 seconds.
Chain 3 finished in 1.7 seconds.
Chain 4 finished in 1.6 seconds.
All 4 chains finished successfully.
Mean chain execution time: 1.7 seconds.
Total execution time: 7.0 seconds.
```

`summary(fit_mac_weibull, prob=0.9, robust=TRUE)`

```
Family: weibull
Links: mu = log; shape = identity
Formula: y | cens(1 - event) ~ 1 + arm + hist1
Data: comb_data (Number of observations: 600)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.12 0.10 1.96 2.30 1.00 4909 3769
armdeltaEffectAvg 0.29 0.18 -0.01 0.59 1.00 4051 2881
armdeltaEffect 0.01 0.11 -0.17 0.18 1.00 4748 3261
armdeltaControl 0.07 0.12 -0.13 0.27 1.00 4343 2769
hist1 -0.20 0.14 -0.43 0.02 1.00 4039 2944
Further Distributional Parameters:
Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
shape 1.00 0.04 0.94 1.07 1.00 4712 3074
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
```

`summary(fit_mac_weibull_ind, prob=0.9, robust=TRUE)`

```
Family: weibull
Links: mu = log; shape = identity
Formula: y | cens(1 - event) ~ 1 + armdeltaEffectAvg + armdeltaEffect + armdeltaControl + hist1
Data: comb_data_ind (Number of observations: 600)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.12 0.10 1.96 2.29 1.00 5075 3553
armdeltaEffectAvg 0.28 0.18 0.00 0.57 1.00 4184 3537
armdeltaEffect 0.01 0.10 -0.17 0.17 1.00 4688 3177
armdeltaControl 0.07 0.11 -0.12 0.25 1.00 5187 3251
hist1 -0.19 0.14 -0.43 0.02 1.00 4334 3137
Further Distributional Parameters:
Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
shape 1.00 0.04 0.93 1.07 1.00 4680 2946
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
```

Note that the equality holds exactly for the entire posterior, since the very same problem (using the same random seed) has been fitted by `brms`

:

```
<- as_draws_matrix(fit_mac_weibull)
post_mac <- as_draws_matrix(fit_mac_weibull_ind)
post_mac_ind all(post_mac == post_mac_ind)
```

`[1] FALSE`

The benefit of this approach is that we can then include external data corresponding to the `controlAverage`

from the second historical study CheckMate 649, which is not one of the treatment groups but the average effect of control groups, \(\frac{1}{2} \, (\lambda_{\mbox{controlChemoA}}+\lambda_{\mbox{controlChemoB}})\).

To now find the indicators needed for the effect of `controlAverage`

in the given parametrization so far, we can obtain this via simple algebra as:

```
<- MASS::fractions(solve(cc_inv))
cc cc
```

```
intercept deltaEffectAvg deltaEffect deltaControl
activeChemoA 1 1/2 1 -1/2
controlChemoA 1 -1/2 0 -1/2
activeChemoB 1 1/2 -1 1/2
controlChemoB 1 -1/2 0 1/2
```

`0.5*(cc["controlChemoA",] + cc["controlChemoB",])`

```
intercept deltaEffectAvg deltaEffect deltaControl
1 -1/2 0 0
```

Hence we must set the `deltaEffectAvg`

column to \(-\frac{1}{2}\) and the others to \(0\) to get the indicator for `controlAverage`

(the overall intercept is always added to the linear predictor as we define `1`

as the first term of the model formula) in order to include the new historical data from CheckMate649.

```
<- rename(hdata2, y=time, event=status) %>%
hdata2_ind mutate(hist1=0, hist2=1, armdeltaEffectAvg=-0.5,
armdeltaEffect=0, armdeltaControl=0)
<- bind_rows(comb_data_ind, hdata2_ind) all_comb_data_ind
```

To include the second historical data, we specify the formula and the priors as we did for one historical study.

```
<- bf(y | cens(1-event) ~ 1 + armdeltaEffectAvg +
model_mac_weibull_ind_all + armdeltaControl + hist1+ hist2,
armdeltaEffect family=weibull())
<- model_weibull_mac_prior +
model_weibull_mac_prior_all prior(normal(0, sdHist), coef=hist2)
```

Fit the model and compare:

```
<- brm(model_mac_weibull_ind_all, data=all_comb_data_ind,
fit_all_mac_weibull_ind prior=model_weibull_mac_prior_all,
stanvars=prior_sd + prior_mean + prior_sdHist,
seed=43534, refresh=0)
```

```
Running MCMC with 4 sequential chains...
Chain 1 finished in 6.0 seconds.
Chain 2 finished in 6.2 seconds.
Chain 3 finished in 6.0 seconds.
Chain 4 finished in 5.8 seconds.
All 4 chains finished successfully.
Mean chain execution time: 6.0 seconds.
Total execution time: 24.3 seconds.
```

` fit_mac_weibull_ind`

```
Family: weibull
Links: mu = log; shape = identity
Formula: y | cens(1 - event) ~ 1 + armdeltaEffectAvg + armdeltaEffect + armdeltaControl + hist1
Data: comb_data_ind (Number of observations: 600)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.12 0.10 1.93 2.33 1.00 5075 3553
armdeltaEffectAvg 0.29 0.17 -0.05 0.63 1.00 4184 3537
armdeltaEffect 0.01 0.10 -0.20 0.21 1.00 4688 3177
armdeltaControl 0.07 0.11 -0.16 0.29 1.00 5187 3251
hist1 -0.20 0.14 -0.47 0.06 1.00 4334 3137
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 1.00 0.04 0.92 1.09 1.00 4680 2946
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
```

` fit_all_mac_weibull_ind`

```
Family: weibull
Links: mu = log; shape = identity
Formula: y | cens(1 - event) ~ 1 + armdeltaEffectAvg + armdeltaEffect + armdeltaControl + hist1 + hist2
Data: all_comb_data_ind (Number of observations: 1389)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.09 0.09 1.92 2.27 1.00 4002 3141
armdeltaEffectAvg 0.23 0.16 -0.09 0.55 1.00 2743 2920
armdeltaEffect 0.00 0.10 -0.19 0.20 1.00 3616 2884
armdeltaControl 0.09 0.12 -0.14 0.32 1.00 3555 3039
hist1 -0.23 0.13 -0.49 0.00 1.00 2287 2849
hist2 0.48 0.12 0.24 0.71 1.00 2108 1868
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 1.12 0.03 1.07 1.18 1.00 4034 3292
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
```

## 14.6 Conclusion

Here we showed how `brms`

enabled us to model time-to-event data with censoring and incorporate historical data into the analysis. By use of custom contrasts we controlled the parametrization of the problem allowing for the specification of priors aligned with the prior knowledge. This is possible in `R`

using the `contrasts`

function to setup custom contrast matrices used to setup the design matrix of the model. Going beyond that we also discussed that in some situations the facilities in `R`

are too limiting and a customized setup of the model design matrix allows for even greater flexibility. This allowed in this case study the inclusion of historical data which corresponds to the average effect of standard of care.

## 14.7 Exercises

In

`brms`

various families can be used to model time-to-event data. Re-fit the model using`exponential`

distribution and compare with the`weibull`

results.Cox regression is implemented in

`brms`

using M-splines, which is a very close approximation to the semi-parametric Cox model. Show that the estimates for the covariate effect we get from`brms`

are quite similar to what we get from`survival`

package. Also compare the estimate of the overall intercept which needs to account for an offset as detailled in section @ref(sec:tte-model-brms).