R Session

To run exercises on your hardware please follow the instructions below.

Stan

Please first install the cmdstanr R package following the instructions from their Getting started with CmdStanR article. This also covers installation of the necessary compiler and accompanying C++ toolchain as well as cmdstan itself.

An alternative is to use rstan as described on the wiki page for rstan to get started. In this case, please use as backend for brms the setting rstan. Note that the use of the cmdstanr backend is recommended over rstan allowing to use more recent versions of Stan. Moreover, the instructions for installing a C++ toolchain are better described as part of the respective cmdstan guide chapter.

C++ toolchain installation

• Windows: cmdstanr can automatically install the RTools based C++ toolchain via

cmdstanr::check_cmdstan_toolchain(fix=T)

• macOS: The project macrtools provides a comprehensive and simple installation of a full toolchain as used by CRAN https://mac.thecoatlessprofessor.com/macrtools/

# install.packages("remotes")
remotes::install_github("coatless-mac/macrtools")
macrtools::macos_rtools_install()

R packages

Please also ensure to have all R packages installed as listed on the next preliminary code section.

# in case packages are missing, please run:
install.packages(c("here", "ggplot2", "dplyr", "knitr", "brms", "posterior", "tidybayes", "bayesplot", "RBesT", "ggrepel", "patchwork", "ggdist", "ggthemes", "withr", "simsurv", "gt"))

Web-site

To download the full web-site to complement the material you may do so using the following R code snippet:

bamdd_zip <- tempfile(fileext=".zip")
download.file("https://github.com/Novartis/bamdd/archive/refs/heads/main.zip", bamdd_zip)
## extracts web-site into the user's home
unzip(bamdd_zip, exdir=normalizePath("~"))
browseURL(normalizePath(file.path("~", "bamdd-main")))
# to install all dependencies needed to build the web-site, please run
source(file.path("~", "bamdd-main", "src", "install_dependencies.R"))

Preliminary code

First load the required packages and set some options as described in the case study.

# uncomment below and set to the filename of your R file for the
# exercises
# here::i_am("your-exercise-file.R")
library(here)
library(ggplot2)
library(dplyr)
library(knitr)
library(brms)
library(posterior)
library(tidybayes)
library(bayesplot)
library(RBesT)
library(ggrepel)
library(patchwork)
library(ggdist)

options(
  # how many processor cores would you like to use?
  mc.cores = 4, 
  # how would you like to access Stan?
  brms.backend = "cmdstanr",
  # cache model binaries
  cmdstanr_write_stan_file_dir=here::here("_brms-cache"),
  # no need to normalize likelihoods
  brms.normalize = FALSE
)

dir.create(here::here("_brms-cache"), FALSE)
set.seed(254335)

# restrict output precision
options(digits=3)

# allow wide printing
options(width=120)

# in case rstan is used as backend, consider to enable model caching
# by enabling the following line
# rstan_options(auto_write = TRUE)

# Set defaults for ggplot2 ----
theme_set(theme_bw(base_size=18))

scale_colour_discrete <- function(...) {
  # Alternative: ggsci::scale_color_nejm(...)
  # scale_colour_brewer(..., palette="Set1")
  ggthemes::scale_colour_colorblind(...)
}
scale_fill_discrete <- function(...) {
  # Alternative: ggsci::scale_fill_nejm(...)
  #scale_fill_brewer(..., palette="Set1")
  ggthemes::scale_fill_colorblind(...)
}
scale_colour_continuous <- function(...) {
  scale_colour_viridis_c(..., option="turbo")
}
update_geom_defaults("point", list(size=2))
update_geom_defaults("line", list(linewidth=1.5))

Resources

• brms
• cmdstanr
• bayesplot
• loo
• posterior
• tidybayes
• mc-stan documentation overview
• Stan user forum

Exercise setup

Create the region dataset used in the example:

withr::with_seed(654873,
                 AS_region <- bind_cols(RBesT::AS, region=sample(c("asia", "europe", "north_america"), 8, TRUE)))
gt::gt(AS_region)

And fit the model as shown in the case study:

model <- bf(r | trials(n) ~ 1 + (1 | study), family=binomial)

model_prior <- prior(normal(0, 2), class=Intercept) +
    prior(normal(0, 1), class=sd, coef=Intercept, group=study)

map_mc_brms  <- brm(model, AS_region, prior = model_prior,
                    seed = 4767, control=list(adapt_delta=0.99),
                    silent = 2, refresh = 0)
map_mc_brms

 Family: binomial 
  Links: mu = logit 
Formula: r | trials(n) ~ 1 + (1 | study) 
   Data: AS_region (Number of observations: 8) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Multilevel Hyperparameters:
~study (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.38      0.20     0.06     0.88 1.00     1088     1254

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -1.10      0.19    -1.47    -0.70 1.00     1398     1575

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).

Exercise 1: Posterior predictive check

pp_count_AS_region <- posterior_predict(map_mc_brms)

pp_rate_AS_region <- sweep(pp_count_AS_region, 2, AS_region$n, "/")

head(pp_rate_AS_region) |> kable()

with(AS_region, ppc_intervals(r / n, pp_rate_AS_region))  +
    scale_x_continuous("Study", breaks=1:nrow(AS_region), labels=AS_region$study) +
    ylab("Responder Rate") +
    coord_flip() +
    theme(legend.position="right",
          # suppress vertical grid lines for better readability of intervals
          panel.grid.major.y = element_blank())

pp_AS_region <- AS_region |>
    add_predicted_rvars(map_mc_brms, value="pp_r") |>
    mutate(pp_rate=pp_r/n)

pp_AS_region

# A tibble: 8 × 6
  study       n     r region               pp_r       pp_rate
  <chr>   <dbl> <dbl> <chr>          <rvar[1d]>    <rvar[1d]>
1 Study 1   107    23 europe         24.2 ± 5.7  0.23 ± 0.053
2 Study 2    44    12 asia           11.5 ± 3.5  0.26 ± 0.080
3 Study 3    51    19 europe         16.1 ± 4.4  0.32 ± 0.086
4 Study 4    39     9 europe          9.5 ± 3.3  0.24 ± 0.085
5 Study 5   139    39 asia           37.8 ± 6.9  0.27 ± 0.050
6 Study 6    20     6 north_america   5.3 ± 2.3  0.27 ± 0.115
7 Study 7    78     9 europe         13.5 ± 4.7  0.17 ± 0.060
8 Study 8    35    10 asia            9.3 ± 3.2  0.27 ± 0.091

with(pp_AS_region, ppc_intervals(r / n, as_draws_matrix(pp_rate)))  +
    scale_x_continuous("Study", breaks=1:nrow(AS_region), labels=AS_region$study) +
    ylab("Responder Rate") +
    coord_flip() +
    theme(legend.position="right",
          # suppress vertical grid lines for better readability of intervals
          panel.grid.major.y = element_blank())

Exercise 2: Fixed vs random effect

AS_region_all <- data.frame(region=c("asia", "europe", "north_america", "other")) |>
    mutate(study=paste("new_study", region, sep="_"), r=0, n=6)

## to get brms to include the other factor level in the model, we have
## to add a fake row with region "other" and n=0
AS_region_2 <- mutate(bind_rows(AS_region, mutate(AS_region_all, n=0)[4,]),
                      region=factor(region, levels=c("asia", "europe", "north_america", "other")))

str(AS_region_2)

'data.frame':   9 obs. of  4 variables:
 $ study : chr  "Study 1" "Study 2" "Study 3" "Study 4" ...
 $ n     : num  107 44 51 39 139 20 78 35 0
 $ r     : num  23 12 19 9 39 6 9 10 0
 $ region: Factor w/ 4 levels "asia","europe",..: 2 1 2 2 1 3 2 1 4

model_fixed <- bf(r | trials(n) ~ 1 + region + (1|study), family=binomial)

get_prior(model_fixed, AS_region_2)

                prior     class                coef group resp dpar nlpar lb ub tag       source
               (flat)         b                                                          default
               (flat)         b        regioneurope                                 (vectorized)
               (flat)         b regionnorth_america                                 (vectorized)
               (flat)         b         regionother                                 (vectorized)
 student_t(3, 0, 2.5) Intercept                                                          default
 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)

model_fixed_prior <- prior(normal(0, 2), class=Intercept) +
    prior(normal(0, 1), class=sd, coef=Intercept, group=study) +
    prior(normal(0, log(2)/1.96), class=b)

fixed_mc_brms  <- brm(model_fixed, AS_region_2, prior=model_fixed_prior,
                      seed=4767,
                      silent = 2, refresh = 0, control=list(adapt_delta=0.99))

fixed_mc_brms

 Family: binomial 
  Links: mu = logit 
Formula: r | trials(n) ~ 1 + region + (1 | study) 
   Data: AS_region_2 (Number of observations: 9) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Multilevel Hyperparameters:
~study (Number of levels: 9) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.38      0.22     0.04     0.92 1.00     1314     1612

Regression Coefficients:
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept              -1.02      0.23    -1.51    -0.57 1.00     1778     1669
regioneurope           -0.17      0.26    -0.66     0.34 1.00     2226     2381
regionnorth_america     0.03      0.31    -0.60     0.63 1.00     3351     2898
regionother             0.00      0.35    -0.68     0.66 1.00     3691     2599

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).

post_regions <- posterior_linpred(fixed_mc_brms,
                                  newdata=AS_region_all,
                                  transform=TRUE,
                                  allow_new_levels=TRUE,
                                  sample_new_levels="gaussian")

head(post_regions) |> kable(digits=3)

colnames(post_regions) <- AS_region_all$region
map_region <- list()
for(r in AS_region_all$region) {
    map_region[[r]] <- mixfit(post_regions[,r], type="beta", Nc=3, constrain_gt1=TRUE)
}

kable(bind_rows(lapply(map_region, summary), .id="MAP"), digits=3)

sapply(map_region, ess)

         asia        europe north_america         other 
         32.6          33.6          18.4          16.3 

Exercise 3: Normal endpoint meta-analysis

crohn <- RBesT::crohn

crohn_sigma <- 88
crohn$y.se <- crohn_sigma/sqrt(crohn$n)

library(RBesT)
set.seed(1234)
rbest_normal_map_mcmc <- gMAP(cbind(y, y.se) ~ 1 | study, 
                              weights=n, data=crohn,
                              family=gaussian,
                              beta.prior=cbind(0, crohn_sigma),
                              tau.dist="HalfNormal",tau.prior=cbind(0,crohn_sigma/2))

model_normal <- bf(y | se(y.se) ~ 1 + (1 | study), family=gaussian())

prior_normal <- prior(normal(0, 88), class=Intercept) +
    prior(normal(0, 88/2), class=sd, coef=Intercept, group=study)

brms_normal_map_mcmc <- brm(model_normal, crohn, prior=prior_normal,
                            seed=4767,
                            silent = 2, refresh = 0, control=list(adapt_delta=0.99))

## comparing the outputs we see that the random effect posterior
## matches...
rbest_normal_map_mcmc

Generalized Meta Analytic Predictive Prior Analysis

Call:  gMAP(formula = cbind(y, y.se) ~ 1 | study, family = gaussian, 
    data = crohn, weights = n, tau.dist = "HalfNormal", tau.prior = cbind(0, 
        crohn_sigma/2), beta.prior = cbind(0, crohn_sigma))

Exchangeability tau strata: 1 
Prediction tau stratum    : 1 
Maximal Rhat              : 1 
Estimated reference scale : 88 

Between-trial heterogeneity of tau prediction stratum
       mean median   sd q2.5  q50 q97.5
tau[1] 14.4   12.7 9.48 1.37 12.7  38.1

MAP Prior MCMC sample
                 mean median   sd  q2.5   q50 q97.5
theta_resp_pred -50.1  -48.6 19.3 -92.8 -48.6 -10.9

brms_normal_map_mcmc

 Family: gaussian 
  Links: mu = identity 
Formula: y | se(y.se) ~ 1 + (1 | study) 
   Data: crohn (Number of observations: 6) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Multilevel Hyperparameters:
~study (Number of levels: 6) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    14.25     10.02     1.07    40.33 1.00      699      885

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   -49.26      9.10   -68.11   -32.85 1.00     1061     1037

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).

brms_normal_map <- posterior_epred(brms_normal_map_mcmc,
                                   newdata=data.frame(study="new", n=1, y=0, y.se=88),
                                   allow_new_levels=TRUE,
                                   sample_new_levels="gaussian")

## ... and the MAP prior is also the same
summarise_draws(brms_normal_map, mean, sd, ~quantile2(., probs = c(0.025, 0.5, 0.975)))

# A tibble: 1 × 6
  variable  mean    sd  q2.5   q50 q97.5
  <chr>    <dbl> <dbl> <dbl> <dbl> <dbl>
1 ...1     -48.9  18.8 -90.4 -48.0 -10.5

Peanut allergy exercise overview

These exercises accompany the dose finding case study.

In a hypothetical phase 2b trial children and adolescents with peanut allergy were randomly assigned to placebo or 5 doses (5 to 300 mg) of a new drug.

After 26 weeks the patients underwent a double-blind placebo-controlled food challenge and the primary endpoint was the proportion of patients that could ingest 600 mg or more of peanut protein without experiencing dose-limiting symptoms.

The plots below show an overview of the data. Note that no placebo patients were considered “responders” per the primary endpoint definition.

peanut <- tibble(
  TRT01P = structure(c(
    6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
    6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 5L, 5L, 5L, 5L, 5L,
    5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
    5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
    4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
    4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
    3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L
  ),levels = c("PBO", "5 mg", "15 mg", "50 mg", "150 mg", "300 mg"), class = "factor"),
  dose = c(
    300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L,
    300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L,
    300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 300L, 150L, 150L, 150L, 150L, 150L, 150L, 150L,
    150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L,
    150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L,
    150L, 150L, 150L, 150L, 150L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L,
    50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L,
    50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 50L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,
    15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,
    15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
    5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
    5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L),
  AVAL = c(
    1000L, 100L, 0L, 600L, 600L, 300L, 300L, 10L, 10L, 100L, 1000L, 300L, 10L, 1000L, 1000L, 1000L, 3L, 3L, 600L, 300L,
    1000L, 1000L, 100L, 1000L, 600L, 600L, 1000L, 1000L, 1000L, 300L, 600L, 1000L, 1000L, 1000L, 1000L, 1000L, 300L,
    1000L, 600L, 1000L, 300L, 1000L, 600L, 600L, 1000L, 1000L, 600L, 1000L, 1000L, 1000L, 1000L, 600L, 1000L, 1000L,
    1000L, 0L, 600L, 300L, 1000L, 1000L, 1000L, 100L, 1000L, 1000L, 600L, 1000L, 1000L, 1000L, 600L, 300L, 600L, 600L,
    100L, 600L, 1000L, 300L, 10L, 3L, 1000L, 300L, 300L, 1000L, 300L, 300L, 10L, 1000L, 1000L, 10L, 10L, 600L, 3L, 10L,
    600L, 600L, 600L, 1000L, 100L, 1000L, 1000L, 10L, 1000L, 3L, 10L, 1000L, 100L, 1000L, 100L, 10L, 300L, 10L, 100L,
    1000L, 0L, 1000L, 100L, 3L, 10L, 100L, 100L, 300L, 1000L, 1000L, 1000L, 10L, 1000L, 3L, 3L, 600L, 600L, 10L, 3L,
    600L, 600L, 300L, 300L, 1000L, 3L, 3L, 1000L, 10L, 1000L, 1000L, 600L, 100L, 300L, 600L, 10L, 100L, 0L, 100L, 3L,
    0L, 10L, 3L, 600L, 300L, 300L, 300L, 600L, 300L, 100L, 3L, 0L, 10L, 600L, 300L, 10L, 300L, 600L, 1000L, 600L, 600L,
    0L, 0L, 600L, 600L, 600L, 0L, 0L, 300L, 100L, 0L, 10L, 300L, 1000L, 300L, 600L, 600L, 300L, 10L, 600L, 100L, 100L,
    300L, 3L, 3L, 300L, 1000L, 10L, 3L, 100L, 3L, 100L, 100L, 300L, 3L, 3L, 600L, 300L, 3L, 3L, 3L, 300L, 3L, 0L, 10L,
    3L, 300L, 10L, 10L, 600L, 0L, 300L, 600L, 0L, 0L, 100L, 100L, 10L, 100L, 10L, 100L, 600L, 0L, 600L, 0L, 10L, 100L,
    0L, 0L, 0L, 0L, 3L, 10L, 3L, 300L, 600L, 0L, 0L, 300L, 10L, 10L, 100L, 300L, 0L, 0L, 0L, 3L, 0L, 0L, 10L, 0L, 300L,
    3L, 0L, 0L, 0L, 0L, 100L, 3L, 0L, 3L, 10L, 10L, 3L, 0L, 3L, 10L, 3L, 0L, 0L, 3L, 100L, 0L, 0L, 300L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 3L, 0L, 3L, 0L, 0L, 0L, 0L),
  PARAM = structure(
    c(
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L),
    levels = "Tolerated amount of peanut protein (mg)", class = "factor"),
  USUBJID = c(
    129L, 103L, 104L, 185L, 34L, 15L, 140L, 151L, 7L, 144L, 200L, 23L, 138L, 115L, 255L, 147L, 224L, 101L, 156L, 284L,
    136L, 248L, 179L, 298L, 168L, 295L, 289L, 241L, 36L, 27L, 123L, 266L, 11L, 291L, 236L, 130L, 173L, 195L, 203L, 160L,
    274L, 167L, 290L, 94L, 9L, 269L, 122L, 135L, 64L, 26L, 95L, 10L, 5L, 234L, 161L, 299L, 88L, 69L, 35L, 233L, 286L,
    85L, 91L, 189L, 80L, 152L, 223L, 287L, 244L, 57L, 108L, 18L, 62L, 157L, 300L, 283L, 164L, 243L, 89L, 220L, 24L, 271L,
    166L, 118L, 201L, 127L, 121L, 41L, 267L, 213L, 49L, 73L, 202L, 134L, 112L, 25L, 227L, 29L, 251L, 273L, 119L, 132L,
    74L, 270L, 83L, 37L, 181L, 258L, 253L, 48L, 120L, 54L, 277L, 176L, 65L, 264L, 107L, 171L, 262L, 162L, 187L, 272L,
    288L, 294L, 245L, 109L, 172L, 204L, 275L, 22L, 66L, 186L, 247L, 17L, 149L, 141L, 177L, 280L, 216L, 40L, 75L, 263L,
    246L, 14L, 81L, 260L, 153L, 45L, 237L, 8L, 117L, 86L, 296L, 146L, 154L, 116L, 38L, 100L, 191L, 175L, 92L, 158L, 192L,
    180L, 256L, 254L, 125L, 222L, 145L, 261L, 155L, 159L, 3L, 206L, 278L, 19L, 31L, 63L, 208L, 55L, 259L, 218L, 111L,
    226L, 33L, 2L, 44L, 297L, 53L, 87L, 16L, 28L, 90L, 207L, 56L, 137L, 128L, 178L, 142L, 143L, 148L, 193L, 229L, 265L,
    97L, 252L, 205L, 150L, 165L, 188L, 52L, 99L, 93L, 1L, 221L, 124L, 210L, 6L, 232L, 21L, 211L, 163L, 96L, 60L, 183L,
    190L, 242L, 42L, 46L, 67L, 126L, 209L, 72L, 194L, 238L, 184L, 39L, 105L, 249L, 61L, 113L, 30L, 77L, 12L, 4L, 51L,
    139L, 20L, 268L, 215L, 292L, 217L, 199L, 32L, 276L, 47L, 225L, 230L, 79L, 71L, 98L, 50L, 13L, 76L, 231L, 250L, 58L,
    68L, 239L, 198L, 293L, 212L, 110L, 59L, 182L, 133L, 170L, 43L, 282L, 281L, 131L, 114L, 196L, 214L, 285L, 70L, 102L,
    279L, 106L, 197L, 257L, 228L, 84L, 235L, 78L, 169L, 240L, 219L, 174L, 82L),
  CRIT1 = structure(c(
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L),
    levels = "Tolerating >=600 mg of peanut protein without dose-limiting symptoms", class = "factor"),
  CRIT1FL = structure(c(
    2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L,
    2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
    2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L,
    1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L,
    1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L,
    2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L,
    2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L,
    1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L),
    levels = c("N", "Y"), class = "factor"),
  CRIT2 = structure(c(
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L),
    levels = "Tolerating >=1000 mg of peanut protein without dose-limiting symptoms", class = "factor"),
  CRIT2FL = structure(c(
    2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L,
    2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L,
    1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L,
    1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L,
    1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L,
    2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L),
    levels = c("N", "Y"), class = "factor" ) )

peanut |>
  ggplot(aes(x=factor(AVAL), fill=TRT01P)) +
  geom_bar( position=position_dodge(preserve = "single")) +
  scale_fill_manual(values=c("royalblue", "#fb6a4a", "#ef3b2c", "#cb181d", "#a50f15", "#67000d")) +
  xlab("Peanut protein tolerated without dose-limiting symptoms (mg)") +
  ylab("Patients") +
  theme(legend.position="right")

peanut |>
  group_by(TRT01P) |>
  summarize(proportion = sum(CRIT1FL=="Y") / n()) |>
  ggplot(aes(x=TRT01P, y=proportion, fill=TRT01P)) +
  geom_bar(stat="identity", position=position_dodge()) +
  geom_text(aes(label=scales::percent(proportion)), size=7, vjust=c(0, rep(1.5, 5))) +
  scale_fill_manual(values=c("royalblue", "#fb6a4a", "#ef3b2c", "#cb181d", "#a50f15", "#67000d")) +
  scale_y_continuous(breaks=seq(0, 0.7, 0.1), labels=scales::percent) +
  xlab("Tolerating >=600 mg of peanut protein without dose-limiting symptoms") +
  ylab("Percentage of patients")

Exercise 1: Emax logistic regression

Fit a sigmoid Emax logistic regression model (i.e. family=binomial(link="logit")). To accelerate likelihood evaluations, first summarize data as “responders” y out of n patients per dose. We tell brms to use this information using y | trials(n) ~ ....

We expect the true placebo proportion to be around 0.05 (-2.944 on the logit scale) with much more than 0.1 or much less than 0.02 considered unlikely. It is a-priori at least possible that there is a huge treatment effect such as 95% versus 5% responders (difference on the logit scale close to 6), but we are at least mildly skeptical.

The dose response is a-priori somewhat likely to follow an Emax curve (with Hill parameter near 1), but we wish to allow for the possibility of a steeper or shallower curve. The dose with half the effect (ED50) might be near 15 mg, but have considerable uncertainty around that.

model1 <- bf( #---- YOUR CODE HERE --- # ~ E0 + Emax * dose^h/(dose^h + ED50^h),
              family= #---- YOUR CODE HERE --- # ,
              nlf(h ~ exp(logh)), nlf(ED50 ~ exp(logED50)),
              E0 ~ 1, logED50 ~ 1, logh ~ 1, Emax ~ 1, 
              nl=TRUE)

priors1 <- prior(normal(-2.944, 0.38), nlpar=E0) +
           prior(normal(0, log(4)/1.96), nlpar=logh) +
           prior(normal(0, 6), nlpar=Emax) +
           prior(normal(2.7, 2.5), nlpar=logED50)

Now fit the model

brmfit1 <- brm(
  formula = model1,
  prior = priors1, 
  data = summarize(peanut, y=sum(CRIT1FL=="Y"), n=n(), .by=c(TRT01P, dose))
)

Exercise 1: reading the brms code

• y | trials(n) ~ ...: a binomial outcome needs the number of trials. Summarizing patients as y responders out of n per dose gives the same posterior as one row per patient, but fits faster.
• nl=TRUE: the right-hand side is a non-linear formula. Every unknown name (E0, Emax, logh, logED50) is a non-linear parameter with its own formula, here ~ 1: one value shared by all patients.
• nlf(h ~ exp(logh)): defines h from other parameters. Estimating logh and logED50 on the whole real line and exponentiating keeps h and ED50 positive.
• prior(..., nlpar=E0): priors refer to non-linear parameters via nlpar. On the logit scale, E0 is the placebo log-odds and Emax the maximal log-odds ratio vs. placebo.

The priors translate the beliefs described above:

• E0 ~ N(-2.944, 0.38): 95% prior interval [-3.69, -2.20], i.e. 2.4% to 10% placebo responders.
• Emax ~ N(0, 6): a log-odds ratio of 6 (5% vs. 95%) is one SD away: possible, but we are mildly skeptical.
• logh ~ N(0, log(4)/1.96): Hill parameter centered on e^0 = 1 (standard Emax), 95% between 1/4 and 4.
• logED50 ~ N(2.7, 2.5): ED50 centered on e^{2.7} \approx 15 mg, with considerable uncertainty.

Exercise 1: predictions

Now use

tibble(dose = seq(0, 300, 1), n=1) |> tidybayes::add_epred_rvars(brmfit1)

to obtain predicted proportions for every dose from 0 to 300 mg.

Then, plot the curve of predicted proportions for each of these doses using ggplot2 and ggdist (using ydist=.epred in the aesthetics and the stat_lineribbon() geom from the ggdist package).

If you prefer, replace the default fill colors e.g. with shades of blue using + scale_fill_brewer(palette="Blues")

Next obtain the posterior distribution for the difference in proportions for each of these dose levels compared with placebo.

Hint: add_epred_rvars() returns rvar columns, each holding all posterior draws; arithmetic on them (e.g. .epred - .epred_placebo) is carried out draw by draw.

Exercise 2: using more of the data

The primary endpoint throws information away: a patient who tolerated 300 mg and a patient who tolerated nothing are both just “non-responders”. The food challenge actually records, for every patient, the highest amount of peanut protein tolerated (AVAL): one of 0, 3, 10, 100, 300, 600 or 1000 mg.

We model this ordered outcome directly, in two different ways:

• 2a: interval censoring. Each patient has an unobserved, continuous tolerance threshold; we only know between which two challenge steps it lies.
• 2b: ordinal regression. The 7 outcomes are ordered categories; we make no assumption about the distances between them.

Both keep the sigmoid Emax curve from Exercise 1. Only the left-hand side of the formula (the response) and the family (the likelihood) change, which is typical for brms: the model for the mean and the outcome distribution are specified separately.

Both models still let us derive the primary endpoint, the proportion of patients tolerating \geq 600 mg, so at the end we can compare all three analyses on the same scale.

Exercise 2a: interval censoring, the idea

Let T_i be the (unobserved) amount of peanut protein patient i could tolerate. The challenge only tells us which interval T_i lies in:

We assume \log T_i \sim \mathrm{N}(\mu(d_i), \sigma^2) with the Emax model for the mean, \mu(d) = E_0 + E_{\max}\, d^h / (d^h + \mathrm{ED}_{50}^h). A patient with \log T_i \in [L_i, U_i) contributes the probability of that interval to the likelihood:

\Pr(L_i \le \log T_i < U_i) = \Phi\left(\frac{U_i - \mu(d_i)}{\sigma}\right) - \Phi\left(\frac{L_i - \mu(d_i)}{\sigma}\right)

In brms, the cens() addition term on the response requests exactly this: lower | cens("interval", upper) ~ ..., with family = gaussian(). On the log scale, the first and last intervals are open-ended (-\infty to \log 3 and \log 1000 to \infty). The first argument of cens() may also be a variable, so each row can be censored differently: "left" (the response is an upper bound), "right" (the response is a lower bound), "interval" (response and second argument are the lower and upper ends) or "none".

The parameters now live on the log-mg scale instead of the logit scale:

• E_0: mean of \log T on placebo, so e^{E_0} is the median tolerated amount (mg).
• E_{\max}: maximal shift on the log scale; the drug multiplies the median tolerated amount by up to e^{E_{\max}}.
• \sigma (new): SD of \log T between patients on the same dose.

Exercise 2a: your task

The data with one censoring type per patient: left-censored below 3 mg, right-censored above 1000 mg, interval-censored otherwise.

challenge_mg <- c(0, 3, 10, 100, 300, 600, 1000)

peanut_cens <- peanut |>
  mutate(censored = case_when(AVAL == 0 ~ "left",
                              AVAL == 1000 ~ "right",
                              TRUE ~ "interval"),
         # for left-censored rows the response is the upper bound log(3)
         logAVAL = ifelse(AVAL == 0, log(3), log(AVAL)),
         logAVALnext = ifelse(censored == "interval",
                              log(challenge_mg[match(AVAL, challenge_mg) + 1]), NA))

model2 <- bf( #---- YOUR CODE HERE --- # ~ E0 + Emax * dose^h/(dose^h + ED50^h),
              family = #---- YOUR CODE HERE --- # ,
              nlf(h ~ exp(logh)), nlf(ED50 ~ exp(logED50)),
              E0 ~ 1, logED50 ~ 1, logh ~ 1, Emax ~ 1,
              nl=TRUE)

priors2 <- prior( #---- YOUR CODE HERE --- # , nlpar=E0) +
           prior(normal(0, log(4)/1.96), nlpar=logh) +
           prior( #---- YOUR CODE HERE --- # , nlpar=Emax) +
           prior(normal(2.7, 2.5), nlpar=logED50) +
           prior( #---- YOUR CODE HERE --- # , class=sigma)

brmfit2 <- brm(formula = model2, prior = priors2, data = peanut_cens)

Hints:

• The response goes left of |; cens() takes the censoring type and the upper interval end.
• get_prior(model2, data = peanut_cens) lists every parameter that needs a prior, including the new sigma.
• logh and logED50 describe the shape of the dose response and need not change. E0 and Emax are now on the log-mg scale: what is a plausible median tolerated amount on placebo, and how many-fold an increase is plausible?
• After fitting, look at summary(brmfit2) and transform E0 and Emax back with exp().

Exercise 2b: ordinal regression, the idea

Treat AVAL as an ordered outcome Y with K = 7 categories, k = 1 (0 mg) to k = 7 (1000 mg). The cumulative logit model (family = cumulative(link = "logit", threshold = "flexible")) is

\Pr(Y_i \le k \mid d_i) = \text{logit}^{-1}\big(\theta_k - \eta(d_i)\big), \quad k = 1, \ldots, K-1

• \theta_1 < \cdots < \theta_6 are ordered thresholds (cut-points). On placebo, \eta(0) = 0, so \Pr(Y \le k) = \text{logit}^{-1}(\theta_k). These are cumulative probabilities; the probability of one category is a difference, \Pr(Y = k) = \Pr(Y \le k) - \Pr(Y \le k - 1). The thresholds appear as Intercept[1], …, Intercept[6] in brms; "flexible" means they are only constrained to be ordered.
• \eta(d) = E_{\max}\, d^h / (d^h + \mathrm{ED}_{50}^h) is the Emax curve without E_0: adding a constant to \eta and to all \theta_k gives the same probabilities, so E_0 could not be estimated. The thresholds take over its role.
• Mind the minus sign: a larger \eta lowers \Pr(Y \le k), i.e. moves patients to higher categories. A positive E_{\max} means the drug helps.
• Proportional odds: the same \eta(d) shifts every threshold, so the log-odds ratio vs. placebo of Y > k is \eta(d) for every cut-point, including “\geq 600 mg” from Exercise 1. This is an assumption we should check.

Latent variable view: this is the model of 2a with a logistic instead of a normal distribution, scale fixed to 1, and unknown cut-points instead of the known challenge amounts. See Bürkner & Vuorre (2019) (preprint) for details.

Exercise 2b: your task

model3 <- bf( #---- YOUR CODE HERE --- # ~ Emax * dose^h/(dose^h + ED50^h),
              family = #---- YOUR CODE HERE --- # ,
              nlf(h ~ exp(logh)), nlf(ED50 ~ exp(logED50)),
              logED50 ~ 1, logh ~ 1, Emax ~ 1,
              nl=TRUE)

priors3 <- prior(normal(0, log(4)/1.96), nlpar=logh) +
           prior(normal(0, 6), nlpar=Emax) +
           prior(normal(2.7, 2.5), nlpar=logED50) +
           prior( #---- YOUR CODE HERE --- # , class=Intercept)

brmfit3 <- brm(formula = model3, prior = priors3,
               data = mutate(peanut, AVAL = ordered(AVAL)))

Hints:

• The response must be an ordered factor, ordered(AVAL), which creates one level per observed value.
• prior(..., class=Intercept) sets a prior for all thresholds; prior(..., class=Intercept, coef=5) for \theta_5 only.
• Translate placebo beliefs into thresholds: if about 5% of placebo patients tolerate \geq 600 mg, then \Pr(Y \le 300\text{ mg}) = 0.95, i.e. \theta_5 = \text{logit}(0.95) = 2.944, the negative of the E0 prior mean in Exercise 1.
• After fitting, add_epred_rvars() returns, for each dose, the probabilities of all 7 categories.

How would you change the model and the prior distributions compared with Exercise 1?

Exercise 1: Emax logistic regression

model1 <- bf( y | trials(n) ~ E0 + Emax * dose^h/(dose^h + ED50^h),
              family=binomial(link="logit"),
              nlf(h ~ exp(logh)), nlf(ED50 ~ exp(logED50)),
              E0 ~ 1, logED50 ~ 1, logh ~ 1, Emax ~ 1, 
              nl=TRUE)

priors1 <- prior(normal(-2.944, 0.38), nlpar=E0) +
           prior(normal(0, log(4)/1.96), nlpar=logh) +
           prior(normal(0, 6), nlpar=Emax) +
           prior(normal(2.7, 2.5), nlpar=logED50)


brmfit1 <- brm(
  formula = model1,
  prior = priors1, 
  data = summarize(peanut, y=sum(CRIT1FL=="Y"), n=n(), .by=c(TRT01P, dose)),
  refresh=0
)

Running MCMC with 4 parallel chains...

Chain 1 finished in 0.2 seconds.
Chain 2 finished in 0.2 seconds.
Chain 3 finished in 0.2 seconds.
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.2 seconds.
Total execution time: 0.3 seconds.

summary(brmfit1)

 Family: binomial 
  Links: mu = logit 
Formula: y | trials(n) ~ E0 + Emax * dose^h/(dose^h + ED50^h) 
         h ~ exp(logh)
         ED50 ~ exp(logED50)
         E0 ~ 1
         logED50 ~ 1
         logh ~ 1
         Emax ~ 1
   Data: summarize(peanut, y = sum(CRIT1FL == "Y"), n = n() (Number of observations: 6) 
  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
E0_Intercept         -3.20      0.33    -3.86    -2.56 1.00     1651     1915
logED50_Intercept     3.48      1.11     1.88     6.24 1.01      927      925
logh_Intercept       -0.55      0.35    -1.18     0.19 1.00     1113     1774
Emax_Intercept        5.26      1.39     3.44     8.97 1.01      928      990

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).

tibble(dose = seq(0, 300, 1), n=1) |> 
  tidybayes::add_epred_rvars(brmfit1) |>
  ggplot(aes(x=dose, ydist=.epred)) +
  stat_lineribbon() +
  scale_fill_brewer(palette="Blues") +
  ylab("Model predicted proportion") +
  xlab("Dose [mg]")

pr <- tibble(dose = seq(0, 300, 1), n=1) |> 
    tidybayes::add_epred_rvars(brmfit1)

left_join(pr, filter(pr, dose == 0)|> 
              dplyr::select(-dose) |> 
              rename(.pbo=.epred), 
          by="n") |>
    mutate(diff = .epred - .pbo) |>
    ggplot(aes(x=dose, ydist=diff)) +
    stat_lineribbon() +
    scale_fill_brewer(palette="Blues") +
    ylab("Model predicted difference\nin proportion vs. placebo") +
    xlab("Dose [mg]")

Exercise 1: linpred vs. epred vs. predicted

tidybayes offers three ways to get posterior draws for new data. They differ in what is predicted, here for a hypothetical arm of n = 50 patients:

tibble(dose = c(0, 50, 300), n = 50) |>
  add_linpred_rvars(brmfit1, value = "linpred") |>
  add_epred_rvars(brmfit1, value = "epred") |>
  add_predicted_rvars(brmfit1, value = "predicted")

# A tibble: 3 × 5
   dose     n       linpred        epred  predicted
  <dbl> <dbl>    <rvar[1d]>   <rvar[1d]> <rvar[1d]>
1     0    50  -3.20 ± 0.33   2.1 ± 0.66    2 ± 1.5
2    50    50  -0.27 ± 0.19  21.7 ± 2.27   22 ± 4.2
3   300    50   0.71 ± 0.25  33.4 ± 2.67   33 ± 4.3

• add_linpred_rvars() (posterior_linpred()): the non-linear predictor \eta(d) = E_0 + E_{\max}\, d^h / (d^h + \mathrm{ED}_{50}^h) on the link scale, here log-odds. Does not depend on n.
• add_epred_rvars() (posterior_epred()): the expected outcome, E[y \mid d] = n \cdot \text{logit}^{-1}(\eta(d)). It reflects parameter uncertainty only; with n = 1 it is the responder proportion, which is why we used n = 1 above.
• add_predicted_rvars() (posterior_predict()): simulated new outcomes, y \sim \text{Binomial}(n, \text{logit}^{-1}(\eta(d))), i.e. parameter uncertainty plus sampling noise. Hence integer counts and the widest intervals; use it for posterior predictive checks or to predict the result of a future trial.

For families with further distributional parameters, such as the Gaussian sigma in Exercise 2a, add_linpred_rvars(dpar = TRUE) returns these per draw as well.

Exercise 1: logistic E_max vs. BLRM

Exercise 2a: interval censoring, solution

model2 <- bf( logAVAL | cens(censored, logAVALnext) ~ E0 + Emax * dose^h/(dose^h + ED50^h),
              family = gaussian(),
              nlf(h ~ exp(logh)), nlf(ED50 ~ exp(logED50)),
              E0 ~ 1, logED50 ~ 1, logh ~ 1, Emax ~ 1, 
              nl=TRUE )

priors2 <- prior(normal(0, 5), nlpar=E0) +             # centered on 1 mg, very uncertain
           prior(normal(0, log(4)/1.96), nlpar=logh) + # same shape as for binomial outcome
           prior(normal(0, 6.5), nlpar=Emax) +         # >= 1000 mg plausible, 300000 conceivable
           prior(normal(2.7, 2.5), nlpar=logED50) +    # ED50 as for binomial outcome
           prior(lognormal(5, 10), class=sigma)        # new: between-patient SD on log scale

brmfit2 <- brm(
  formula = model2,
  prior = priors2,
  data = peanut_cens,
  refresh=0
  )

Running MCMC with 4 parallel chains...

Chain 1 finished in 4.2 seconds.
Chain 4 finished in 4.2 seconds.
Chain 2 finished in 4.4 seconds.
Chain 3 finished in 5.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 4.5 seconds.
Total execution time: 5.3 seconds.

summary(brmfit2)

 Family: gaussian 
  Links: mu = identity 
Formula: logAVAL | cens(censored, logAVALnext) ~ E0 + Emax * dose^h/(dose^h + ED50^h) 
         h ~ exp(logh)
         ED50 ~ exp(logED50)
         E0 ~ 1
         logED50 ~ 1
         logh ~ 1
         Emax ~ 1
   Data: peanut_cens (Number of observations: 300) 
  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
E0_Intercept          0.83      0.43    -0.02     1.65 1.00     2111     2354
logED50_Intercept     3.16      0.94     1.86     5.61 1.00     1383     1207
logh_Intercept       -0.56      0.31    -1.15     0.08 1.00     1306     1869
Emax_Intercept        7.99      1.81     5.58    12.58 1.00     1206     1257

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     2.64      0.15     2.36     2.94 1.00     2274     2139

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).

Exercise 2a: reading priors and posterior

The priors, translated back to milligrams:

• E0 ~ N(0, 5): placebo median tolerated amount centered on e^0 = 1 mg, with a 95% prior interval of e^{\pm 9.8}, i.e. 0.00006 to 18,000 mg: deliberately vague.
• Emax ~ N(0, 6.5): a one-SD effect multiplies the tolerated amount by e^{6.5} \approx 665, moving a placebo patient from about 1 mg to beyond 600 mg; two SD (e^{13}) remain conceivable.
• logh, logED50: unchanged, because the shape of the dose response does not depend on how we measure the outcome.
• sigma ~ lognormal(5, 10): intended as “allow a large SD on the log scale”. We revisit this choice at the end.

as_draws_df(brmfit2) |>
  mutate_variables(median_pbo_mg = exp(b_E0_Intercept),
                   max_fold_increase = exp(b_Emax_Intercept),
                   ED50_mg = exp(b_logED50_Intercept)) |>
  subset_draws(variable = c("median_pbo_mg", "max_fold_increase", "ED50_mg", "sigma")) |>
  summarise_draws(median, ~quantile2(.x, probs = c(0.025, 0.975)))

# A tibble: 4 × 4
  variable           median    q2.5     q97.5
  <chr>               <dbl>   <dbl>     <dbl>
1 median_pbo_mg        2.30   0.982      5.21
2 max_fold_increase 2021.   266.    291245.  
3 ED50_mg             19.1    6.43     272.  
4 sigma                2.63   2.36       2.94

Emax is the effect at an infinitely high dose, far beyond the 300 mg tested, so the data constrain it only weakly: note the long upper tail of max_fold_increase.

Exercise 2b: ordinal regression, solution

# no E0: the thresholds (Intercept[1], ..., Intercept[6]) describe the placebo outcome
model3 <- bf(AVAL ~ Emax * dose^h/(dose^h + ED50^h),
             family = cumulative(link = "logit", threshold = "flexible"),
             nlf(h ~ exp(logh)),
             nlf(ED50 ~ exp(logED50)),
             logED50 ~ 1,
             logh ~ 1,
             Emax ~ 1, 
             nl=TRUE)

priors3 <- prior(normal(0, log(4)/1.96), nlpar=logh) +
           prior(normal(0, 6), nlpar=Emax) +
           prior(normal(2.7, 2.5), nlpar=logED50) +
           prior(student_t(3, 0, 2.5), class=Intercept) +
           prior(student_t(3, -0.2, 1), class=Intercept, coef=1) +
           prior(student_t(3, 2.197, 1), class=Intercept, coef=5)

brmfit3 <- brm(
  formula = model3,
  prior = priors3,
  data = mutate(peanut, AVAL = ordered(AVAL)),
  refresh=0
  )

Running MCMC with 4 parallel chains...

Chain 1 finished in 10.1 seconds.
Chain 3 finished in 10.5 seconds.
Chain 4 finished in 10.8 seconds.
Chain 2 finished in 11.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 10.6 seconds.
Total execution time: 11.3 seconds.

summary(brmfit3)

 Family: cumulative 
  Links: mu = logit 
Formula: AVAL ~ Emax * dose^h/(dose^h + ED50^h) 
         h ~ exp(logh)
         ED50 ~ exp(logED50)
         logED50 ~ 1
         logh ~ 1
         Emax ~ 1
   Data: mutate(peanut, AVAL = ordered(AVAL)) (Number of observations: 300) 
  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[1]          0.05      0.25    -0.43     0.57 1.00     3574     2409
Intercept[2]          1.01      0.27     0.50     1.58 1.00     4195     2469
Intercept[3]          1.74      0.29     1.19     2.33 1.00     4024     2637
Intercept[4]          2.29      0.30     1.73     2.89 1.00     3930     2519
Intercept[5]          3.04      0.30     2.47     3.65 1.00     4137     2731
Intercept[6]          4.03      0.33     3.40     4.69 1.00     4124     2825
logED50_Intercept     3.82      1.09     2.23     6.38 1.00     1304     1171
logh_Intercept       -0.60      0.30    -1.13     0.00 1.00     1498     2102
Emax_Intercept        5.45      1.50     3.47     9.20 1.00     1225     1119

Further Distributional Parameters:
     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
disc     1.00      0.00     1.00     1.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).

Exercise 2b: reading priors and posterior

• logh, logED50: unchanged (shape of the dose response).
• Emax ~ N(0, 6): under proportional odds Emax is the log-odds ratio for every cut-point, including \geq 600 mg, so the prior from Exercise 1 carries over directly.
• Intercept ~ student_t(3, 0, 2.5): weakly informative default for all thresholds; the ordering constraint does the rest.
• Intercept[1] ~ student_t(3, -0.2, 1): \text{logit}^{-1}(-0.2) \approx 0.45, i.e. about 45% of placebo patients are expected to fail at the first 3 mg step.
• Intercept[5] ~ student_t(3, 2.197, 1): \text{logit}^{-1}(2.197) = 0.9, i.e. 10% placebo responders, whereas Exercise 1 assumed 5%. We revisit this at the end.

Predicted category probabilities on placebo, and the observed proportions:

tibble(dose = 0) |> add_epred_rvars(brmfit3) |> pull(.epred)

rvar<4000>[1,7] mean ± sd:
     0               3               10              100             300             600             1000           
[1,] 0.513 ± 0.0624  0.218 ± 0.0339  0.117 ± 0.0244  0.058 ± 0.0149  0.048 ± 0.0135  0.029 ± 0.0087  0.018 ± 0.0060 

peanut |> filter(dose == 0) |> count(AVAL) |> mutate(observed = n / sum(n)) |> kable(digits = 2)

Comparing the three models: take-aways

• All models agree on the overall picture: responder rates rise steeply up to about 50 mg and flatten afterwards.
• The ordinal model reproduces the binomial dose response with somewhat narrower intervals: it uses all 7 outcome levels, at the price of the proportional odds assumption.
• The censored model predicts fewer responders at high doses, reflecting the stronger assumption that \log T is normally distributed.
• “More efficient” analyses buy precision with assumptions. Check them, e.g. with posterior predictive checks, before trusting the narrower intervals, and pre-specify the primary analysis.

Discussion: the sigma prior (2a)

sigma ~ lognormal(5, 10) was meant to “allow a large SD on the log scale”. What does it actually say?

qlnorm(c(0.05, 0.5, 0.95), meanlog = 5, sdlog = 10)

[1] 1.07e-05 1.48e+02 2.07e+09

plnorm(10, meanlog = 5, sdlog = 10, lower.tail = FALSE)

[1] 0.606

Its median is an SD of about 148 on the log scale, and more than half of its mass lies above 10. Yet already \sigma = 3 means that the middle 95% of patients on one dose span a factor of e^{2 \cdot 1.96 \cdot 3} \approx 10^5. A half-normal prior with scale 3 (normal(0, 3); brms bounds sigma at 0) puts 95% of its mass below 5.9 and is still generous:

priors2_new <- priors2[priors2$class != "sigma", ] + prior(normal(0, 3), class = sigma)
brmfit2_new <- update(brmfit2, prior = priors2_new, refresh = 0)

Running MCMC with 4 parallel chains...

Chain 1 finished in 4.1 seconds.
Chain 2 finished in 4.1 seconds.
Chain 3 finished in 4.4 seconds.
Chain 4 finished in 4.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 4.2 seconds.
Total execution time: 4.6 seconds.

compare_fits <- function(fits, variable) {
  lapply(fits, \(fit) summarise_draws(as_draws_df(fit, variable = variable),
                                      mean, sd, ~quantile2(.x, probs = c(0.025, 0.975)))) |>
    bind_rows(.id = "fit") |>
    arrange(variable)
}

compare_fits(list(original = brmfit2, `half-normal` = brmfit2_new),
             c("b_E0_Intercept", "b_Emax_Intercept", "sigma")) |>
  kable(digits = 2)

The posterior hardly changes: with 300 patients the data dominate sigma. The original prior was harmless here, but would not be in a small trial or a prior predictive check. Always translate priors back to the outcome scale.

Discussion: the threshold prior (2b)

Exercise 1 assumed about 5% placebo responders (E0 ~ N(-2.944, 0.38) on the logit scale), but the ordinal model centered Intercept[5] on 10%. The same belief about the same trial should give the same prior, i.e. \theta_5 centered on \text{logit}(0.95) = 2.944:

priors3_new <- priors3[!(priors3$class == "Intercept" & priors3$coef == "5"), ] +
  prior(student_t(3, 2.944, 1), class = Intercept, coef = 5)
brmfit3_new <- update(brmfit3, prior = priors3_new, refresh = 0)

Running MCMC with 4 parallel chains...

Chain 3 finished in 9.7 seconds.
Chain 2 finished in 9.7 seconds.
Chain 4 finished in 11.3 seconds.
Chain 1 finished in 12.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 10.7 seconds.
Total execution time: 12.2 seconds.

compare_fits(list(original = brmfit3, `5% responders` = brmfit3_new),
             c("b_Intercept[5]", "b_Emax_Intercept")) |>
  kable(digits = 2)

Again little changes: the 50 placebo patients, and through proportional odds all other patients, inform the thresholds well. Still, keep priors consistent across alternative analyses of the same trial; otherwise differences between the analyses may be driven by the priors rather than by the models.