Using a Mixture Prior for E0 in the Bayesian Emax Model

Overview

We demonstrate how to compute and apply a mixture prior for placebo response using summary data from internal historical studies or publications. The derived prior is used to supplement information for \(E_0\) in Emax model fitting. The examples have two parts:

The mixture prior is developed using RBesT package, and then it is applied to \(E_0\) using the clinDR function fitEmaxB. There are examples with binary and continuous responses. Note that covariate adjustment is not implemented with mixture PBO prior distributions.

Binary Response

Deriving a Mixture Prior Using RBesT

This example is from Getting Started with RBesT (binary). Additional details for the mixture prior are in this primary reference.

Data

The dataset has control group information from historical trials in Ankyloising Spondelitis [1]. The primary efficacy endpoint is a binary response indicating >=20% change from baseline in the Assessment of Spondelitis Arthritis International Society criteria for improvement (ASAS20) at Week 6:

kable(AS, digits=2, caption="Literature Data: ASAS20 Responder Summary at Week 6") %>% 
  kable_styling(full_width = FALSE, position="center")
Literature Data: ASAS20 Responder Summary at Week 6
study n r
Study 1 107 23
Study 2 44 12
Study 3 51 19
Study 4 39 9
Study 5 139 39
Study 6 20 6
Study 7 78 9
Study 8 35 10

Computing a MAP Prior MCMC Sample using gMAP

An MCMC sample from the Meta Analytic Predictive (MAP) prior can be generated using the gMAP() function in RBesT. For binary endpoints, the prior distribution is specified on the logit scale. A conservative choice of the between trial heterogeneity is Half-Normal(0, 1). A Normal(0,2) prior distribution is specified for the mean of the logit rates. The prior mean of the logit rates may be assigned a more substantive value than 0, which corresponds to a rate of 50 percent. For more information about gMAP use ?gMAP(), and for a detailed description of the statistical methodology refer to [2].

set.seed(34563)
map_mcmc <- gMAP(cbind(r, n - r) ~ 1 | study,
  data = AS,
  tau.dist = "HalfNormal",
  tau.prior = 1,
  beta.prior = 2,
  family = binomial
)
print(map_mcmc)
## Generalized Meta Analytic Predictive Prior Analysis
## 
## Call:  gMAP(formula = cbind(r, n - r) ~ 1 | study, family = binomial, 
##     data = AS, tau.dist = "HalfNormal", tau.prior = 1, beta.prior = 2)
## 
## Exchangeability tau strata: 1 
## Prediction tau stratum    : 1 
## Maximal Rhat              : 1 
## 
## Between-trial heterogeneity of tau prediction stratum
##         mean median    sd   q2.5   q50 q97.5
## tau[1] 0.387  0.359 0.215 0.0399 0.359 0.897
## 
## MAP Prior MCMC sample
##                  mean median     sd  q2.5   q50 q97.5
## theta_resp_pred 0.255  0.247 0.0864 0.106 0.247 0.461

Parametric Approximation Derived from the MCMC Sample

The current representation in clinDR allows a normal mixture for the prior for \(E_0\). This is a prior for placebo response on the logit scale. Therefore, it is required to extract the MAP MCMC sample on the logit scale.

The next step is to convert the MCMC sample to a parametric representation with the automixfit() function that computes the optimal number of mixture components using the AIC. In practice, the optimal number of mixture components is often large. Therefore, we recommend using mixfit() with the number of components selected manually. For the protocol and Statistical Analysis Plan (SAP), 2-3 component mixtures are preferred over the optimal selection strategy. The optimal and manual selection performance can be visually compared.

# Extract posterior samples on the logit scale
post_samp <- as.matrix(map_mcmc)[,'theta_pred']

map <- mixfit(post_samp,Nc=2)

print(map)
## EM for Normal Mixture Model
## Log-Likelihood = -2464.067
## 
## Univariate normal mixture
## Mixture Components:
##   comp1      comp2     
## w  0.6293999  0.3706001
## m -1.1044228 -1.1491835
## s  0.2706598  0.7031189
## check accuracy of mixture fit
plot(map)$mix
plot of chunk mixmap
plot of chunk mixmap

Robustification and Effective Sample Size

We recommend robustifying the MAP prior, which helps to protect against type-I error inflation in case of prior-data conflict. The unit-information prior \(logit(p)\) [3] is: \[logit(p) \sim Normal(mean=logit(p_0), SD=[p_0 (1-p_0)]^{-1/2})\] where \(p_0\) is the mean of MAP samples.

p0 <- plogis(mean(post_samp))
sdp0 <- (p0*(1-p0))^(-1/2)
map_robust <- robustify(map, weight = 0.5, mean = qlogis(p0), sigma=sdp0)
print(map_robust)
## Univariate normal mixture
## Mixture Components:
##   comp1      comp2      robust    
## w  0.3146999  0.1853001  0.5000000
## m -1.1044228 -1.1491835 -1.1210111
## s  0.2706598  0.7031189  2.3224781

The effective sample size (ESS) contributed by the historical data can be calculated with the ess() function.

round(ess(map, method = 'elir', sigma=sdp0))
## [1] 36

Fit Binary Bayesian Emax Model with Mixture Prior for E0

A simulated data set is created to illustrate the method.

# Dose level
dose <- c(0, 5, 25, 50 , 100)
# Sample size
ss <- c(38, 77, 79, 82, 80)
resplev<-c(0.25, 0.40, 0.65, 0.72, 0.75) 
y<-rbinom(length(ss),ss,resplev)
y
## [1] 12 32 58 60 61

The mixture prior for \(E_0\) is specified in the function emaxPrior.control()[4][5][6]:

Other non-default parameters (difTargetmu, difTargetsca, dTarget, p50) must also be provided. Any values specified for epmu and epsca are ignored if mixP>=1.

To fit the binary \(E_{max}\) model using fitEmaxB(), y must be 0/1 and counts must be supplied for each 0/1 value.

# y in grouped data 0/1
cnt<-c(ss-y,y)
y <- c(rep(0, length(dose)), rep(1, length(dose)))
# Dose specification for grouped data
dose <- rep(dose,2)

ed50 <- 15
dTarget <- max(dose)
difTarget <- qlogis(0.75) - qlogis(0.25)

# E0 prior and other hyper priors
prior_bin <- emaxPrior.control(epmu=-1.1, epsca=1,
                               mixP=3, 
                               w_ep= map_robust[1,],
                               mu_ep=map_robust[2,],
                               sd_ep = map_robust[3,],
                               difTargetmu = difTarget,
                               difTargetsca= 1,
                               dTarget = dTarget,
                               p50 = ed50,
                               binary = TRUE)

# mcmc setup
mcmc<-mcmc.control(chains=3,warmup=1000,iter=7000,seed=53453,
                   propInit=0.5,adapt_delta = 0.8)

# Now fit emax model

fitbin <- fitEmaxB(y=y, dose=dose, count=cnt, prior=prior_bin, mcmc=mcmc, binary=TRUE)

summary(fitbin)
## Inference for Stan model: mrmod.
## 3 chains, each with iter=7000; warmup=1000; thin=1; 
## post-warmup draws per chain=6000, total post-warmup draws=18000.
## 
##              mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
## led50        2.54    0.01 0.81    1.49    2.03    2.38    2.82    4.73  3510    1
## lambda       1.19    0.01 0.48    0.49    0.86    1.11    1.42    2.39  6923    1
## emax         2.64    0.02 0.91    1.70    2.17    2.45    2.83    5.00  3202    1
## e0[1]       -0.96    0.00 0.25   -1.44   -1.12   -0.97   -0.81   -0.46  7919    1
## difTarget    2.22    0.00 0.31    1.61    2.01    2.22    2.42    2.82  7972    1
## loglambda    0.10    0.01 0.39   -0.72   -0.15    0.10    0.35    0.87  5798    1
## lp__      -215.14    0.03 1.74 -219.59 -215.96 -214.73 -213.88 -213.02  3944    1
## 
## Samples were drawn using NUTS(diag_e) at Sun Sep 20 20:22:21 2026.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

A plot the fitted dose response model:

plot(fitbin)
plot of chunk plotEmaxBayes
plot of chunk plotEmaxBayes

Continuous Response

Deriving a Mixture Prior Using RBesT

This example is from Application of RBesT in a Normal Endpoint. Additional details for the mixture prior are in this primary reference. The computations are very similar to the binary example.

Data

The primary endpoint is the change from baseline in Crohn’s Disease Activity Index (CDAI), which is assumed to be normally distributed. Note that for CDAI, an improved outcome corresponds to a negative change from baseline. The estimated standard deviation from historical studies is \(\sigma=88\).

h.data <- crohn
crohn_sigma <- 88
h.data$y.se <- crohn_sigma / sqrt(h.data$n)
kable(h.data, digits=2, caption = "CDAI from Published Studies")%>%
  kable_styling(full_width = FALSE, position = "center")
CDAI from Published Studies
study n y y.se
Gastr06 74 -51 10.23
AIMed07 166 -49 6.83
NEJM07 328 -36 4.86
Gastr01a 20 -47 19.68
APhTh04 25 -90 17.60
Gastr01b 58 -54 11.55

Computing a MAP Prior MCMC Sample using gMAP

An MCMC sample from the Meta Analytic Predictive (MAP) prior can be computed using the gMAP() function in RBesT. The between trial heterogeneity is a Half-Normal distribution with scale parameter \(\tau\).
A conservative specification is \(\tau= \sigma/2\). A normal prior distribution is specified for the mean placebo change from baseline with a prior mean of \(0.0\) and a unit-information prior SD.

map_mcmc <- gMAP(cbind(y, y.se) ~ 1 | study,
  weights = n, data = h.data,
  family = gaussian,
  beta.prior = cbind(0, crohn_sigma),
  tau.dist = "HalfNormal", tau.prior = cbind(0, crohn_sigma / 2)
)
print(map_mcmc)
## Generalized Meta Analytic Predictive Prior Analysis
## 
## Call:  gMAP(formula = cbind(y, y.se) ~ 1 | study, family = gaussian, 
##     data = h.data, 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.3   12.3 9.71 1.51 12.3  39.2
## 
## MAP Prior MCMC sample
##                  mean median   sd  q2.5   q50 q97.5
## theta_resp_pred -49.6  -48.4 18.8 -91.1 -48.4 -12.7

Parametric Approximation derived from the MCMC Sample

As with the binary example, a parametric approximation is computed using mixfit() with 2 components selected manually.

map <- mixfit(map_mcmc, Nc=2)
print(map)
## EM for Normal Mixture Model
## Log-Likelihood = -16937.04
## 
## Univariate normal mixture
## Reference scale: 88
## Mixture Components:
##   comp1       comp2      
## w   0.7405834   0.2594166
## m -48.3825486 -53.0860498
## s  10.6234701  31.9008848
## check accuracy of mixture fit
plot(map)$mix
plot of chunk mixmapc
plot of chunk mixmapc

Robustification and Effective Sample Size

The robustification of the mixture prior is unchanged from the binary example:

map_robust <- robustify(map, weight = 0.5, mean = -50)
print(map_robust)
## Univariate normal mixture
## Reference scale: 88
## Mixture Components:
##   comp1       comp2       robust     
## w   0.3702917   0.1297083   0.5000000
## m -48.3825486 -53.0860498 -50.0000000
## s  10.6234701  31.9008848  88.0000000

The effective sample size (ESS) also computed using ess() function:

round(ess(map, method = 'elir'))
## [1] 40

Bayesian Emax Model with a Mixture Prior for E0

Simulated dose response study

The methods are illustrated with simulated data from a dose response study with placebo and four doses ranging from 5 mg to 100 mg. The simulated responses are normally distributed with means determined by a 3-parameter Emax model with \(E_0=-50\), a difference with placebo at the \(100\) mg dose of \(difTarget= -1SD\), with the \(SD=75\), and the \(ED_{50}=30\)mg. The simulated data are consistent with the mixture prior.

The reduced sample size in the placebo group matches the sample size gain from the informative mixture prior.

# Dose groups
doselev<-c(0, 5, 25, 50, 100) 
# Sample size
ss<-c(40, 80, 80, 80, 80)
## Population Parameters
e0<--50
ed50<- 30
dtarget<- 100
diftarget<- -75
emax<-solveEmax(diftarget,dtarget,log(ed50),1,e0)
sdy<-75
pop<-c(log(ed50),emax,e0)

# Mean for each dose level
mu<-emaxfun(doselev,pop)

# Generate single study data
dose_data <- rep(doselev, times=ss)
mu_data <- rep(mu, times=ss)
set.seed(20260217)
y <- rnorm(n=sum(ss), mean=mu_data, sd=sdy)

study.data <- data.frame(y=y, dose= dose_data)

# Create summary data
summ.data <- study.data %>% group_by(dose) %>%
  summarise(
    n=n(),
    mean = mean(y),
    std=sd(y)
  )

kable(summ.data, digits=2, caption="Summary at each Dose Level") %>%
  kable_styling(full_width = FALSE, position = "center")
Summary at each Dose Level
dose n mean std
0 40 -62.38 78.34
5 80 -67.24 79.20
25 80 -106.93 75.80
50 80 -111.03 84.94
100 80 -120.21 74.33

Model fit using fitEmaxB

To set up a mixture prior for E0, we use emaxPrior.control() similar to before without specifying binary=TRUE:

prior.emax <- emaxPrior.control(epmu=0,epsca=100,mixP=3, 
                                mu_ep=c(-48.6, -53.3, -50),sd_ep = c(10.9,32.7, 88.0),
                                w_ep = c(0.35, 0.15, 0.5), difTargetmu = -150, difTargetsca =75,
                                dTarget = 100, p50=30, sigmalow =10, sigmaup = 500)

msSat <- sum((summ.data$n-1)*(summ.data$std)^2)/(sum(summ.data$n)-length(summ.data$n))
mcmc <- mcmc.control(chains=3)

## Fit Emax Model

fitout<-fitEmaxB(summ.data$mean,summ.data$dose,prior.emax,modType=4,prot=rep(1, nrow(summ.data)),
                 count=summ.data$n,msSat=msSat,mcmc=mcmc)
summary(fitout)
## Inference for Stan model: mrmod.
## 3 chains, each with iter=4333; warmup=1000; thin=1; 
## post-warmup draws per chain=3333, total post-warmup draws=9999.
## 
##               mean se_mean    sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
## led50         3.30    0.02  1.18     1.75     2.54     3.04     3.77     6.42  2685    1
## lambda        1.07    0.01  0.46     0.45     0.76     0.99     1.29     2.22  4990    1
## emax       -101.04    1.50 67.98  -266.13  -106.73   -84.20   -70.04   -50.45  2051    1
## e0[1]       -54.53    0.11  8.41   -71.39   -60.15   -54.50   -48.98   -38.11  5812    1
## sigma[1]     78.80    0.04  2.97    73.24    76.75    78.71    80.79    84.71  7011    1
## difTarget   -67.61    0.15 11.13   -89.23   -75.20   -67.66   -60.38   -45.49  5428    1
## loglambda    -0.01    0.01  0.40    -0.80    -0.28    -0.01     0.25     0.80  4790    1
## lp__      -1863.67    0.03  1.84 -1868.27 -1864.60 -1863.26 -1862.31 -1861.27  3026    1
## 
## Samples were drawn using NUTS(diag_e) at Sun Sep 20 20:22:31 2026.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Finally, the plot function displays the fitted Emax curve

plot(fitout)
plot of chunk plotc
plot of chunk plotc

References

[1] Baeten D. et al., The Lancet, 2013, (382), 9906, p 1705

[2] Neuenschwander B et. al, Clin Trials. 2010; 7(1):5-18

[3] Kass RE, Wasserman L, J Amer Statist Assoc; 1995, 90(431):928-934.

[4] Thomas N et. al, Statistics in Biopharmaceutical Research, 2014, 302-317

[5] Thomas N et. al, Statistics in Biopharmaceutical Research, 2016, 302-317

[6] Wu, J.et. al., Statistical Methods in Medical Research, 2017