```
library(dplyr)
library(brms)
library(emmeans)
library(tidyr)
library(readr)
library(purrr)
library(ggplot2)
library(patchwork)
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(8904568)
<- list(adapt_delta = 0.95) control_args
```

# 12 Longitudinal data

## 12.1 Background

Many clinical trials assess efficacy and other endpoints at numerous timepoints post-baseline. For example, consider the case where a continuous response endpoint is assessed every 2 weeks post baseline (Week 2, Week 4, …, Week 12). Although clinical interest may focus on the response rate at a particular visit (e.g. at Week 12), efficacy data at other visits is of course still of interest, either in its own right, to build understanding of the full time vs. efficacy profile, or as a means to increase statistical power for estimation at Week 12.

There are multiple approaches for proceeding with estimation of the Week 12 treatment effect. Adopting a *cross-sectional approach* , one could ignoring all other post-baseline assessments except Week 12. If certain patients missed the Week 12 assessment, this missing data would need to be handled using an appropriate pre-defined strategy (e.g. for response rate endpoints, imputing missing outcomes with nonresponder status, carrying the last observation forward, multiple imputation approaches, or even dropping missing data altogether). Such approaches may lose signifcant power or incur bias if there is substantial missing data.

Alternatively *longitudinal models* incorporate the full post-baseline set of assessments to facilitate modeling of efficacy at all visits. In some cases, these may substantially increase statistical power.

The `brms`

package contains a deep suite of modelling tools for longitudinal data, including many options for modelling the mean trajectory across visits/time, and for modelling autocorrelated errors within patients across time. In this section we illustrate some potential uses of these tools in a clinical trial setting.

```
<- readr::read_csv(here::here("data", "longitudinal.csv"), show_col_types = FALSE) %>%
adpasi ::filter(TRT01P %in% c("PBO", "TRT")) %>%
dplyrmutate(AVISIT = factor(AVISIT, paste("Week", c(1, 2 * (1:6)))),
TRT01P = factor(TRT01P, c("PBO", "TRT")))
<- filter(adpasi, PARAMCD == "PASITSCO") pasi_data
```

## 12.2 Data

The example involves simulated results of a hypothetical Phase-II study of an experimental treatment for Psoriasis. Synthetic data are generated using the `mmrm`

package. We consider a subset of the study involving 100 patients, 50 of whom were randomized to receive placebo, and 50 of whom received treatment.

Efficacy was assessed using the Psoriasis Area and Severity Index (PASI), a numerical score which measures the severity and extent of psoriasis. This is assessed at baseline, and again at 7 post-baseline timepoints.

Two endpoints of interest based on the PASI score are (1) PASI change from baseline, and (2) the binary endpoint PASI 75, which defines a responder as any patient with at least a 75% change from baseline in PASI.

The data has been transformed to follow a typical CDISC Analysis Data Model (ADaM) format.

## 12.3 Models

There are a few key ingredients to a longitudinal model for the PASI score outcomes.

A

*mean model*which describes the expected value of the response over time, across treatments, and across values of any other relevant covariates.A

*correlation model*which describes the correlation structure of the error terms.

`brms`

offers many modelling options for each component.

### 12.3.1 Mean models

#### 12.3.1.1 Cell-means model

The most common approach for modeling the mean in clinical trial practice is to adopt a specification which allows the mean to vary freely across visits, without any parametric specification of the trajectory of the mean over time. Here we will call this the “cell-means” model, where levels of the treatment group and visit comprise the cells.

In this approach, we include all treatment-by-visit interactions. In `brms`

, the formula specification would be:

`CHG ~ BASE + TRT01P * AVISIT`

#### 12.3.1.2 Linear-in-time model

A far stronger assumption would be to assume the mean response is linear in time (i.e. the number of weeks post baseline for visit). In this case, the formula specification would be

`CHG ~ BASE + TRT01P * AVISITN`

Note the key difference is the use of `AVISITN`

(a numeric variable indicating the number of weeks post baseline) rather than `AVISIT`

(a factor variable for the visit id).

Such a model should be considered only for exploratory modelling purposes, and not for confirmatory analyses, the main reason being that the linearity assumption cannot be assessed at the trial design stage before having collected the data.

Even if the linearity assumption appeared reasonable over the time range explored in the trial, it should not be used to extrapolate outside the observed time period.

#### 12.3.1.3 Quadratic-in-time model

In the previous section, we saw there was a hint of curvature in the sample mean response over time:

Hence one might consider a mean model that assumes the mean response is quadratic in time rather than linear:

`CHG ~ BASE + TRT01P * AVISITN + TRT01P * AVISITN^2`

This model should similarly not be used to extrapolate outside the week 1-12 window.

#### 12.3.1.4 Gaussian process prior

An interesting nonparametric alternative to the just-discussed specifications is based on a Gaussian process prior for the mean across visits within treatments. As in the case of the cell-means model, this model does not assume any parametric shape for the mean response over time. It assumes only that the mean response over time follows a continuous curve which is assigned a Gaussian process (GP) prior.

While the GP itself is a continuous-time process, the joint distribution of a realization at any discrete set of timepoints is multivariate normal.

In `brms`

, the formula specification would be as follows:

`CHG ~ BASE + TRT01P + gp(AVISITN, by = TRT01P)`

The implication of this model is that there is a correlation between the mean response at any collection of visits, and the correlation between any pair of visits decays as the time between them increases, according to an exponential covariance function. See `?gp`

for additional details on the `brms`

implementation.

### 12.3.2 Correlation models

The repeated measurement of an endpoint over time on the same patient induces *autocorrelation* between the within-patient measurements. To illustrate, consider if we ignored the correlations and fit a model with uncorrelated errors. Below are scatterplots and correlation coefficients for the residuals of such a fitted model. The correlations are quite strong, especially between visits that are closer in time.

```
<- filter(adpasi, PARAMCD == "PASITSCO")
pasi_data <- lm(CHG ~ BASE + TRT01P * AVISIT, data = pasi_data)
lm_fit <- residuals(lm_fit)
res %>%
pasi_data mutate(res = res) %>%
select(SUBJID, AVISIT, res) %>%
pivot_wider(id_cols = SUBJID, names_from = AVISIT, values_from = res) %>%
select(-SUBJID) %>%
::ggpairs() GGally
```

Failing to model these within-subject correlations will result in a loss of statistical power (wider confidence intervals, less powerful tests). `brms`

offers several options for modelling within-subject correlations.

#### 12.3.2.1 Subject-level random effects

One way to achieve within-subject correlation is to use a model that includes a subject-level random effect. When the resulting mean response function is averaged over the distribution of the random effect terms (“integrating them out”), a uniform correlation is induced between all measurements on the same patient. The magnitude of the correlation is determined by the relative size of the random-effect variance and the error variance. (See exercise 1.)

In `brms`

, a subject-level random intercept can be easily added in the formula specification. For example:

`CHG ~ BASE + TRT01P * AVISIT + (1 | SUBJID)`

#### 12.3.2.2 Autoregressive correlation structures

A very common model for serial correlation is based on an *autoregressive process*. Under such a process, the error term \(\varepsilon_t\) at a timepoint \(t\) is explicitly dependent on some collection of preceeding error terms. Under a first-order autoregressive process, for example, it is dependent only on one of its predecessors: \[ \varepsilon_t = \alpha\varepsilon_{t-1} + Z_t,\] for \(t\geq 2\), where \(Z_t\) are iid \(\mathrm N(0, \sigma^2)\). The process is initialized with \(\varepsilon_1 \sim \mathrm N(0, \sigma^2 / (1 - \alpha^2))\).

The resulting covariance matrix of a collection \((\varepsilon_1,\ldots,\varepsilon_T)\) has a simple form; the reader is referred to the SAS paper “Guidelines for Selecting the Covariance Structure in Mixed Model Analysis” by Chuck Kincaid for detail on the AR structure (c.f. page 2), and other covariance structures.

The `autocor`

argument of `brm()`

and `brmsformula()`

is used to set an autoregressive correlation strucutre. For order-1 autoregressive. Below is a choice of AR(1) autocorrelation that we might consider in the psoriasis example:

`~ar(time = AVISIT, gr = SUBJID, p = 1)`

This choice implies that autocorrelation exists across visits within subjects.

#### 12.3.2.3 Compound symmetry

Another choice offered by `brms`

is that of compound symmetry. The reader is again referred to the SAS paper linked above and `?cosy`

for more information.

`~cosy(time = AVISIT, gr = SUBJID)`

#### 12.3.2.4 Other choices

`brms`

offers several other choices for autocorrelation models. See `?'autocor-terms'`

for more a listing.

We note that currently, *unstructured* correlation models (a standard choice for MMRM specifications in clinical trial protocols) re not supported by `brms`

.

However, fixed correlation structures are, and one could consider plugging in for this an unstructured correlation estimate from a frequentist MMRM. For example,

```
# fit a MMRM using gls
<- nlme::gls(CHG ~ BASE + TRT01P * AVISIT,
gls_fit data = pasi_data,
correlation = nlme::corSymm(form = ~ 1 | SUBJID))
# estimated correlation matrices by subject
<- nlme::corMatrix(gls_fit$modelStruct[[1]])
Sig_subj 1] Sig_subj[
```

```
$`1`
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.00000000 0.896714971 0.5883171 0.2669088 -0.03659905 -0.05666876
[2,] 0.89671497 1.000000000 0.6313654 0.3841036 0.03162521 0.03550198
[3,] 0.58831712 0.631365431 1.0000000 0.4509002 0.37604603 0.21075194
[4,] 0.26690877 0.384103637 0.4509002 1.0000000 0.63152510 0.47401182
[5,] -0.03659905 0.031625213 0.3760460 0.6315251 1.00000000 0.77205541
[6,] -0.05666876 0.035501982 0.2107519 0.4740118 0.77205541 1.00000000
[7,] -0.08775131 0.008196384 0.2465146 0.5921762 0.79671624 0.74110361
[,7]
[1,] -0.087751306
[2,] 0.008196384
[3,] 0.246514616
[4,] 0.592176181
[5,] 0.796716238
[6,] 0.741103610
[7,] 1.000000000
```

A block-diagonal matrix containing these estimated within-subject correlations could be plugged in as the `M`

argument of the `fcor()`

function in `brms`

.

### 12.3.3 Cross-sectional approaches

For modeling a continuous endpoint such as PASI, the cross-sectional analogue of the longitudinal models we’ve discussed is the *ANCOVA* model, in which a linear regression model (conditioning as before on the baseline PASI covariate) is fit to the Week-12 cross section of data.

```
# analysis data includes only the Week 12 cross section
<- filter(adpasi, PARAMCD == "PASITSCO", DROPFL, AVISITN == 12)
ancova_data
# formula specification for ANCOVA
<- bf(
ancova_formula ~ BASE + TRT01P,
CHG family = gaussian(),
center = FALSE,
nl = FALSE
)
<- get_prior(
ancova_prior
ancova_formula,data = ancova_data
)
<- brm(
ancova_fit
ancova_formula,prior = ancova_prior,
data = ancova_data,
seed = 46474576,
control = control_args,
refresh = 0
)
```

In the next section, this model is also explored for the purposes of comparing to the longitudinal approaches

## 12.4 Results

A typical estimand for longitudinal or ANCOVA involve Least Squares means (LS means) or estimated marginal means (EMM). Inference for EMMs is extremely convenient with `brms`

due to its integration with the `emmeans`

R package. The EMM can be roughly understood as the average, stratified by treatment group and visit, of the mean prediction across subjects in the trial.

We begin by briefly illustrating how `emmeans`

can be used with a `brmsfit`

object to estimate LS means and their contrasts.

```
<- filter(adpasi, PARAMCD == "PASITSCO", DROPFL)
analysis_data
<- bf(
brms_formula ~ BASE + TRT01P * AVISIT,
CHG autocor = ~ cosy(time = AVISIT, gr = SUBJID),
family = gaussian(),
center = FALSE,
nl = FALSE
)
<- get_prior(
prior
brms_formula,data = analysis_data
)
<- brm(
fit
brms_formula,prior = prior,
data = analysis_data,
seed = 46474576,
control = control_args,
refresh = 0
)
```

```
Running MCMC with 4 sequential chains...
Chain 1 finished in 12.1 seconds.
Chain 2 finished in 11.9 seconds.
Chain 3 finished in 12.5 seconds.
Chain 4 finished in 11.4 seconds.
All 4 chains finished successfully.
Mean chain execution time: 12.0 seconds.
Total execution time: 48.3 seconds.
```

` fit`

```
Family: gaussian
Links: mu = identity; sigma = identity
Formula: CHG ~ BASE + TRT01P * AVISIT
autocor ~ cosy(time = AVISIT, gr = SUBJID)
Data: analysis_data (Number of observations: 684)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Correlation Structures:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
cosy 0.41 0.04 0.33 0.51 1.00 3163 2356
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 12.95 1.87 9.29 16.56 1.00 2673
BASE -0.70 0.08 -0.85 -0.54 1.00 3934
TRT01PTRT -3.56 1.59 -6.62 -0.45 1.00 1884
AVISITWeek2 -0.16 1.16 -2.47 2.11 1.00 2034
AVISITWeek4 -2.49 1.23 -4.87 0.02 1.00 2158
AVISITWeek6 -0.87 1.31 -3.46 1.70 1.00 2142
AVISITWeek8 1.10 1.34 -1.55 3.77 1.00 2250
AVISITWeek10 -1.22 1.32 -3.78 1.43 1.00 2184
AVISITWeek12 -3.66 1.17 -6.01 -1.41 1.00 2064
TRT01PTRT:AVISITWeek2 -4.13 1.75 -7.59 -0.71 1.00 1942
TRT01PTRT:AVISITWeek4 -6.15 1.79 -9.73 -2.59 1.00 2173
TRT01PTRT:AVISITWeek6 -9.30 1.83 -12.89 -5.67 1.00 1916
TRT01PTRT:AVISITWeek8 -12.69 1.85 -16.38 -9.16 1.00 1981
TRT01PTRT:AVISITWeek10 -10.77 1.89 -14.41 -7.06 1.00 2080
TRT01PTRT:AVISITWeek12 -7.38 1.73 -10.82 -4.09 1.00 1890
Tail_ESS
Intercept 2488
BASE 3210
TRT01PTRT 2654
AVISITWeek2 2992
AVISITWeek4 2789
AVISITWeek6 2597
AVISITWeek8 2734
AVISITWeek10 2506
AVISITWeek12 2704
TRT01PTRT:AVISITWeek2 2370
TRT01PTRT:AVISITWeek4 2821
TRT01PTRT:AVISITWeek6 2290
TRT01PTRT:AVISITWeek8 2465
TRT01PTRT:AVISITWeek10 2343
TRT01PTRT:AVISITWeek12 2542
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 8.21 0.32 7.65 8.89 1.00 3278 2296
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).
```

```
# estimate the EMMs
<- emmeans(fit, c("TRT01P", "AVISIT"), nesting = list())
emm emm
```

```
TRT01P AVISIT emmean lower.HPD upper.HPD
PBO Week 1 -1.214 -3.38 0.8924
TRT Week 1 -4.782 -6.98 -2.3953
PBO Week 2 -1.378 -3.45 0.7517
TRT Week 2 -9.065 -11.30 -6.7377
PBO Week 4 -3.718 -6.02 -1.4343
TRT Week 4 -13.426 -15.69 -11.3379
PBO Week 6 -2.099 -4.40 0.3942
TRT Week 6 -14.976 -17.23 -12.6760
PBO Week 8 -0.144 -2.63 2.2711
TRT Week 8 -16.363 -18.77 -14.1551
PBO Week 10 -2.449 -4.94 0.0976
TRT Week 10 -16.780 -19.11 -14.3734
PBO Week 12 -4.874 -7.17 -2.7912
TRT Week 12 -15.795 -18.10 -13.5977
Point estimate displayed: median
HPD interval probability: 0.95
```

```
# estimate EMM contrasts
<- contrast(emm, method = "revpairwise", by = "AVISIT")
emm_contrasts emm_contrasts
```

```
AVISIT = Week 1:
contrast estimate lower.HPD upper.HPD
TRT - PBO -3.58 -6.53 -0.391
AVISIT = Week 2:
contrast estimate lower.HPD upper.HPD
TRT - PBO -7.67 -10.81 -4.698
AVISIT = Week 4:
contrast estimate lower.HPD upper.HPD
TRT - PBO -9.70 -12.87 -6.692
AVISIT = Week 6:
contrast estimate lower.HPD upper.HPD
TRT - PBO -12.88 -16.12 -9.505
AVISIT = Week 8:
contrast estimate lower.HPD upper.HPD
TRT - PBO -16.24 -19.47 -12.831
AVISIT = Week 10:
contrast estimate lower.HPD upper.HPD
TRT - PBO -14.33 -17.57 -10.757
AVISIT = Week 12:
contrast estimate lower.HPD upper.HPD
TRT - PBO -10.93 -14.22 -8.029
Point estimate displayed: median
HPD interval probability: 0.95
```

### 12.4.1 Fitting several models

In order to understand the impact of possible choices for the mean model and correlation model, we fit a series of models in loop.

Code to prepare a list of models to be fit:

## Code

```
<- function(adpasi){
setup_analyses
<- tribble(
formulas ~endpoint, ~formula_name, ~paramcd, ~family, ~formula, ~correl, ~longitudinal,
"PASI", "cell-means", "PASITSCO", gaussian(), CHG ~ BASE + TRT01P * AVISIT, TRUE, TRUE,
"PASI", "linear", "PASITSCO", gaussian(), CHG ~ BASE + TRT01P * AVISITN, TRUE, TRUE,
"PASI", "quadratic", "PASITSCO", gaussian(), CHG ~ BASE + TRT01P * AVISITN + TRT01P * AVISITN ^ 2, TRUE, TRUE,
"PASI", "gp", "PASITSCO", gaussian(), CHG ~ BASE + TRT01P + gp(AVISITN, by = TRT01P), TRUE, TRUE,
"PASI", "raneff", "PASITSCO", gaussian(), CHG ~ BASE + TRT01P * AVISIT + (1 | SUBJID), FALSE, TRUE,
"PASI", "cross-section", "PASITSCO", gaussian(), CHG ~ BASE + TRT01P, FALSE, FALSE
)
<- tribble(
datasets ~paramcd, ~longitudinal, ~missing_approach, ~analysis_data,
"PASITSCO", TRUE, "drop", filter(adpasi, PARAMCD == "PASITSCO", DROPFL),
"PASITSCO", FALSE, "drop", filter(adpasi, PARAMCD == "PASITSCO", DROPFL, AVISITN == 12)
)
<- tribble(
autocor ~correl, ~autocor_name, ~autocor,
TRUE, "AR1", ~ ar(time = AVISIT, gr = SUBJID, p = 1),
TRUE, "COSY", ~ cosy(time = AVISIT, gr = SUBJID)
)
<- formulas %>%
analyses full_join(datasets, c("paramcd", "longitudinal"), multiple = "all") %>%
full_join(autocor, c("correl"), multiple = "all") %>%
replace_na(list(autocor_name = "none"))
analyses
}
```

Code to fit the models in a loop using `clustermq`

:

## Code

```
::i_am("src/longitudinal/fit_models.R")
herelibrary(dplyr)
library(brms)
library(tidyr)
library(ggplot2)
library(readr)
library(here)
library(emmeans)
library(clustermq)
library(purrr)
source(here("src", "longitudinal", "setup_analyses.R"))
# iter <- Sys.getenv("BRMS_ITER")
# chains <- Sys.getenv("BRMS_CHAINS")
<- readr::read_csv(here("data", "longitudinal.csv")) %>%
adpasi mutate(AVISIT = factor(AVISIT, paste("Week", c(1, 2 * (1:6)))),
TRT01P = factor(TRT01P, c("PBO", "TRT")))
<- setup_analyses(adpasi)
analyses <- filter(analyses, endpoint == "PASI") %>%
analyses mutate(id = 1:n())
<- function(id, analysis_data, formula, formula_name,
fit_model
family, correl, autocor, autocor_name,save_individuals = FALSE){
# set up the formula ---------------------------------------------------------
if(correl){
<- autocor
autocor else{
} <- NULL
autocor
}
<- bf(
brms_formula
formula,autocor = autocor,
family = family,
center = FALSE,
nl = FALSE
)
# fit the model --------------------------------------------------------------
<- get_prior(
prior
brms_formula,data = analysis_data
)
cat("\n*** Fitting model", id, ": ", formula_name, ", ", autocor_name, "***\n")
<- brm(
fit
brms_formula,prior = prior,
cores = 4,
backend = "rstan",
data = analysis_data
)
if(save_individuals) saveRDS(fit, here::here("reports", paste0("longitudinal_fit_", id, ".rds")))
# estimate marginal means ----------------------------------------------------
cat("\n*** Estimating marginal means", id, ": ", formula_name, ", ", autocor_name, "***\n")
<- analysis_data %>%
base select(SUBJID, BASE) %>%
distinct()
<- analysis_data %>%
visits select(AVISIT, AVISITN) %>%
distinct()
<- analysis_data %>%
treatments select(TRT01P) %>%
distinct()
<- bind_rows(lapply(
emm split(visits, 1:nrow(visits)),
function(visit){
<- visit %>%
nd crossing(treatments) %>%
crossing(base) %>%
split(.$TRT01P)
<- map(nd, ~ posterior_linpred(fit, transform = TRUE, newdata = .))
lp
<- lapply(lp, rowMeans) %>%
marginal_mean as_tibble() %>%
mutate(diff = .[[2]] - .[[1]]) %>%
setNames(c(names(nd),
paste(rev(names(nd)), collapse = " - ")))
bind_cols(
visit,summarise_draws(as_draws_df(marginal_mean))
%>%
) mutate(
TRT01P = factor(variable, levels(analysis_data$TRT01P)),
contrast = case_when(
is.na(TRT01P) ~ variable,
TRUE ~ NA_character_
)
)
}
))
if(save_individuals) saveRDS(emm, here::here("reports", paste0("longitudinal_emm", id, ".rds")))
return(emm)
}
<- clustermq::Q_rows(
emm select(analyses, id, analysis_data, formula, formula_name, family, correl, autocor, autocor_name),
fit_model,const = list(save_individuals = FALSE),
n_jobs = nrow(analyses),
pkgs = c("brms", "emmeans", "tidyr", "dplyr", "purrr", "posterior"),
template = list(
walltime = 120,
job_name = "longitudinal_child",
log_file = here("reports", "longitudinal_child_%I.log"),
memory = 3000,
cores = 4
),job_size = 1
)
saveRDS(emm, file = here("reports", "longitudinal_fits.rds"))
saveRDS(analyses, file = here("reports", "longitudinal_analyses.rds"))
$emm <- emm analyses
```

### 12.4.2 PASI change from baseline: EMMs by visit

```
`geom_path()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
```

### 12.4.3 PASI change from baseline: EMM contrasts with placebo by visit

```
<- analyses %>%
contrast_results select(formula_name, autocor_name, emm) %>%
unnest(emm) %>%
::filter(!is.na(contrast)) %>%
dplyrmutate(formula_name = factor(formula_name, unique(formula_name)),
autocor_name = factor(autocor_name, unique(autocor_name)),
contrast = factor(contrast, unique(contrast))) %>%
rename(
mean_model = formula_name,
corr_model = autocor_name
)
<- mutate(slice(group_by(contrast_results, mean_model, corr_model), 1), estimate = 0)
blank
ggplot(
data = contrast_results,
mapping = aes(x = AVISIT, group = mean_model, color = mean_model,
y = median, ymin = q5, ymax = q95)
+
) geom_pointrange(position = position_dodge(0.5)) +
geom_path(position = position_dodge(0.5)) +
labs(x = "Visit", y = "PASI change from baseline\nTreatment - Control\nDifference in Estimated Marginal Means") +
geom_blank(data = blank) +
facet_grid(. ~ corr_model, labeller = label_both) +
scale_x_discrete(labels = levels(pasi_data$AVISIT)) +
scale_color_discrete("Mean model") +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1),
legend.position = "bottom")
```

### 12.4.4 Estimates of Week 12 efficacy

```
%>%
contrast_results ::filter(AVISITN == 12,
dplyr!mean_model %in% c("linear", "quadratic"),
%in% c("COSY", "none")) %>%
corr_model mutate(method = paste("mean model:", mean_model, "\ncorrelation model:", corr_model)) %>%
ggplot(aes(x = method, y = median, ymin = q5, ymax = q95)) +
geom_pointrange(position = position_dodge(0.6)) +
labs(x = "Modeling approach",
y = "PASI change from baseline\nTreatment - Control\nDifference in Estimated Marginal Means") +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1),
legend.position = "bottom")
```

### 12.4.5 Assessing model fit

In this section we illustrate the use of side-by-side posterior predictive checks to compare the fit of two models: one with a GP prior for the mean structure and one with a linear-in-time model for the mean, each with compound-symmetric autocorration models.

The checks are done across timepoint (# of weeks since baseline), by treatment group and quartile of the distribution of baseline PASI.

```
# Gaussian process model -------------------------------------
<- bf(
gp_formula ~ BASE + TRT01P + gp(AVISITN, by = TRT01P),
CHG autocor = ~ cosy(time = AVISIT, gr = SUBJID),
family = gaussian(),
center = FALSE,
nl = FALSE
)
<- get_prior(
gp_prior
gp_formula,data = analysis_data
)
<- brm(
gp_fit
gp_formula,prior = gp_prior,
data = analysis_data,
seed = 46474576,
control = control_args,
refresh = 0
)
```

`Start sampling`

```
Running MCMC with 4 sequential chains...
Chain 1 finished in 34.7 seconds.
Chain 2 finished in 38.5 seconds.
Chain 3 finished in 35.7 seconds.
Chain 4 finished in 38.4 seconds.
All 4 chains finished successfully.
Mean chain execution time: 36.8 seconds.
Total execution time: 147.7 seconds.
```

```
Warning: 10 of 4000 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.
```

```
# Linear mean model -------------------------------------
<- bf(
lm_formula ~ BASE + TRT01P * AVISITN,
CHG autocor = ~ cosy(time = AVISIT, gr = SUBJID),
family = gaussian(),
center = FALSE,
nl = FALSE
)
<- get_prior(
lm_prior
lm_formula,data = analysis_data
)
<- brm(
lm_fit
lm_formula,prior = lm_prior,
data = analysis_data,
seed = 46474576,
control = control_args,
refresh = 0
)
```

`Start sampling`

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

```
<- pasi_data %>%
base_quartiles select(SUBJID, BASE) %>%
distinct() %>%
pull(BASE) %>%
quantile(c(0, 0.25, 0.5, 0.75, 1))
<- (base_quartiles[1:4] + base_quartiles[2:5]) / 2
quartile_center
<- pasi_data %>%
full_newdata select(SUBJID, BASE) %>%
distinct() %>%
mutate(base_catn = as.numeric(cut(BASE, base_quartiles)),
base_cat = paste("Baseline PASI quartile", base_catn),
BASE = quartile_center[base_catn])
<- lapply(
pp_checks split(full_newdata, full_newdata$base_catn),
function(nd){
<- brms::pp_check(gp_fit, type = "ribbon_grouped", group = "TRT01P", x = "AVISITN",
p1 newdata = inner_join(nd, select(pasi_data, -BASE), "SUBJID", multiple = "all"),
y_draw = "points") +
labs(x = "Weeks post baseline",
y = "Week 12 PASI change from baseline",
title = paste("Quartile", unique(nd$base_catn), "of Baseline PASI")) +
theme(legend.position = "bottom") +
ylim(-50, 30)
<- brms::pp_check(lm_fit, type = "ribbon_grouped", group = "TRT01P", x = "AVISITN",
p2 newdata = inner_join(nd, select(pasi_data, -BASE), "SUBJID", multiple = "all"),
y_draw = "points") +
labs(x = "Weeks post baseline",
y = "Week 12 PASI change from baseline",
title = paste("Quartile", unique(nd$base_catn), "of Baseline PASI")) +
theme(legend.position = "bottom") +
ylim(-50, 30)
list(
gp = p1,
lm = p2
)
} )
```

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

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

Now we visualize the posterior predictive distribution against the observed data, first for the *Gaussian process* model:

## Code

`1]]$gp / pp_checks[[2]]$gp / pp_checks[[3]]$gp / pp_checks[[4]]$gp pp_checks[[`

and next for the *linear-in-time* model:

## Code

`1]]$lm / pp_checks[[2]]$lm / pp_checks[[3]]$lm / pp_checks[[4]]$lm pp_checks[[`

Another useful visualization for assessing model involves approximate leave-one-out (LOO) cross-validation techniques.

One can numerically compare the approximate expected log predictive density (ELPD) for holdout observations using `loo_compare`

.

```
<- loo(gp_fit)
loo_gp <- loo(lm_fit)
loo_lm loo_compare(loo_gp, loo_lm)
```

```
elpd_diff se_diff
gp_fit 0.0 0.0
lm_fit -20.8 5.2
```

This suggests the Gaussian-process based model has a superior model fit.

To visualize the predictive densities versus observed values, `pp_check()`

with type `"loo_pit"`

. In these visualizations, the observed outcomes are compared with their respective leave-one-out predictive distributions using the probability integral transformation (PIT). In an ideal model fit, the resulting PIT-transformed values would be uniformly distributed: i.e. the blue points and dashed lines would align with the distribution function of a Uniform(0,1) variable (solid line).

```
<- pp_check(gp_fit, type = "loo_pit") + geom_abline(intercept = 0, slope = 1) + ggtitle("LOO PIT: GP model")
gp_loo <- pp_check(lm_fit, type = "loo_pit") + geom_abline(intercept = 0, slope = 1) + ggtitle("LOO PIT: linear-in-time model")
lm_loo + lm_loo gp_loo
```

## 12.5 Exercises

- Using
`brms`

, fit a cross-sectional (non longitudinal) logistic regression model to the binary PASI 75 endpoint at Week 12. The only covariate term should be the treatment effect. Use the following analysis data (which uses only the Week 12 outcomes, and uses nonresponder imputation for any missing outcomes at Week 12). Use the`emmeans`

package and function to estimate the marginal mean Week 12 response rates by treatment. Use`emmeans::contrast`

to estimate the treatment-vs-control contrasts in response rates.

```
<- readr::read_csv(here("data", "longitudinal.csv"), show_col_types = FALSE) %>%
adpasi ::filter(TRT01P %in% c("PBO", "TRT")) %>%
dplyrmutate(AVISIT = factor(AVISIT, paste("Week", c(1, 2 * (1:6)))),
TRT01P = factor(TRT01P, c("PBO", "TRT")))
<- filter(adpasi, PARAMCD == "PSRS75", NRFL, AVISIT == "Week 12") analysis_data1
```

## Code

```
# solution
<- brm(AVAL ~ TRT01P, family = bernoulli(), data = analysis_data1,
fit1 silent = 2, refresh = 0)
<- emmeans(fit1, c("TRT01P"), transform = "response")
emm1
emm1contrast(emm1, method = "revpairwise")
```

- Fit a longitudinal model to the binary PASI 75 endpoint. Use a Gaussian process prior (stratified by treatment arm) for the mean across weeks and an AR(1) process to model autocorrelation in the residuals. Use the following analysis data (in which any missing assessments are not imputed). Use the
`emmeans`

package and function to estimate the marginal mean response rates by treatment and visit. Use`emmeans::contrast`

to estimate the treatment-vs-control contrasts in response rates by visit. How does the inference for Week 12 response rate difference compare to the cross-sectional model fit from 1?

`<- filter(adpasi, PARAMCD == "PSRS75", !MISSFL) analysis_data2 `

## Code

```
# solution
<- brm(
fit2 bf(
~ TRT01P + gp(AVISITN, by = TRT01P),
AVAL autocor = ~ cosy(time = AVISIT, gr = SUBJID)
),data = analysis_data2,
silent = 2,
refresh = 0
)<- emmeans(fit2, c("TRT01P", "AVISITN"), cov.keep = c("TRT01P", "AVISITN"), nesting = list("AVISTN" = "AVISIT"), transform = "response")
emm2
emm2contrast(emm2, method = "revpairwise", by = "AVISITN")
```

- Use Leave-One-Out cross validation to compare the model fits to the observed Week 12 data. (Hint: use
`loo(fit, newdata = newdata)`

with the below choice of newdata, then use`loo_compare`

). Which model has better predictive power?

`<- filter(adpasi, PARAMCD == "PSRS75", !MISSFL, AVISIT == "Week 12") newdata `

## Code

```
# solution
<- loo(fit1, newdata = newdata)
loo1 <- loo(fit2, newdata = newdata)
loo2 loo_compare(loo1, loo2)
```

- (Advanced) In this exercise, we will use a model fit to the continuous PASI change from baseline to do inference on the binary 75 response rate. First, fit a GP-based model to the PASI endpoint using the following code:

```
<- brm(
pasi_fit bf(
~ BASE + TRT01P + gp(AVISITN, by = TRT01P),
CHG autocor = ~ cosy(time = AVISIT, gr = SUBJID)
),data = pasi_data,
family = gaussian()
)
```

Next, create two sets of covariate values that includes “counterfactuals” for each patient as if they received both treatment and control. The following code is convenient:

```
<- mutate(filter(pasi_data, AVISIT == "Week 12"),
x_treatment TRT01P = factor("TRT", levels(pasi_data$TRT01P)))
<- mutate(filter(pasi_data, AVISIT == "Week 12"),
x_control TRT01P = factor("PBO", levels(pasi_data$TRT01P)))
```

Next, for each of `x_treatment`

and `x_control`

:

- use
`posterior_predict`

using the`newdata`

argument to sample from the posterior predictive distribution for PASI change from baseline at each level`x_treatment`

and`x_control`

, respectively. The result of`posterior_predict`

whose columns are posterior predictions of PASI change from baseline for individual patients. - Divide these predictions by the baseline PASI for the respective patients to obtain predictive draws for PASI % change from baseline.
- Convert the result to binary indicators of PASI % change from baseline being below \(-75\%\).
- For each MCMC iteration, compute the percentage of responders to obtain a vector of posterior draws for the marginal mean response rates for
`x_treatment`

and`x_control`

, respectively. - How do the median and credible intervals for the marginal mean response rates compare to the results of
`emmeans`

from exercises 1 and 2? - Graph the posterior density for the difference in marginal mean response rates for treatment minus control.

## Code

```
# solution
# predicted change from baseline
<- posterior_predict(pasi_fit, newdata = x_treatment)
pasi_chg_treatment <- posterior_predict(pasi_fit, newdata = x_control)
pasi_chg_control
# predicted % change from baseline
<- sweep(pasi_chg_treatment, 2, x_treatment$BASE, '/')
pasi_pchg_treatment <- sweep(pasi_chg_control, 2, x_control$BASE, '/')
pasi_pchg_control
# predicted PASI 75 response
<- pasi_pchg_treatment < -0.75
pasi_rr_treatment <- pasi_pchg_control < -0.75
pasi_rr_control
# marginal PASI 75 response rates
<- rowMeans(pasi_rr_treatment)
marginal_rr_treatment <- rowMeans(pasi_rr_control)
marginal_rr_control
# How do the median and credible intervals compare to emm1 and emm2?
apply(cbind(trt = marginal_rr_treatment,
ctrl = marginal_rr_control), 2, median)
::HPDinterval(coda::as.mcmc(cbind(trt = marginal_rr_treatment,
codactrl = marginal_rr_control)))
emm1
emm2
# visualize posterior distribution for marginal mean response rates by treatment
qplot(
x = cbind(marginal_rr_control, marginal_rr_treatment),
group = cbind(rep("control", length(marginal_rr_treatment)),
rep("treatment", length(marginal_rr_treatment))),
color = cbind(rep("control", length(marginal_rr_treatment)),
rep("treatment", length(marginal_rr_treatment))),
geom = "density"
+ scale_color_discrete(NULL) + theme(legend.position = "bottom") +
) labs(x = "Marginal mean response rate",
y = "Posterior density")
# visualize posterior distribution for difference in response rates
qplot(
x = marginal_rr_treatment - marginal_rr_control,
geom = "density"
+
) labs(x = "Difference in marginal mean response rate",
y = "Posterior density")
```