Exercises Hierarchical Models in Preclincal Research
+
+
+
+
+
+
+
+
+
+
Setup
+
+
+
+
R Session
+
+
Note: All shown code can be copied by simply clicking the top right icon of the code block. The code blocks from the main presentation can also be copied by the same icon, which is visible when hovering over the respective block with your mouse.
+
+
For the course a RStudio server has been setup by the University of Göttingen: http://134.76.20.1:8787
+
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 is better described as part of the respective cmdstan guide chapter.
+
C++ toolchain installation
+
+
Windows: cmdstanr can automatically install the RTools based C++ toolchain via
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", "RBesT", "ggrepel", "patchwork", "ggdist", "withr"))
+
+
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 users 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 you 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)
+
+# 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(size=1.5))
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
+
Create a posterior predictive check based on the predictive distribution for the response rate.
+
Steps:
+
+
Use posterior_predict to create samples from the predictive distribution of outcomes per trial.
+
Use sweep(predictive, 2, AS_region$n, "/") to convert these samples from the predictive distribution of the outcome counts to samples from the predictive distribution for responder rates.
+
Use ppc_intervals from bayesplot to create a plot showing your results. Consider using the with function to refer more easily to the data columns of the analysis data set.
+
Redo the above using from the tidybayes package the add_predicted_rvars function adding the predictions directly to the analysis data set and then redo task 2 and 3. To convert a rvars column to the respective draws format of bayesplot you can use the as_draws_matrix from the posterior package.
+
+
+
+
Exercise 2: Fixed vs random effect
+
Redo the analysis with region, but treat region as a fixed effect. Evaluate the informativeness of the obtained MAP priors. The model formula for brms should look like
+
+
region_model_fixed <-bf(r |trials(n) ~1+ region + (1| study), family=binomial)
+
+
Steps:
+
+
Consider the prior for the region fixed effect first. The reference region is included in the intercept. The reference region is implicitly defined by the first level of the variable region when defined as factor.
+
+
Define asia to be the reference region in the example. Also include a level other in the set of levels.
+
Assume that an odds-ratio of \(2\) between regions can be seen as very large such that a prior of \(\text{N}(0, (\log(2)/1.96)^2)\) for the region main effect is adequate.
+
+
Obtain the MAP prior for each region by using the AS_region_all data frame defined below and apply posterior_linpred as shown in the case study.
+
Convert the MCMC samples from the MAP prior distribution into mixture distributions with the same code as in the case study.
+
Calculate the effective sample size (ESS) for each prior distribution with the ess function from the RBesT R package.
Run the analysis for the normal endpoint in the crohn data set of RBesT. Refer to the RBesT vignette for a normal endpoint on more details and context.
+
Steps:
+
+
Use as family=gaussian and use the se response modifier in place of trials to specify a known standard error. More details on the additional response information specifiers can be found for the documentation of brmsformula.
+
Use the same priors as proposed in the vignette from RBesT.
+
Compare the obtained MAP prior (in MCMC sample form) from RBesT and brms.
# 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.3 ± 5.7 0.23 ± 0.053
+2 Study 2 44 12 asia 11.4 ± 3.6 0.26 ± 0.081
+3 Study 3 51 19 europe 16.1 ± 4.4 0.32 ± 0.086
+4 Study 4 39 9 europe 9.3 ± 3.3 0.24 ± 0.085
+5 Study 5 139 39 asia 37.8 ± 6.9 0.27 ± 0.049
+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)
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 sd 2.5% 50% 97.5%
+14.60 9.65 1.47 12.70 39.10
+
+MAP Prior MCMC sample
+ mean sd 2.5% 50% 97.5%
+-49.6 19.6 -91.0 -48.4 -11.0
+
+
brms_normal_map_mcmc
+
+
Family: gaussian
+ Links: mu = identity; sigma = 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)))
This exercises accompanies 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.
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 log-scale close to 6), but we are at least mildly skeptical.
+
The dose response is a-priori somewhat likely to follow a 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.
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.
+
+
+
Excercise 2: ordinal outcome and interval censored
+
Try to perform a more efficient statistical analysis by either using an ordinal outcome or treating the data as interval censored.
+
Interval censored
+
Try treating the data about the logarithm of the amount of protein tolerated as interval censored (\(-\infty\) to \(< \log(3)\), \(>= \log(3)\) to \(< \log(10)\), …) using logAVAL | cens("interval", logAVALnext) as endpoint and with family = gaussian(). For the latter approach you might have to set log(0) to a very low number such as -28 (approximate \(\log(\text{weight})\) of one single peanut protein molecule) and assume the upper interval end when a patient could tolerate 1000 mg to, say, \(\log(30000)\) (approximately 100 peanut kernels). This achieves approximately the same thing as setting these to \(-\infty\) and \(\infty\), but is handled better by brms/Stan.
+
How would you change the prior distributions in this case?
+
Ordinal outcome aproach
+
Now use an ordinal outcome via family = cumulative(link = "logit", threshold = "flexible"). Details of this model are described in a journal article, for which there is also a PsyArXiv preprint in case you do not have access to the journal.
+
Note that the placebo response is described by a series of \(K-1\) ordered thresholds \(\theta_1, \ldots, \theta_{K-1}\) on the logit scale for \(K\) ordered categories. For the placebo group, the probability of being in categories \(k=1, \ldots, K-1\) is be given by \(\text{logit}^{-1} \theta_k\) and for category \(K\) by \(1-\text{logit}^{-1} \theta_{K-1}\). In this case, these parameters will appear in summary(brmfit3) as Intercept[1], …, Intercept[6].
+
How would you change the model and the prior distributions in this case?
+
Note, if we had a baseline amount of protein tolerated, we could treat this as a monotonic covariate e.g. using mo(BASE) or assume a particular functional form.
Family: gaussian
+ Links: mu = identity; sigma = identity
+Formula: logAVAL | cens("interval", logAVALnext) ~ E0 + Emax * dose^h/(dose^h + ED50^h)
+ h ~ exp(logh)
+ ED50 ~ exp(logED50)
+ E0 ~ 1
+ logED50 ~ 1
+ logh ~ 1
+ Emax ~ 1
+ Data: mutate(peanut, logAVAL = ifelse(AVAL == 0, -28, lo (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.90 0.40 0.14 1.67 1.00 2174 2058
+logED50_Intercept 3.11 1.01 1.75 5.66 1.00 993 926
+logh_Intercept -0.63 0.34 -1.24 0.07 1.00 1007 1582
+Emax_Intercept 7.60 1.73 5.28 11.99 1.00 764 1066
+
+Further Distributional Parameters:
+ Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
+sigma 2.43 0.12 2.19 2.68 1.00 2200 2326
+
+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).
+
+
+
Ordinal outcome aproach
+
+
# Note that we omit E0 from the model, because the expected placebo outcome
+# is already described by the intercept parameters of the distribution
+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, 1), nlpar=logh) +# No particular reason to
+# change prior here
+prior(normal(0, 6), nlpar=Emax) +# Prior could be
+# different (odds from
+# one category to the
+# next not necessarily
+# the same as for top two
+# vs. rest)
+prior(normal(2.7, 2.5), nlpar=logED50) +# No particular
+# reason to change
+# prior here
+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 4 finished in 9.9 seconds.
+Chain 2 finished in 10.3 seconds.
+Chain 3 finished in 11.8 seconds.
+Chain 1 finished in 12.7 seconds.
+
+All 4 chains finished successfully.
+Mean chain execution time: 11.2 seconds.
+Total execution time: 12.8 seconds.
+
+
summary(brmfit3)
+
+
Family: cumulative
+ Links: mu = logit; disc = identity
+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.07 0.25 -0.41 0.54 1.00 3641 2673
+Intercept[2] 1.02 0.26 0.53 1.54 1.00 4164 2782
+Intercept[3] 1.75 0.27 1.22 2.30 1.00 4376 2983
+Intercept[4] 2.29 0.28 1.75 2.85 1.00 4357 3066
+Intercept[5] 3.05 0.29 2.49 3.63 1.00 4469 3287
+Intercept[6] 4.04 0.32 3.43 4.67 1.00 4795 2844
+logED50_Intercept 3.96 1.16 2.26 6.76 1.00 1356 1160
+logh_Intercept -0.64 0.30 -1.18 -0.03 1.00 1496 2133
+Emax_Intercept 5.66 1.58 3.57 9.62 1.00 1314 1195
+
+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).
The opinions & statements expressed in this presentation and on the following slides are solely of the presenters, and not necessarily those of Novartis.
+
+
+
Learning Goals
+
After this course, you should:
+
+
Be familiar with brms syntax and workflow
+
Recognize its versatility for statistical modelling in drug development
+
Have hands-on experience with the package from two guided exercises
+
+
and of course:
+
+
Feel empowered to use brms going forward!
+
+
+
+
Housekeeping
+
+
Q&A: you may raise your hand at any time, or hold for Q&A sessions at the end of each section
+
Laptop charging: we recommend conserving battery by keeping your laptop powered down except during the hands-on exercises
Internally write Stan code that is readable yet fast
+
Provide an easy interface for defining priors & models
+
Facilitate post-processing
+
+
+
+
+
+
+
Some Highlights of brms
+
+
Flexible hierarchical (random effects) modeling
+
Both built-in and user-defined likelihoods
+
Explicit and implicit non-linear modeling
+
Distributional regression
+
Within-chain parallelization
+
Posterior and prior predictions
+
Model checking
+
Highly dense feature matrix
+
+
+
+
Model specification via formula
+
Varying intercept model for a single grouping factor:
+
+
formula <- y ~1+ x + (1| g)
+
+
Varying intercept-slope model for a single grouping factor:
+
+
formula <- y ~1+ x + (1+ x | g)
+
+
Advanced non-linear terms such as Gaussian processes:
+
+
formula <- y ~1+gp(x) + (1+ x | g)
+
+
… or approximate Gaussian processes:
+
+
formula <- y ~1+gp(x, k=9) + (1+ x | g)
+
+
+
+
Model specification via formula
+
Linear formulas for multiple distributional parameters (e.g., predict mean and overdispersion of negative binomial):
+
+
formula <-bf(
+ y ~1+ x + (1| g) + ...,
+ par2 ~1+ x + (1| g) + ...,
+ par3 ~1+ x + (1| g) + ...,
+)
+
+
Non-linear formula for a single distributional parameter:
+
+
formula <-bf(
+ y ~fun(x, nlpar1, nlpar2),
+ nlpar1 ~1+ x + (1| g) + ...,
+ nlpar2 ~1+ (1| g) + ...,
+nl =TRUE
+)
+
+
+
+
Model likelihood via family
+
General structure:
+
+
family <-brmsfamily(
+family ="<family>", link ="<link>",
+ more_link_arguments
+)
+
+
+
Gaussian likelihood (default):
+
+
family <-brmsfamily(family ="gaussian", link ="identity",
+link_sigma ="log")
+
+
Poisson likelihood:
+
+
family <-brmsfamily(family ="poisson", link ="log")
+
+
See also vignette("brms_families") for details on the families.
+
+
+
+
Recommended global brms settings
+
+
+
+
Ensure to have at least 8 GB of RAM available for compilation step
+
Preferable to use the cmdstanr backend and set cmdstanr_write_stan_file_dir to a model binary caching directory. This avoids model recompilation with and between R sessions.
+
Use of the here R package is encouraged to manage file paths. This makes all file paths relative to a common project root.
+
+
+
+
library(brms)
+here::i_am("relative_path_in_project_to/rscript.R")
+library(here)
+
+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
+)
+# create cache directory if not yet available
+dir.create(here::here("_brms-cache"), FALSE)
+
+
+
+
+
+
Part 1: historical control data
+
+
+
+
Use of historical control data in clinical trials
+
+
Goal: reduce control group sample size while maintaining power
+
+
+
+
Collect historical data from relevant literature systematically
+
Evaluate heterogeniety of historical data
+data quality, patient population, trial design
Borrow from historical data under an exchangeability assumption
+
Discount/down-weight historical data through between-trial heterogeneity \tau
+
Meta-analytic-predictive (MAP) prior through generative nature of hierarchical models
+
+
+
+
Pre-requisites
+
+
essentially the same endpoint
+
similar treatment procedures
+
similar patient eligibility criteria
+
similar patient population characteristics
+
+
+
+
+
+
+
+
+
+
flowchart LR
+ M((β,τ)) --> A
+ M --> B
+ M --> C
+ M -->|prediction| T
+ A((θ<sub>1</sub>)) --> YA[Y<sub>1</sub>]
+ B((θ<sub>2</sub>)) --> YB[Y<sub>2</sub>]
+ C((θ<sub>3</sub>)) --> YC[Y<sub>3</sub>]
+ T((θ<sub>*</sub>)) -->|new trial| YT[Y<sub>*</sub>]
+
Historical control data
+8 studies with 513 patients
+
Bayesian design
+
+
P(\pi_t - \pi_p > 0 | y) > 0.95
+
Placebo: (“MAP prior” from historical data)
+
+
\pi_{p} \sim \text{Beta}(11,32)
+
+
Active:
+
+
\pi_{t} \sim \text{Beta}(0.5,1)
+
+
+
Randomization ratio 4:1 (24 vs 6)
+
+
+
+
+
+
+
+
+
+
+
+
MAP prior:
+
+
+
+
+
+
+
+
+
mean
+
sd
+
2.5%
+
50%
+
97.5%
+
+
+
+
+
0.26
+
0.09
+
0.11
+
0.25
+
0.47
+
+
+
+
+
+
+
+
+
+
+
Generalized Meta-Analytic-Predictive model
+
Y is the (control) group summary data for H historical trials \begin{align*}
+Y_{h}|\theta_{h} &\sim f(\theta_{h}) & \forall \; h \in [1,H] \\
+Y_\ast|\theta_\ast &\sim f(\theta_\ast) & \text{for new trial}
+\end{align*}
Region specific MAP Priors via varying regions model
+
+
form_AS_region <-bf(r |trials(n) ~1+ (1| region/study),
+family = binomial, center=FALSE)
+
+# show which priors brms defines for this model
+get_prior(form_AS_region, AS_region)
+
+bprior_AS_region <-prior(normal(0, 2), class=b, coef=Intercept) +
+prior(normal(0, 0.50), class=sd, coef=Intercept, group=region) +
+prior(normal(0, 0.25), class=sd, coef=Intercept, group=region:study)
+
+fit_AS_region <-brm(
+ form_AS_region, data = AS_region, prior = bprior_AS_region, seed =29856341,
+control=list(adapt_delta=0.99))
+
+
+
(1 | region/study) specifies a nested random effect structure as in this example each study is run in a single region. It is equivalent to
+(1 | region) + (1 | region:study).
+
+
+
+
Summary region model
+
+
summary(fit_AS_region)
+
+
Family: binomial
+ Links: mu = logit
+Formula: r | trials(n) ~ 1 + (1 | region/study)
+ Data: AS_region (Number of observations: 8)
+
+Multilevel Hyperparameters:
+~region (Number of levels: 3)
+ Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
+sd(Intercept) 0.27 0.22 0.01 0.82 1.00 1570 1731
+
+~region:study (Number of levels: 8)
+ Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
+sd(Intercept) 0.23 0.13 0.01 0.51 1.00 1675 1548
+
+Regression Coefficients:
+ Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
+Intercept -1.08 0.26 -1.63 -0.55 1.00 1914 1889
pe_mix_region <- RBesT::automixfit(pe_region[, 1], type ="beta")
+plot(pe_mix_region)$mix
+
+
+
+
+
Summary historical control data
+
+
Bayesian random-effects meta-analysis models can be used to derive Meta-Analytic-Predictive (MAP) priors
+
Predictions for new study mean inform the MAP prior
+
Specification and inference for MAP models is simple in brms
+
Not covered: Parametric MAP mixtures can be used as priors in brms using the RBesT::mixstanvar adapter
+
+
+
+
Hands-on exercises: historical control data
+
+
+
+
+
+
Priors
+
+
+
+
+
The prior can only be understood in the context of the model 1
+
+
+
+
+
+
+
+
+
+
+
+
Commonly heard reasons to be afraid of priors
+
“Priors are inherently subjective”
+
“Priors bias your analysis”
+
“You should let the data speak for themselves”
+
“Non-Bayesian models don’t have priors and don’t need them”
+
“I have no idea how to set appropriate priors”
+
+
+
Are priors inherently subjective?
+
Consider using the alternative terms “contextual” or “conditional”1
+
+
on the person or people or agent performing the analysis
+
on results of previous studies
+
on structural assumptions
+
…
+
+
+
+
Do priors bias your analysis?
+
Definition:
+
An estimator \hat{\theta} of the parameter \theta is called unbiased if \text{E}_{y}(\hat{\theta}) = \theta.
+
Corrollary:
+
If an estimator \hat{\theta} is unbiased for a given prior, it will be biased for most other priors.
+
Corrollary:
+
If an estimator \hat{\theta} is unbiased at a given scale, it will be biased for most non-linear transformations f(\hat{\theta}) (Jensen’s inequality):
+\text{E}_y(f(\hat{\theta})) \neq f(\theta)
+
+
Unbiasedness is fragile and often only holds in the limit
+
+
+
Is bias even the right metric?
+
A sampling distribution is a distribution of a data-dependent variable (a “statistic”) as implied by variation in data y across the true data generating process
+
Consider two estimators \hat{\theta}_1 and \hat{\theta}_2 of \theta with sampling distributions given by
The RMSE can be defined as \begin{align*}
+\text{RMSE}_y(\hat{\theta}, \theta) & = \sqrt{\int (\hat{\theta}(y) - \theta)^2 \; p(y) \; d y} \\
+& = \sqrt{\text{Var}_y{\hat{\theta}} + \text{Bias}_y(\hat{\theta}, \theta) ^2}
+\end{align*}
+
It expresses a trade-off between variance and bias
+
A lot of priors will help you win the bias-variance trade-off, not the race for minimal bias!
+
+
+
Non-Bayesian models don’t have priors?
+
Basic multilevel model with J groups:
+
+
y ~normal(b0[j] + b1 * x, sigma)
+b0[j] ~normal(b0, sd0)
+
+
Is b0[j] ~ normal(b0, sd0) a likelihood, prior, or something else?
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Some reasons for using priors
+
+
Make a-priori implausible values unlikely (weakly informative priors)
+
Incorporate specific expert information into the model (“subjective” priors)
Biologically plausible values: approx. .25 < h < 4 (see prior interval probabilities).
+The above prior puts roughly 1/3 above 4 and 1/3 below 0.25!
+
Re-writing the fraction to \left(\frac{d}{\text{ED}_{50}}\right)^h / \left(1+\left(\frac{d}{\text{ED}_{50}}\right)^h\right) makes it obvious that for h \rightarrow 0, this goes to 1/2, and for h \rightarrow \infty, either to 0 (d < \text{ED}_{50}) or 1 (d > \text{ED}_{50}).
+Computed from 4000 by 4 log-likelihood matrix.
+
+ Estimate SE
+elpd_loo 0.1 0.7
+p_loo 2.1 0.5
+looic -0.3 1.3
+------
+MCSE of elpd_loo is NA.
+MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 0.6]).
+
+Pareto k diagnostic values:
+ Count Pct. Min. ESS
+(-Inf, 0.7] (good) 2 50.0% 150
+ (0.7, 1] (bad) 2 50.0% <NA>
+ (1, Inf) (very bad) 0 0.0% <NA>
+See help('pareto-k-diagnostic') for details.
+
+
+
+
brms interfaces with the loo R package that performs efficient approximate leave-one-out cross-validation (LOO) for Bayesian models by Pareto smoothed importance sampling (PSIS). See Vethari et al. (2017) and the loo documentation on CRAN.
+Based on 4-fold cross-validation.
+
+ Estimate SE
+elpd_kfold -2.1 1.7
+p_kfold 4.4 2.2
+kfoldic 4.2 3.4
+
+
+
+
elpd_kfold: CV estimate of the expected log pointwise predictive density for a new dataset
+
p_kfold: CV estimate of difference between elpd_loo and the non-cross-validated log posterior predictive density (describes how much more difficult it is to predict future data than the observed data).
+
kfoldic: -2 * elpd_kfold (on the conventional scale of deviance)
loo_compare computes a paired estimate of the difference in ELPD to take advantage of the fact that the same set of data points was used to fit both models (and reduce the standard error). See Vethari et al. (2017) and the loo R package manual for details.
Mixed Effects Model for Repeated Measures (MMRM) commonly used for longitudinal data (each patient measured at multiple visits)
+
Direct likelihood analysis that can address hypothetical estimand of all patients completing the study on treatment without doing missing data imputation first
+
+
+
+
Commonly no structure assumed for correlations between visits and different variance allowed for different visits (to avoid unnecessary assumptions)
+
Convergence issues with REML fit common, especially for small sample sizes, which is alleviated by Bayesian inference with (weakly-)informative priors
+
+
+
+
+
Bayesian approach allows us to incorporate prior information and historical data, which is very interesting for Phase I studies
+
brms lets us easily add more components & structure to the model
+
+
+
+
+
What do our data look like?
+
+
+
+
Analysis Data Model (ADaM) Basic Data Structure (BDS)
+
+
+
+
+
+
+
+
+
USUBJID
+
TRT01P
+
AVISIT
+
ADY
+
AVAL
+
CHG
+
BASE
+
+
+
+
+
3
+
10
+
visit1
+
14
+
1.32
+
0.54
+
0.78
+
+
+
3
+
10
+
visit2
+
28
+
1.30
+
0.52
+
0.78
+
+
+
3
+
10
+
visit3
+
56
+
−0.24
+
−1.02
+
0.78
+
+
+
3
+
10
+
visit4
+
84
+
1.40
+
0.63
+
0.78
+
+
+
9
+
10
+
visit1
+
14
+
0.95
+
0.84
+
0.11
+
+
+
9
+
10
+
visit2
+
28
+
1.46
+
1.35
+
0.11
+
+
+
9
+
10
+
visit3
+
56
+
−0.97
+
−1.08
+
0.11
+
+
+
9
+
10
+
visit4
+
84
+
1.21
+
1.10
+
0.11
+
+
+
13
+
10
+
visit1
+
14
+
2.35
+
0.45
+
1.90
+
+
+
13
+
10
+
visit2
+
28
+
2.04
+
0.14
+
1.90
+
+
+
+
+
+
+
+
+
Patients can miss visits or drop-out of the study prior to the final visit.
+
+
+
+
Model Specification: Informal (popular choice)
+
+
+
Fixed effects:
+
+
baseline (pre-treatment) value of the endpoint
+
visit as a factor
+
treatment as a factor
+
treatment by visit interaction
+
visit by baseline value interaction
+
+
+
Random effects:
+
+
Random subject effect on the visit main effect
+
… or equivalently correlated residual error terms within subjects
+
+
+
+
+
Model Specification: Formal
+
Formally, let us assume that there are V visits. We assume that the V-dimensional response \boldsymbol{Y}_i for patient i satisfies
+\boldsymbol{Y}_i = \boldsymbol{X}_i \boldsymbol{\beta} + \boldsymbol{Z}_i \boldsymbol{b}_i + \boldsymbol{\epsilon}_i
+
+
with \boldsymbol{b}_i \sim \text{MVN}(\boldsymbol{0}, \boldsymbol{D}) and \boldsymbol{\epsilon}_i \sim \text{MVN}(\boldsymbol{0}, \boldsymbol{\Sigma}), where \boldsymbol{\Sigma} is a diagonal matrix. This implies
+\boldsymbol{Y}_i \sim \text{MVN}(\boldsymbol{X}_i \boldsymbol{\beta}, \boldsymbol{V}_i),
+ where \boldsymbol{V}_i = \boldsymbol{Z}_i \boldsymbol{D} \boldsymbol{Z}_i^T + \boldsymbol{\Sigma}. Model the correlated Y_{ij} either by
+
+
marginalizing out random effects & account for them with correlated residual errors (residual covariance matrix V_i), or
+
conditionally on (V-dimensional) random effects \boldsymbol{b}_i with residual errors \boldsymbol{\epsilon}_i independent (once we condition on \boldsymbol{b}_i).
+
+
+
+
MMRMs in SAS
+
Widely-used high-quality reference implementation
+
+
PROC MIXED DATA=simulated_data;
+ CLASS TRT01P AVISIT USUBJID;
+ MODEL CHG ~ TRT01P AVISIT BASE TRT01P*AVISIT AVISIT*BASE
+ / SOLUTION DDFM=KR ALPHA = 0.05;
+ REPEATED AVISIT / TYPE=UN SUBJECT = USUBJID R Rcorr GROUP=TRT01P;
+ LSMEANS TRT01P*AVISIT / DIFFS PDIFF CL OM E;
+RUN;
+
+
+
+
MMRMs with the mmrm R package
+
Fit the model in R in a frequentist framework
+
+
library(mmrm)
+mmrm_fit <-mmrm(
+formula = CHG ~ TRT01P + AVISIT + BASE + AVISIT:TRT01P +
+ AVISIT:BASE +us(AVISIT | TRT01P / USUBJID),
+method ="Kenward-Roger",
+vcov ="Kenward-Roger-Linear", # to match SAS
+data =mutate(simulated_data, USUBJID=factor(USUBJID))
+)
To work with MCMC samples of the expected marginal means & to summarize them exactly as we want (e.g. quantile credible intervals instead of highest posterior density intervals) we can use:
More on estimands, parametrization, contrasts & priors
+
Estimating average differences across visits
+
Meta-analytic combined (MAC) using historical data
+
Robustifying MAC via a “slab-and-spike”-type prior
+
Non-linear functions over time & doses in MMRMs
+
Excercises
+
+
+
+
+
Part 4: Longitudinal data
+
+
+
+
+
Longitudinal Data: Background
+
+
Example: simulated study of a drug for the treatment of Psoriasis
+
112 patients, 59 randomized to placebo, 53 to drug.
+
Efficacy was assessed using the Psoriasis Area and Severity Index (PASI), a numerical score which measures the severity and extent of psoriasis (0: no disease, 72: maximum).
+
Assessed at baseline and 7 post-baseline timepoints.
pp_check in brms performs posterior predictive checks with the help of the bayesplot package (shows dens_overlay by default).
+
+
Can also be called with type loo_pit (probability integral transformation). 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). See chapter 12 of BAMDD for details.
+
+
+
+
+
Mean Model: LOO-CV results
+
+
loo(fit_long_mean)
+
+
+Computed from 4000 by 684 log-likelihood matrix.
+
+ Estimate SE
+elpd_loo -2285.7 33.1
+p_loo 101.9 9.5
+looic 4571.4 66.3
+------
+MCSE of elpd_loo is NA.
+MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.8]).
+
+Pareto k diagnostic values:
+ Count Pct. Min. ESS
+(-Inf, 0.7] (good) 681 99.6% 103
+ (0.7, 1] (bad) 3 0.4% <NA>
+ (1, Inf) (very bad) 0 0.0% <NA>
+See help('pareto-k-diagnostic') for details.
+
+
+
+
+
While this is quick, is this what we are interested in?
+
+
If we’re interested in the prediction for a new (unseen) patient, we should be running leave-one-patient-out CV (not leave-one-observation-out)…
+See kfold_split_grouped in the loo package.
+
For the Emax model data: not possible (IPD needed!)
A more formal model comparison may be theoretically advantageous but given that we add only a single parameter (the AR1 effect), we should probably gain more than we loose by including it.
+
+
+
Gaussian Process AR1 Model: Visualization of effects
