```
library(ggplot2)
library(dplyr)
library(knitr)
library(brms)
library(posterior)
library(bayesplot)
library(RBesT)
library(GGally)
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)
set.seed(57339)
# allow for wider and prettier printing
options(width=120, digits=2)
```

# 6 Use of historical control data with stratification

Here we will demonstrate the use of historical control data as an example for a meta-analytic predictive (MAP) prior approach for the case of data being divided intro strata which are modelled using covariates.

This case study demonstrates

- setting up a random effect meta-analysis with covariates for data strata
- how to summarize with a multi-variate normal mixture the MAP prior with covariates
- how to use the multi-variate normal mixture MAP prior as a prior for the main analysis using
`brms`

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

## 6.1 Background

Assume we wish to evaluate in a Phase III study, the effectiveness of an investigational treatment by comparing it to a placebo in a population with a rare disease. The study focused on evaluating the treatment’s efficacy in different age groups from the age range 2 to under 18 years old. The primary endpoint of the study is the change in a Patient Reported Outcome (PRO) total score at the end of the follow-up period (12 months).

The statistical hypothesis being tested for the primary endpoint is that there is no difference in the change from baseline in the PRO total score between the treatment group and the placebo group.

A Bayesian linear regression model will be applied, adjusting for the age classified into 3 groups (2 to under 5, 5 to under 13, and 13 to under 18) and baseline score. The age group of 5 to under 13 years old will serve as the reference category, allowing the use of priors on the intercept and coefficients for treatment to represent the mean change in the score at 12 months in this subgroup of patients.

In order to inform the model, informative priors will be included for the mean change in the PRO score at 12 months in each age subgroup. These priors will be derived from historical studies.

## 6.2 Data

At the design stage of the trial, data from two historical studies were available for the control group.

- Study 1: A natural history study of patients with the disease.
- Study 2: Control arm of a randomized study evaluating an alternative therapy.

These historical studies provide valuable information for the primary evaluation of the treatment effect in the trial of interest. Casting the historical data into an informative prior on the control response is expected to increase the precision in the estimation of the treatment effect for the Phase III study of interest.

The data set obtained from these two studies are summary statistics for the PRO score at baseline (BASE), the mean change from baseline (CHGm), and is standard error (CHGse) stratified by age group (AGEGRP).

```
<- tibble::tribble(
hdata ~Study , ~AGEGRPN, ~AGEGRP, ~n, ~BASE, ~CHGm, ~CHGse,
"1. STUDY 1" , "1" , "Age >=2 to <5" , 28, 16.8, -0.7, 1.00,
"1. STUDY 1" , "2" , "Age >=5 to <13" , 50, 26.7, -1.1, 0.73,
"1. STUDY 1" , "3" , "Age >=13 to <18" , 23, 18.5, -1.1, 1.27,
"2. STUDY 2" , "1" , "Age >=2 to <5" , 42, 19.9, 0.2, 0.66)
<- with(hdata, round(weighted.mean(BASE, n), 2))
base_mean <- hdata |>
hdata mutate(
# make the age group a factor and set the middle age group to be the reference
age_group = factor(AGEGRPN, c("2", "1", "3")),
# add a centered baseline covariate
cBASE = BASE - base_mean,
# back-calculate the sd
CHGsd = CHGse * sqrt(n),
# make Study a factor for better handling
Study = factor(Study),
# make a variable 'label' for later use in plots
label = paste0(Study, "/", age_group),
label2 = paste0(Study, "/", age_group, ". ", AGEGRP),
# include active treatment indicator (in this case, is set to 0 as all patients are 'control')
active = 0L)
|>
hdata select(Study, AGEGRP, n, BASE, cBASE, CHGm, CHGse, CHGsd) |>
kable()
```

Study | AGEGRP | n | BASE | cBASE | CHGm | CHGse | CHGsd |
---|---|---|---|---|---|---|---|

1. STUDY 1 | Age >=2 to <5 | 28 | 17 | -4.7 | -0.7 | 1.00 | 5.3 |

1. STUDY 1 | Age >=5 to <13 | 50 | 27 | 5.2 | -1.1 | 0.73 | 5.2 |

1. STUDY 1 | Age >=13 to <18 | 23 | 18 | -3.0 | -1.1 | 1.27 | 6.1 |

2. STUDY 2 | Age >=2 to <5 | 42 | 20 | -1.6 | 0.2 | 0.66 | 4.3 |

## 6.3 Model description

The primary endpoint will be analyzed using a Bayesian regression model with the mean change from baseline in the PRO score as the response \(Y\) and with treatment \(Z\) (\(Z = 1\) is treatment, \(Z = 0\) is control). The basic analysis model for the Phase III model can be written as

\[Y|\mu,\sigma \sim \mbox{Normal}(\mu, \sigma^2)\]

\[\mu = \beta_0 + \beta_1 \, x_1 + \beta_2 \, x_2 + \beta_3 \, cBASE + \theta \, Z\]

where:

- \(\beta_0\) is the mean change from baseline in the score on the control arm for the 5 to 13 years old age group,
- \(x_1,x_2\) are indicators for the Age >=2 to <-5 and Age >=13 to <18 groups,
- \(\beta_3\) is the covariate for the centered baseline score,
- \(\theta\) is the treatment effect comparing mean change from baseline for the treatment arm vs control arm,
- \(\sigma\) is the sample standard deviation.

Inference for the primary analysis focuses on the posterior distributions of the treatment effect comparing mean change from baseline in the PRO score on the treatment arm vs control arm, \(\theta\).

To form priors, the mean change from baseline as well the standard errors from each historical study and age subgroup are combined using the meta-analytic predictive approach (MAP, Neuenschwander et al 2010). While the MAP approach is mostly used for an intercept only model, the use of covariates complicates the analysis here. However, this complication can be resolved by considering the equivalent meta-analytic (MAC) combined approach. Recall, that whereas the MAP approach is a 2-step procedure (derived MAP prior from historical data, then apply to current data), the MAC approach performs a joint analysis of the historical data and the current data. Both approaches result in exactly the same results (Schmidli et al 2014). Thus, it is useful to first consider the respective MAC model. For this we expand the basic model shown above to also include a random intercept effect, which accounts for between-study heterogeneity. Denoting with \(i\) the study, we thus model the data as

\[Y_i|\mu_i,\sigma \sim \mbox{Normal}(\mu_i, \sigma^2)\]

\[\mu_i = b_i + \beta_1 \, x_1 + \beta_2 \, x_2 + \beta_3 \, cBASE + \theta \, Z\]

\[b_i|\beta_0,\tau \sim \mbox{Normal}(\beta_0, \tau^2).\]

The key assumption of this model is that the covariate effects for age, centered baseline score and treatment effect are *pooled* accross all studies. That is, no between-study heterogeneity is admitted in these effects; heterogeneity is only accounted for in the intercept term representing the mean change from baseline for control of the reference age group with a baseline value as used for centering (21.45). Note that since we do only include data on control for the historical studies, we will infer the treatment effect \(\theta\) for the new Phase III study without any influence from the other studies.

However, in practice we will still want to run the primary analysis in a 2-step MAP approach. The desire here is to be able to evaluate the properties of the MAP prior derived from the historical data in advance. The challenge becomes then that we have to use the information from fitting the above model to the historical data for the *entire model posterior* and not merley the intercept only. Therefore, the MAP prior becomes multi-variate, since the MAP prior includes the covariate effects.

Hence, the model is fitted first to the historical data under the assumption of a known sampling standard deviation (accounting for the reported standard errors).

As priors for \(\beta_0,\beta_1,\beta_2,\beta_3,\theta\) weakly informative \(\mbox{Normal}(0,s^2)\) distributions are used, with \(s = 5\). The choice of \(5\) follows from the observation that the historical data has approximatly a sampling standard deviation of \(5\) (the pooled estimate from the historical data is 5.11 ). For the between-trial standard deviation \(\tau\) a half-normal prior with scale \(s / 4\) is used. This is considered a plausible conservative, weakly-informative prior, as values of the ratio \(\tau/\sigma \ge 1/4\) are generally considered to represent substantial between-trial variability (Neuenschwander et al 2010) and the proposed prior will generally place more than half of the prior probability on this area.

## 6.4 Implementation

Using `brms`

we now specify the MAP model step by step. We first define the model:

```
# Set up a model formula for use in 'brms'
<- bf(CHGm | se(CHGse) ~ 1 + age_group + cBASE + active + (1 | Study),
hist_model family=gaussian, center=FALSE)
# Display prior model across all parameters
get_prior(hist_model, hdata)
```

```
prior class coef group resp dpar nlpar lb ub source
(flat) b default
(flat) b active (vectorized)
(flat) b age_group1 (vectorized)
(flat) b age_group3 (vectorized)
(flat) b cBASE (vectorized)
(flat) b Intercept (vectorized)
student_t(3, 0, 2.5) sd 0 default
student_t(3, 0, 2.5) sd Study 0 (vectorized)
student_t(3, 0, 2.5) sd Intercept Study 0 (vectorized)
```

**Note** that we do set explicitly the option `center=FALSE`

in the function call to `bf`

. This avoids the otherwise automatic centering done by `brms`

, which is needed here as we require full control over the parametrization of the model.

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

```
## Fit historical data with between-study heterogeneity -----------------------
# For β_0,β_1,β_2,β_3,θ weakly informative N(0,5^2) distributions are used.
# For the group level effect of the study, we use τ^2=(5.0/4)^2.
# Set then sd = 5 for later use
<- 5.0
sd <- sd/4
sd_tau <- prior_string(glue::glue("normal(0, {sd})"), class="b") +
hist_prior prior_string(glue::glue("normal(0, {sd_tau})"), class="sd", coef="Intercept", group="Study")
```

Now we are ready to run the model in `brms`

:

```
# Apply brms model using historical data hdata and the model defined in hist_model
<- brm(hist_model,
map_mc_brms data = hdata,
prior = hist_prior,
seed = 234324, refresh = 0,
control = list(adapt_delta = 0.99))
```

`Start sampling`

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

` map_mc_brms`

```
Family: gaussian
Links: mu = identity; sigma = identity
Formula: CHGm | se(CHGse) ~ 1 + age_group + cBASE + active + (1 | Study)
Data: hdata (Number of observations: 4)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~Study (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.80 0.63 0.03 2.33 1.00 1452 969
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1.13 1.92 -4.95 2.77 1.00 1208 1621
age_group1 1.13 2.71 -4.40 6.44 1.00 1225 1611
age_group3 0.39 2.64 -4.73 5.71 1.00 1369 1837
cBASE 0.06 0.31 -0.55 0.69 1.00 1258 1761
active -0.08 4.94 -9.91 9.51 1.00 3201 2483
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.00 0.00 0.00 0.00 NA NA NA
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).
```

## 6.5 Results

Once the MAP prior is obtained in MCMC form a model check for it is recommended. In `brms`

model diagnostic functions are directly available and essentially expose the functionality found in the `bayesplot`

R package. A suitable `bayesplot`

plot in this situation could be an `intervals`

plot as:

```
## Posterior predictive check for the mean change in Score at 12mo by age group -------
pp_check(map_mc_brms, type="intervals") +
scale_x_continuous("Study / Strata",
breaks = 1:nrow(hdata),
labels = hdata$label2,
trans = "reverse") +
ylab("Mean change from baseline") +
xlab(NULL) +
coord_flip() +
hline_0(linetype=2) +
ggtitle("Posterior predictive check") +
theme(plot.background = element_rect(fill = "white"),
legend.position="right",
# suppress vertical grid lines for better readability of intervals
panel.grid.major.y = element_blank())
```

`Using all posterior draws for ppc type 'intervals' by default.`

The call of the `pp_check`

function is forwarded to the respective `ppc_*`

functions for posterior predictive checks from `bayesplot`

(depending on the `type`

argument). The plots are designed to compare the *posterior predictive* distribution to the observed data rather than comparing mean estimates to one another. Thus, the outcome of each study/age group in the historical data set is sampled according to the fitted model and the resulting predictive distribution of the outcome (change in score) is compared to the observed outcome. The `intervals`

predictive probability check summarises the predictive distributions using a light color for an outer credible interval range and a darker line for an inner credible interval. The outer defaults to a 90% credible interval (`prob_outer`

argument) while the inner uses a 50% credible interval (`prob`

argument). The light dot in the middle is the median of the predictive distribution and the dark dot is the outcome \(y\). As we can observe, the outcomes \(y\) for all age groups in historical studies all are contained within outer credible intervals of the predictive distributions for the simulated replicate data \(y_{rep}\). We can conclude that the model is consistent with the data.

Once the model has been checked for plausibility, we can proceed and derive the main target of the MAP analysis, which is the MAP prior in *parametric* form. The MAP priors are stored in the `post_coefs_mix`

object. Note that since we will be predicting for a *new* study, we sample a random intercept using the between study variation.

```
# For the simple case of this case study for which only a random
# effect is assigned to the intercept, we can use the brms prediction
# functions. We are seeking a sample of the intercept with
# between-trial heterogeneity, but without having other terms
# included in the linear predictor. Hence, we can predict the mean
# change from baseline for a new study arm which is has the covariate
# values set to the reference
<- tibble(Study="new_study", age_group="2", cBASE=0, CHGse=0, active=0L)
new_study_ref <- rvar(posterior_epred(map_mc_brms,
rv_study_sim_inter newdata=new_study_ref,
allow_new_levels=TRUE,
sample_new_levels="gaussian"))
rv_study_sim_inter
```

```
rvar<4000>[1] mean ± sd:
[1] -1.1 ± 2.2
```

```
# Get random draws from the posterior of each parameter
# after fitting historical data
<- as_draws_rvars(map_mc_brms)
hist_post_draws
# Whenever one also models the covariate effects to be varying between
# studies, one needs to smaple the random effect directly, since brms
# always simulates the linear predictor by convention. In this case
# one needs to sample the random intercept directly from the posterior.
# Add variability to the intercept to capture random between study variation
<- with(hist_post_draws,
rv_study_sim_inter_alt rvar_rng(rnorm, 1, b_Intercept, sd_Study__Intercept, ndraws=ndraws(hist_post_draws)))
# Get draws for the fixed effects from the model posterior
<- subset_draws(hist_post_draws,
hist_post_coefs variable="^b_",
regex=TRUE) |>
as_draws_matrix()
# Replace the intercept draws with the intercept and added variability
"b_Intercept"] <- as_draws_matrix(rv_study_sim_inter)[,1]
hist_post_coefs[,
# Get a descriptive statistics for posterior draws by coefficient
summarize_draws(hist_post_coefs) |> kable()
```

variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
---|---|---|---|---|---|---|---|---|---|

b_Intercept | -1.11 | -1.12 | 2.16 | 1.98 | -4.56 | 2.53 | 1 | 1395 | 2185 |

b_age_group1 | 1.13 | 1.15 | 2.71 | 2.64 | -3.34 | 5.54 | 1 | 1225 | 1611 |

b_age_group3 | 0.39 | 0.36 | 2.64 | 2.59 | -3.92 | 4.79 | 1 | 1369 | 1837 |

b_cBASE | 0.06 | 0.06 | 0.31 | 0.31 | -0.46 | 0.57 | 1 | 1258 | 1761 |

b_active | -0.08 | -0.04 | 4.94 | 4.86 | -8.19 | 7.98 | 1 | 3201 | 2483 |

```
# Summarize as mixture multi-variate normal MAP distribution
<- mixfit(sample = hist_post_coefs, type = "mvnorm", Nc = 3)
post_coefs_mix print(post_coefs_mix)
```

```
EM for Multivariate Normal Mixture Model
Log-Likelihood = -31931
Multivariate normal mixture
Outcome dimension: 5
Mixture Components:
comp1 comp2 comp3
w 0.4068 0.3653 0.2279
m[b_Intercept] -1.6714 -1.0509 -0.1817
m[b_age_group1] 1.9530 0.9640 -0.0592
m[b_age_group3] 1.0223 0.1458 -0.3458
m[b_cBASE] 0.1528 0.0212 -0.0556
m[b_active] 0.2407 -0.3209 -0.2476
s[b_Intercept] 1.9893 1.5619 2.8476
s[b_age_group1] 2.6607 2.2126 2.9926
s[b_age_group3] 2.7526 2.4530 2.4833
s[b_cBASE] 0.3133 0.2761 0.3143
s[b_active] 5.1088 4.9198 4.6240
rho[b_age_group1,b_Intercept] -0.8653 -0.9432 -0.5990
rho[b_age_group3,b_Intercept] -0.7765 -0.7996 -0.5198
rho[b_cBASE,b_Intercept] -0.8167 -0.8602 -0.5286
rho[b_active,b_Intercept] 0.0156 0.0289 -0.0121
rho[b_age_group3,b_age_group1] 0.8778 0.7894 0.8255
rho[b_cBASE,b_age_group1] 0.9196 0.9125 0.9271
rho[b_active,b_age_group1] -0.0151 -0.0144 -0.0285
rho[b_cBASE,b_age_group3] 0.8637 0.8029 0.8447
rho[b_active,b_age_group3] 0.0317 -0.0054 -0.0147
rho[b_active,b_cBASE] 0.0027 -0.0248 -0.0355
```

```
# notice the considerable correlations present in the posterior
round(cor(hist_post_coefs), 2)
```

```
b_Intercept b_age_group1 b_age_group3 b_cBASE b_active
b_Intercept 1.00 -0.80 -0.70 -0.74 0.00
b_age_group1 -0.80 1.00 0.84 0.92 0.00
b_age_group3 -0.70 0.84 1.00 0.85 0.02
b_cBASE -0.74 0.92 0.85 1.00 0.00
b_active 0.00 0.00 0.02 0.00 1.00
```

And a comparison of the fitted density vs the histogram of the MCMC sample is available as:

`plot(post_coefs_mix)`

```
Diagnostic plots for mixture multivariate normal densities are experimental.
Please note that these are subject to changes in future releases.
```

`$mixpairs`

```
# Create summary statistics for the posteriors draws of coefficients
<- as_draws_df(hist_post_coefs)
hist_post_coefs_df 1:5] hist_post_coefs_df[,
```

`Warning: Dropping 'draws_df' class as required metadata was removed.`

```
# A tibble: 4,000 × 5
b_Intercept b_age_group1 b_age_group3 b_cBASE b_active
<dbl> <dbl> <dbl> <dbl> <dbl>
1 -2.80 3.09 2.93 0.424 -5.78
2 -1.61 -1.35 -0.711 0.154 2.07
3 -0.352 0.419 -3.26 -0.0461 -3.52
4 1.97 -2.83 -2.24 -0.465 1.13
5 -0.154 -1.13 -1.40 -0.116 -0.732
6 -2.38 4.09 3.57 0.346 1.80
7 -3.08 4.67 2.69 0.145 -4.13
8 0.652 -0.224 -0.208 0.303 -7.01
9 -0.913 1.68 0.173 0.225 -7.50
10 -0.900 0.535 0.821 0.0112 -1.57
# ℹ 3,990 more rows
```

`ggpairs(hist_post_coefs_df[1:1000,1:5])`

`Warning: Dropping 'draws_df' class as required metadata was removed.`

These MAP priors summaries for the coefficients of the model are displayed in the table below. Note that since the historical data did not contain data in both the treatment and control arms, the parameter related to the treatment effect (`b_active`

) will remain with a weakly informative prior.

```
<- hist_post_coefs_df |>
map_coef_table summarise_draws(median, ~quantile2(.x, probs = c(0.025, 0.975))) |>
rename(lower = q2.5,
upper = q97.5)
```

Often it is of interest to evaluate the properties of the MAP prior at the design stage. This is complicated in this example due to the requirement to assume a certain baseline value for the used covariates. However, note that these evaluations can be repeated once the trial has enrolled the study population given that only baseline data is needed for this. Here we assume a baseline score of 20:

```
<- c("Age >=2 to <5",
agegroups "Age >=5 to <13",
"Age >=13 to <18")
# Set a value for the baseline score
<- 20
base_score_new # Create summary statistics for the posteriors draws of means across subgroups
<- tibble(Study="new_study",
new_study_base20 age_group=c("1", "2", "3"),
cBASE=base_score_new-base_mean,
CHGse=0,
active=0L)
<- new_study_base20 |>
map_resp_table ::add_epred_rvars(map_mc_brms, allow_new_levels=TRUE, sample_new_levels="gaussian", value="MAP") |>
tidybayesmutate(summarise_draws(MAP, ~quantile2(.x, probs = c(0.025, 0.5, 0.975))),
MAP_ess=ess(mixfit(draws_of(MAP)[,1], Nc=3), sigma=5),
treatment = "Control") |>
rename(median= q50,
lower = q2.5,
upper = q97.5) |>
select(MAP, treatment, MAP_ess, median, lower, upper)
|> kable() map_resp_table
```

MAP | treatment | MAP_ess | median | lower | upper |
---|---|---|---|---|---|

-0.088 ± 1.4 | Control | 19 | -0.07 | -3.0 | 2.8 |

-1.222 ± 2.5 | Control | 19 | -1.25 | -6.2 | 3.9 |

-0.831 ± 1.8 | Control | 19 | -0.84 | -4.2 | 2.7 |

At this point we can use these MAP priors in the analysis of the pivotal trial as outlined in the following section.

## 6.6 Analysis of a new study

Assume now that the trial reads out and you have to analyze the study using the MAP prior that were derived. For this exercise will use a simulated data set, see code below for reference.

## Show the code

```
# Simulate summary data for 1 trial
#
# @param n1 The sample size in the 3 strata of the treatment arm. A vector of length 3
# @param n0 The sample size in the 3 strata of the control arm. A vector of length 3
# @param e1 The mean change from baseline at month 12 for the 3 strata of the treatment arm. A vector of length 3
# @param e0 The mean change from baseline at month 12 for the 3 strata of the control arm. A vector of length 3
# @param s1 The se of the response in the 3 strata of the treatment arm. A vector of length 3
# @param s0 The se of the response in the 3 strata of the control arm. A vector of length 3
# @param s1 The baseline score in the 3 strata of the treatment arm. A vector of length 3
# @param s0 The baseline score in the 3 strata of the control arm. A vector of length 3
#
# @return A data frame with simulated trial data
#
<- function(n1 = c(39, 18, 18),
simulate_1_trial_ipd n0 = c(26, 12, 12),
e1 = c( 5, 5, 5),
e0 = c( 0, 0, 0),
s1 = rep(5, 3),
s0 = rep(5, 3),
bs1 = c( 33,33,33),
bs0 = c( 33,33,33),
base_center=base_mean) {
# Arm 1: Treatment
# Arm 0: Control
# AGEGRPN 1: Age >=2 to <5
# AGEGRPN 2: Age >=5 to <13
# AGEGRPN 3: Age >=13 to <18
# Simulate individual patient data -------
<- sapply(1:3, function(i) rnorm(n = n1[i], mean = e1[i], sd = s1[i]))
d1 <- sapply(1:3, function(i) rnorm(n = n0[i], mean = e0[i], sd = s0[i]))
d0
# Simulate individual patient data for the baseline score -------
<- sapply(1:3, function(i) rnorm(n = n1[i], mean = bs1[i], sd = s1[i]))
b1 <- sapply(1:3, function(i) rnorm(n = n0[i], mean = bs0[i], sd = s0[i]))
b0
<- c("Age >=2 to <5", "Age >=5 to <13", "Age >=13 to <18")
agegr <- c(rep(agegr, n1), rep(agegr, n0))
AGEGROUP <- c(rep(1:3, n1), rep(1:3, n0))
strataN
<- c(unlist(b1), unlist(b0))
BASE <- c(unlist(d1), unlist(d0))
y
### Put generated values into a data frame -----------------------------------
<- data.frame(AGEGRPN = strataN,
out AGEGRP = AGEGROUP,
age_group = relevel(as.factor(strataN), "2"),
active = rep(c(1, 0), c(sum(n1), sum(n0))),
BASE = BASE,
cBASE= BASE - base_center,
CHG = y
)
out }
```

```
<- simulate_1_trial_ipd(
new_dat n1 = c(39, 18, 18),
n0 = c(26, 12, 12),
e1 = c( 2, 1, 0),
e0 = c(-4, -3, -2),
s1 = rep(5, 3),
s0 = rep(5, 3),
bs1 = c( 33,33,33),
bs0 = c( 33,33,33),
base_center=base_mean)
<- new_dat |>
new_dat mutate(Study = "New")
head(new_dat) |> kable()
```

AGEGRPN | AGEGRP | age_group | active | BASE | cBASE | CHG | Study |
---|---|---|---|---|---|---|---|

1 | Age >=2 to <5 | 1 | 1 | 27 | 5.8 | 2.1 | New |

1 | Age >=2 to <5 | 1 | 1 | 27 | 5.1 | 2.0 | New |

1 | Age >=2 to <5 | 1 | 1 | 34 | 12.4 | 8.5 | New |

1 | Age >=2 to <5 | 1 | 1 | 42 | 20.8 | 6.7 | New |

1 | Age >=2 to <5 | 1 | 1 | 30 | 9.0 | 8.4 | New |

1 | Age >=2 to <5 | 1 | 1 | 39 | 17.6 | -2.2 | New |

We first calculate simple summary statistics for the simulated data:

```
|>
new_dat summarise(BASE = mean(BASE),
cBASE = mean(cBASE),
CHGm = mean(CHG),
CHGsd = sd(CHG),
n = n(),
CHGse = CHGsd / sqrt(n),
.by=c(Study, AGEGRPN, AGEGRP, age_group, active)) |>
kable()
```

Study | AGEGRPN | AGEGRP | age_group | active | BASE | cBASE | CHGm | CHGsd | n | CHGse |
---|---|---|---|---|---|---|---|---|---|---|

New | 1 | Age >=2 to <5 | 1 | 1 | 34 | 12 | 2.76 | 5.7 | 39 | 0.91 |

New | 2 | Age >=5 to <13 | 2 | 1 | 34 | 13 | 2.22 | 4.2 | 18 | 1.00 |

New | 3 | Age >=13 to <18 | 3 | 1 | 35 | 13 | 0.91 | 6.3 | 18 | 1.50 |

New | 1 | Age >=2 to <5 | 1 | 0 | 33 | 12 | -3.43 | 4.5 | 26 | 0.89 |

New | 2 | Age >=5 to <13 | 2 | 0 | 34 | 12 | -3.04 | 5.2 | 12 | 1.51 |

New | 3 | Age >=13 to <18 | 3 | 0 | 32 | 11 | -0.61 | 4.4 | 12 | 1.26 |

We must define the model to be used. Note that now, since we have individual patient level data, we use only `CHG`

as response variable. Hence, the sampling standard deviation is estimated from the data (requiring us to set a prior for this parameter as well). It is obviously important that we must use the very same model specification as before, since otherwise we could not use the derived MAP prior easily:

```
<- bf(CHG ~ 1 + age_group + cBASE + active,
joint_model_ipd family=gaussian, center=FALSE)
```

And then use `brm`

to put it altogether. To use the previously derived priors, we need to use two options in the brm funciton. The `prior`

option indicates that a multivariate normal prior will be used, while the `stanvars`

will pass the actual prior distribution, which is a mixture distribution.

```
# Fit with historical data priors -------------------------------------------
<- 573391
mcmc_seed # use the MAP prior for the coefficients and encode with the gamma
# prior on the sigma that we anticipate an sd of 5
<- prior(mixmvnorm(prior_mix_w, prior_mix_m, prior_mix_sigma_L), class=b) +
map_prior prior(gamma(5,1), class=sigma)
<- brm(formula = joint_model_ipd,
bfit_stage2_ipd prior = map_prior,
stanvars = mixstanvar(prior_mix = post_coefs_mix),
data = new_dat,
seed = mcmc_seed,
refresh = 0,
backend = "cmdstanr")
```

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

` bfit_stage2_ipd`

```
Family: gaussian
Links: mu = identity; sigma = identity
Formula: CHG ~ 1 + age_group + cBASE + active
Data: new_dat (Number of observations: 125)
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.16 1.36 -5.04 0.16 1.00 1693 2140
age_group1 0.37 0.83 -1.26 2.01 1.00 2174 2478
age_group3 0.00 0.98 -1.94 1.91 1.00 2581 2772
cBASE -0.00 0.07 -0.14 0.14 1.00 1967 2630
active 4.09 0.92 2.35 5.88 1.00 3438 3067
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 5.31 0.35 4.67 6.02 1.00 3551 2696
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).
```

Having fit the model, the last step is to use it to test the statistical hypothesis:

`hypothesis(bfit_stage2_ipd, "active > 0", alpha = 0.05)`

```
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (active) > 0 4.1 0.92 2.6 5.6 3999 1 *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
```

## 6.7 Conclusion

This case study demonstrated how a random-effects meta-analysis model has been implemented with `brms`

. The model was implemented in a two-stage approach where first the model is fit with historical data to obtain meta-analytic-predictive priors. Later, when data from the study of interest is available, the model is fit using the study data and MAP priors.

The model enables borrowing between historical studies and age groups and hence a more informative MAP prior.

## 6.8 Exercises

Fit a model where only weakly informative priors are considered to analyze the data from the new study. This will help providing insights of the impact of using the historical data.

To add robustification to the MAP prior, non-informative components can added to the mixture distributions to protect against prior-data conflicts (Schmidli et al 2014). Fit a model where the MAP priors are robustified where the weight for the

non-informative component is 20%. The non-informative component is normally distributed with the mean 0 and variance \(\sigma^2\).

## 6.9 Solution to exercises

```
<- function(h1side, h2side) {
get2sidedCI = h2side$hypothesis |>
xx select(Hypothesis, CI.Lower, CI.Upper) |>
rename(CI.95.Lower=CI.Lower,
CI.95.Upper=CI.Upper)
data.frame(h1side$hypothesis, xx)
}= function(brm_res) {
get1sidehyp2sidedCI = hypothesis(brm_res, "active > 0", alpha=0.025)
h1side = hypothesis(brm_res, "active > 0", alpha=0.05)
h2side
= h2side$hypothesis |>
xx select(Hypothesis, CI.Lower, CI.Upper) |>
rename(CI.95.Lower=CI.Lower,
CI.95.Upper=CI.Upper)
data.frame(h1side$hypothesis, xx)
}
```

## Show the code

```
## Fit Non-Informative Prior using default priors -----------------------
head(new_dat)
<- brm(joint_model_ipd,
bfit_trial_ipd_noninf_flat seed = mcmc_seed,
refresh = 0,
data = new_dat)
## Fit Non-Informative Prior using wide Normal priors ------------------------
<- mixmvnorm(flat = c(1, rep(0, 5), diag(sd^2, 5)))
trial_prior_mix <- bf(CHG ~ 1 + age_group + BASE + active,
joint_model_ipd family=gaussian, center=FALSE)
<- brm(joint_model_ipd,
bfit_trial_ipd_noninf_norm prior = mix_prior,
stanvars = mixstanvar(prior_mix = trial_prior_mix),
seed = mcmc_seed, refresh = 0,
backend = "cmdstanr",
data = new_dat)
## Fit with historical data MAP Priors -------------------------------------------
<- brm(joint_model_ipd,
bfit_stage2_ipd prior = mix_prior,
stanvars = mixstanvar(prior_mix = post_coefs_mix),
seed = mcmc_seed, refresh = 0, backend = "cmdstanr")
= new_dat)
data
## Fit with historical data MAP Priors with robustification ----------------------
= seq(0, 1, 0.2)
ps = hn_ipd = list()
hp_ipd for (i in 1:length(ps)){
= ps[i]
p <- mixcombine(flat = trial_prior_mix,
rob_post_coefs_mix MAP = post_coefs_mix,
weight=c(p, 1 - p))
<- brm(joint_model_ipd,
bfit_stage2_rob_ipd prior = mix_prior,
stanvars = mixstanvar(prior_mix = rob_post_coefs_mix),
seed = mcmc_seed, refresh = 0, backend = "cmdstanr")
= new_dat)
data = get1sidehyp2sidedCI(bfit_stage2_rob_ipd)
hp_ipd[[i]] = glue::glue("MAP prior {p*100}% robust")
hn_ipd[[i]]
}
= get1sidehyp2sidedCI(bfit_trial_ipd_noninf_flat)
h1_ipd_noninf_flat = get1sidehyp2sidedCI(bfit_trial_ipd_noninf_norm)
h1_ipd_noninf_norm = get1sidehyp2sidedCI(bfit_stage2_ipd)
h1_ipd_map
= unlist(hn_ipd)
h1n_ipd = do.call(rbind, lapply(hp_ipd, function(x) x))
h1_ipd_rob
= data.frame(type = c("Non-informative prior. brm default",
res "Non-informative prior. wide normal prior",
"MAP prior",
h1n_ipd),bind_rows(
h1_ipd_noninf_flat,
h1_ipd_noninf_norm,
h1_ipd_map,
h1_ipd_rob
)) res
```