Data

We extract the data for compound taid=23 and Protocol=A3921019 from the metaData repository which contains hundreds of compounds/protocols.

# Access metaData

data("metaData")

# Extract study data for taid=23 and protno='A3921019'

dat1019 <- metaData[metaData$taid==23 & metaData$protno=="A3921019",]

# Extract necessary columns

dat1019 <- dat1019 %>% select(c("protno","dose", "rslt","se","sampsize","indication", "endpointShort", "primtype", "regimen",))

kable(dat1019, row.names=FALSE, align="c") %>% 
  kable_styling(full_width=FALSE,
                bootstrap_options = c("hover"),
                position = "center")

Key study design components:

• Placebo and 3 actives BID doses - 10mg , 30mg and 60mg in total daily dose
• Parallel arm design
• Primary endpoint is ACR20 at week 6. ACR20 is a binary endpoint (Responder / Non Responder)
• Targeted efficacy: A 0.3 point improvement over placebo is commercially meaningful.
• An estimated improvement of 0.2 or lower versus placebo (even if p-value is significant) is considered failure for Phase 3

A quick visualization of the dose response shape uses clinDR:plotD:

attach(dat1019)

plotD(rslt, dose, meansOnly = TRUE, sem=se, ylab="ACR20 Proportion")
#> $ggp

plot of chunk plotraw

#> 
#> $means
#> [1] 0.3387097 0.7241379 0.7941176 0.8095238
#> 
#> $se
#> [1] 0.06010548 0.05868710 0.04903402 0.04947262

Pre-specified alpha-controlled hypothesis testing: MCP-Mod Trend Test

Pre-specified alpha-controlled tests are not required but they can be useful as reviewers often expect/require p-value reporting. Once a successful p-value is reported, the primary objectives of dose-response analyses can be pursued. Obtaining a successful p-value is typically a minimal requirement. A signal so weak it cannot be differentiated from noise will not contribute much to the differentiation of active doses for dose selection.

Selecting MCP-Mod contrasts created from Emax shapes is a natural first step in creating a prior distribution for Bayesian Emax modeling. A detailed overview of the implementation of MCP-Mod is available in the package vignettes MCP-Mod. Some default values values for candidate model selection that have been found to be useful are:

• Power (Hill) Parameter = 1
• ED50: half- of the lowest active dose , median dose, and average of the two highest doses. If there are only 3 active doses, then the contrast based on the median dose is omitted.

The candidate dose response models for MCP-Mod in the example are:

Ndose <- length(dose)
#ED50 candidates
ed50_cnd <- c((dose[1] + dose[2])/2, (dose[Ndose-1]+dose[Ndose])/2)
#Lambda = 1
lambda_cnd <- rep(1,2)

parms_cnd <- cbind(ed50_cnd, lambda_cnd)

# Candidate Models
testMods <- Mods(sigEmax=parms_cnd, doses=dose, placEff = 0, maxEff = 1)

plot(testMods)

plot of chunk mcpmod_candidate

MCP-Mod handles the binary data through the asymptotic normal approximation of logit transformed rates:

\[Var(logit(p)) = \frac{1}{n*p*(1-p)}\]

# Logit transformed rates
ylogit <- qlogis(rslt)
# Asymptotic variance covariance matrix
VARCOV <- diag(1/(rslt*(1-rslt)*sampsize))
#Optimal Contrasts
contMat <- optContr(testMods, S=VARCOV)

MCTtest(dose, ylogit, contMat = contMat, S=VARCOV, type='general')$tStat
#> sigEmax1 sigEmax2 
#> 6.069198 5.415233 
#> attr(,"pVal")
#> [1] 1.170114e-09 5.422183e-08

The p-values are already alpha adjusted for the multiple candidate dose response shapes. The overall p-value is the lowest one. Clearly, there is a strong dose response signal.

Setting up Prior for DR analysis

The prior distribution for the Emax model parameters is specified in package clinDR using the emaxPrior.control() function. For binary data like the example, the proportion of responders is modeled using on the logit scale, so the prior for the Emax parameters is specified on the logit scale. The clinDR package allows us to specify priors for the Emax model parameters using the emaxPrior.control() function. The response parameter is modeled using the logit link function, so the prior for the Emax parameter is specified on the logit scale.

• \(E_0\) (placebo response) is expected to be around 0.15, so we set the prior mean for \(E_0\) to logit(0.15). The historical placebo response is highly variable and changing over time, so we use a diffuse prior around the mean

\[E_0 \sim t_5(epmu=logit(0.15),epsca= 4)\]

• \(difTarget\) (placebo adjusted effect at the maximum dose on the logit scale) We use a diffuse prior around 0, i.e., no effect, which is common in early dose response studies

• \(ED_{50}\) (dose at which 50% of the maximum effect is achieved) Based on pre-clinical data and pharmacology models, the \(p_{50}\), the preliminary prediction of the \(ED_{50}\) is 10 mg

For the other parameters, default priors included in emaxPrior.control() are used.

A description of prior and model is available through clinDR:print

# Prior control object
# defaults
# effDF=5,parmDF=5,             
# loged50mu=0.0,loged50sca=1.73,
# loglammu=0.0,loglamsca=0.425,parmCor=-0.45,
# lowled50=log(0.001),highled50=log(1000),
# lowllam=log(0.3),highllam=log(4.0)

prior <- emaxPrior.control(epmu=qlogis(0.15), epsca=4, dTarget=max(dose),
                           difTargetmu=0, difTargetsca=4, p50=10,
                           binary=TRUE)


print(prior, doc=TRUE, diffuse=c(pbo=TRUE,eff=TRUE), 
      docType='sap')
#> [1] "sap/binary/diffuse/t"

Fit Bayesian Dose Response

The Bayesian Emax model is fitted using clinDR:fitEmaxB().

• \(y:\) Response data in 1/0 format
• \(dose:\) Dose levels, needs to be in the same order as response data
• \(count:\) Sample size of responder (i.e 1) and non-responder (i.e 0), needs to be in the same order as dose

### convert proportions to 1/0 count data
ndose<-length(dose)
ybin<-c(rep(1,ndose),rep(0,ndose))
nbin<-c(round(sampsize*rslt),sampsize-round(sampsize*rslt))
dbin<-c(dose,dose)

cat(
  "Response data in 1/0 format: ", ybin, "\n",
      "Dose levels: ", dbin, "\n",
    "Count of responders and non-responders: ", nbin, "\n"
)
#> Response data in 1/0 format:  1 1 1 1 0 0 0 0 
#>  Dose levels:  0 10 30 60 0 10 30 60 
#>  Count of responders and non-responders:  21 42 54 51 41 16 14 12

When patient level data are available, the 0/1 value and dose of each patient can be entered with a sample size of 1 for each 0 or 1 value. In addition to study data:

• \(prior:\) Prior control object created using emaxPrior.control()
• \(mcmc:\) MCMC control parameters. There are default values for all required inputs. Consideration should be given to modifying these settings (e.g. warmup, iterations, chains, thin etc)

mcmc<-mcmc.control(chains=3,thin=3,iter=3333*3+1000, warmup=1000)

fitout<-fitEmaxB(ybin,dbin,prior,count=nbin, mcmc=mcmc,binary=TRUE)

A modified maximum likelihood fit can also be computed using clinDR:fitEmax. This estimation is not illustrated in the current example.

Posterior Diagnostics

Convergence of MCMC should always be checked. The rstan objects are stored under fitEmaxB object

names(fitout)
#>  [1] "estanfit"  "y"         "dose"      "prot"      "count"     "nbase"     "xbase"     "dimFit"    "vcest"    
#> [10] "modType"   "binary"    "pboAdj"    "msSat"     "prior"     "mcmc"      "localParm"

Traceplot: This helps in assessing the convergence of the chains.

traceplot(fitout$estanfit,pars=c('led50','loglambda','e0[1]','difTarget'))

plot of chunk traceplot

pair plot: Pair plot visualizes posterior marginal distributions and correlation between parameters. It also helps diagnosing mixing of chains, identifiability issues, etc.

pairs(fitout$estanfit,pars=c('led50','loglambda','e0[1]','difTarget'))

plot of chunk pairplot

Finally, clinDR provides a convenient function to check the goodness of fit of the model. The function bpchkMonoEmax() checks for monotonicity of the Emax model fit comparing the best response from lower doses to the response from highest dose. The function returns a posterior predictive check designed to be sensitive to non-monotonicity. A small predictive probability indicates that the model fit is inconsistent with a monotonic Emax dose response relationship.

bpchkMonoEmax(fitout)
#> [1] 0.4922492

Posterior Inference

Posterior summaries of the Emax parameters are displayed by:

summary(fitout,pars=c('led50','lambda','emax','e0','difTarget','loglambda'))
#> Inference for Stan model: mrmod.
#> 3 chains, each with iter=10999; warmup=1000; thin=3; 
#> 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        1.55    0.02 1.54   -1.72    0.80    1.58    2.28    4.94  7406    1
#> lambda       1.08    0.01 0.54    0.37    0.70    0.98    1.32    2.48  7889    1
#> emax         2.67    0.02 1.36    1.50    2.02    2.37    2.88    5.77  6886    1
#> e0[1]       -0.64    0.00 0.27   -1.19   -0.82   -0.64   -0.47   -0.13  7682    1
#> difTarget    2.13    0.00 0.37    1.42    1.88    2.13    2.37    2.88  7707    1
#> loglambda   -0.04    0.01 0.48   -1.00   -0.35   -0.02    0.28    0.91  7998    1
#> lp__      -140.27    0.02 1.87 -144.99 -141.18 -139.85 -138.88 -137.91  5859    1
#> 
#> Samples were drawn using NUTS(diag_e) at Sun Sep 20 20:21:39 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).

Although they are easy to summarize and extract for further calculations, the draws from the posterior distribution of the model parameters should be viewed with caution and skepticism because they are usually poorly identified by the data in typical dose response studies. The parameter estimates are most often highly correlated and unstable. The resulting estimated dose response within the observed dosing range, however, is much better determined. The plot() function can be utilized to visualize the fitted dose-response curve along with the observed data. The black bars represents posterior credible interval, while the grey ones represent posterior predictive intervals for future samples rates in a similar study. The observed data points are marked using red asterisks

plot(fitout,xlab='Total Daily Dose',ylab='ACR20', clev=0.9)

plot of chunk plot_fit

The placebo adjusted dose response curve is plotted by setting plotDif=TRUE.

plot(fitout,xlab='Total Daily Dose',ylab='Placebo Adjusted ACR20', 
     clev=0.9, dref=0, plotDif=TRUE)

plot of chunk plot_fit_dif

Assessing the Need For a Second Study and Its Potential Design

The evaluated dose range in the first study spanned only a six‑fold difference between the lowest and highest doses, prompting a key question:

• Are all active doses on a plateau?

This is a critical concern given the compound’s immunosuppressive mechanism and intended chronic usage, particularly in the presence of dose‑related safety risks at higher doses. A strong and defensible dose selection rationale is therefore essential:

• How far can the dose be reduced while maintaining clinically meaningful efficacy?

To assess the location of the plateau, the lowest active dose is directly compared with the middle dose, evaluating whether additional exposure provides incremental benefit.

# Posterior 
pred10v30 <- predict(fitout,dosevec=c(10,30),dref=30,clev=0.90)
names(pred10v30)
#>  [1] "pred"      "predMed"   "lb"        "ub"        "se"        "fitdif"    "fitdifMed" "lbdif"     "ubdif"    
#> [10] "sedif"     "simResp"   "sigsim"

A 90% CI of difference between 10 mg vs 30 mg

dif10v30 <-round(pred10v30$fitdif[1],3)
lb10v30 <- round(pred10v30$lbdif[1],3)
ub10v30 <- round(pred10v30$ubdif[1],3)

cat("Difference in response between 10mg and 30mg: ", dif10v30, "\n",
    "90% CI: [", lb10v30, ", ", ub10v30, "]\n")
#> Difference in response between 10mg and 30mg:  -0.071 
#>  90% CI: [ -0.147 ,  -0.002 ]

The 90% credible interval confirms that the response of the 10 mg dose may be close to the plateau. The first study failed to convincingly identify the dosing range where the dose-response is rapidly changing. The 10 mg dose was chosen as the lowest dose for the first study because only 5 and 10 mg tablets were manufactured (Recall the dosing is BID). While the first study was accruing, a 1 mg tablet was manufactured so it is now possible to study much lower doses. A practical dosing constraint based on safety assessments from the first study is that the doses above 30 mg cannot be included in future studies.

A key design question is whether a 2mg dose will display lower efficacy that provides a lower bound on useful doses. If not, a third study will be necessary. Should we proceed with the 2mg low dose or wait for a lower dose to be manufactured?

Predicting the response to a 2mg based based on current study data:

# Simulate 2 mg and 30 mg response based on current data
pred2v30 <- predict(fitout,dosevec=c(2,30),dref=30,clev=0.90)

dif2v30 <-round(pred2v30$fitdif[1],3)
lb2v30 <- round(pred2v30$lbdif[1],3)
ub2v30 <- round(pred2v30$ubdif[1],3)

cat("Difference in response between 10mg and 30mg: ", dif2v30, "\n",
    "90% CI: [", lb2v30, ", ", ub2v30, "]\n")
#> Difference in response between 10mg and 30mg:  -0.238 
#>  90% CI: [ -0.403 ,  -0.031 ]

While not certain, it is likely the 2mg dose will have appreciably less efficacy than the highest 30 mg dose. However, reviewers focus more on sample proportions rather than theoretical population rates estimated from a fitted model. To account for the naive assessment likely applied to the second study results by key decision makers, we will compute the following criteria for demonstrating the lowest dose lacks sufficient efficacy:

• C1: What is the chance that the observed sample proportion for 2mg is 0.3 point better improvement compared to placebo
• C2: What is the chance that the observed sample rate for 2mg is 0.1 point (remember 0.2 point is the minimal difference) or less worse than the observed rate for 30 mg

The future study will include approximately 50 patients per dose group, based largely on practical operational considerations:

### obtain simulated pbo response
pred2v0<-predict(fitout,dosevec=c(2.0,0),dref=0)

nsim<-length(pred2v30$simResp[,1])
nnew<-50 # Sample size in New study

ynew2<-rbinom(nsim,nnew,pred2v30$simResp[,1])   ## 2 mg dose new study
ynew30<-rbinom(nsim,nnew,pred2v30$simResp[,2])  ## 30 mg highest new study
ynew0<-rbinom(nsim,nnew,pred2v0$simResp[,2])  ## placebo in new study

# Estimated sample proportion
pnew2<-ynew2/nnew
pnew30<-ynew30/nnew
pnew0<-ynew0/nnew

# Both C1 and C2 satisfied
pjoint<-mean(pnew30-pnew2<0.1 & 
                            pnew2-pnew0>0.3)
pjoint<-round(pjoint,2)

pjoint
#> [1] 0.12

detach(dat1019)

The refined joint assessment predicting results for observed rates in a future study show that it is unlikely (though not certain) that the 2mg dose will display adequate efficacy similar to the 30 mg dose. Based on PK-PD simulation and these results, a new study was proposed with 0mg, 2mg, 6mg, 10mg, 20mg and 30mg doses.

Bayesian Emax Model fitting

The Bayesian Emax fitting is the same as before, except we fit both the studies together by adding variable indicating the protocol for each dose/response input. The prior distribution is unchanged except there are now 2 independent prior distributions for the 2 placebo responses.

attach(dat23)

ndose<-length(dose)
ybin<-c(rep(1,ndose),rep(0,ndose))
nbin<-c(round(sampsize*rslt),sampsize-round(sampsize*rslt))
dbin<-c(dose,dose)
# protocol indicator
prot<-c(as.character(protno),as.character(protno))

# fitEmaxB
fitout2<-fitEmaxB(ybin,dbin,prior,count=nbin,prot=prot,
                                 mcmc=mcmc,binary=TRUE)

Posterior Diagnostics

Note: e0[1] and e0[2] represent the placebo response for the two studies, while the other parameters are shared across studies.

Traceplot

traceplot(fitout2$estanfit,pars=c('led50','loglambda','e0','difTarget'))

plot of chunk study2_traceplot

Pairplot

pairs(fitout2$estanfit,pars=c('led50','loglambda','e0','difTarget'))

plot of chunk study2_pairplot

Posterior Inference

A downward shift in response rates across dose groups was observed in the second trial and can be evaluated through the model parameter summaries. After accounting for this shift, treatment effects relative to placebo remain consistent.

This pattern has been commonly observed in dose-response meta‑analyses, including studies with laboratory endpoints, when multiple trials are combined.

summary(fitout2,pars=c('led50','lambda','emax','e0','difTarget','loglambda'))
#> Inference for Stan model: mrmod.
#> 3 chains, each with iter=10999; warmup=1000; thin=3; 
#> 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        2.21    0.01 0.84    1.06    1.72    2.04    2.48    4.52  5790    1
#> lambda       1.10    0.00 0.42    0.48    0.81    1.04    1.31    2.14  7538    1
#> emax         2.82    0.02 1.12    1.80    2.24    2.56    3.03    5.77  5198    1
#> e0[1]       -0.66    0.00 0.22   -1.09   -0.80   -0.65   -0.51   -0.24  8650    1
#> e0[2]       -0.90    0.00 0.22   -1.35   -1.05   -0.90   -0.75   -0.47  8761    1
#> difTarget    2.27    0.00 0.30    1.71    2.07    2.26    2.47    2.89  7679    1
#> loglambda    0.03    0.00 0.37   -0.73   -0.21    0.04    0.27    0.76  7164    1
#> lp__      -319.86    0.02 1.83 -324.41 -320.79 -319.46 -318.51 -317.46  5876    1
#> 
#> Samples were drawn using NUTS(diag_e) at Sun Sep 20 20:21:48 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).

The int = 0 option in the plotting function displays model fits for both protocols side by side, while int = 1 or int = 2 restricts the plot to the first or second protocol, respectively.

plot(fitout2,xlab='Total Daily Dose',ylab='ACR20', clev=0.9, int=0)

plot of chunk plot_study2

Dose Selection

The clinical team has clear dose selection criteria:

• A 0.3 point improvement versus placebo is the commercial target efficacy
• An improvement of <=0.2 point versus placebo is considered failure for Phase 3 even if its p-value is significant
• Safety is assessed separately and overlaid with efficacy to form a therapeutic index. A major challenge is that the safety issues are mostly observed only with long term dosing

The probability that the population response rates achieved the clinical criteria are computed over a range of potential doses based on the analysis of the two combined studies:

dvec<-c(2*(1:10),30)

# Simulate response for dose range dvec
pout2<-predict(fitout2,dvec)$simResp ### by default, pbo rate from first study

# Simulate placebo response
pbo2<-as.vector(predict(fitout2,0)$simResp)

# Median difference between dose and placebo
med23<-round(apply(pout2-pbo2,2,median),2)
# At least 0.2 improvement over placebo
p20<-round(apply(pout2-pbo2>0.2,2,mean),2)
# At least 0.3improvement over placebo
p30<-round(apply(pout2-pbo2>0.3,2,mean),2)

pop<-cbind(dvec,med23,p20,p30)
pop
#>    dvec med23  p20  p30
#> 2     2  0.12 0.10 0.00
#> 4     4  0.21 0.57 0.07
#> 6     6  0.27 0.89 0.31
#> 8     8  0.32 0.98 0.61
#> 10   10  0.35 1.00 0.80
#> 12   12  0.37 1.00 0.91
#> 14   14  0.39 1.00 0.95
#> 16   16  0.40 1.00 0.97
#> 18   18  0.41 1.00 0.98
#> 20   20  0.42 1.00 0.99
#> 30   30  0.46 1.00 1.00

The 10 mg dose meets the criteria for dose selection.

This assessment is for population rates. The Phase 3 observed rates will be used for the final regulatory and commercial decisions. These observed rates are more uncertain and the rate differences will be attenuated by missing data. At the time of the planned analyses, regarding dropouts as failures was the accepted approach for missing data adjustments, so is the approach assessed here using simulation generated from MCMC parameter draws from the fitted combined model. The Phase 3 design will be 2:1 with approximately 200 patients on drug, 100 on PBO. The team anticipated a 15% dropout rate. The probability of achieving the 2 clinical criteria over a range of doses:

nsim<-length(pbo2)
np<-ncol(pout2)
med23obs<-numeric(np)
p20obs<-numeric(np)
p30obs<-numeric(np)

### placebo simulations
y0sim<-rbinom(nsim*100,1,pbo2)

# Impute Missing data as non-response
y0sim[sample(1:(100*nsim),0.15*100*nsim)]<-0
y0sim<-matrix(y0sim,ncol=100)
p0sim<-apply(y0sim,1,sum)/100   

### doses
for(i in 1:np){
    ## simulate 200 patients for each mcmc prob for dose i
    ydsim<-rbinom(nsim*200,1,pout2[,i])
    ## random dropouts assigned non-response
    ydsim[rbinom(nsim*200,1,.15)==1]<-0
    ## row has 9999 obs for each mcmc-generated prob for dose i
    ydsim<-matrix(ydsim,ncol=200)
    pdsim<-apply(ydsim,1,sum)/200
    
    med23obs[i]<-round(median(pdsim-p0sim),2)
    p20obs[i]<-round(mean(pdsim-p0sim>0.2),2)
    p30obs[i]<-round(mean(pdsim-p0sim>0.3),2)
}
cbind(pop,med23obs,p20obs,p30obs)
#>    dvec med23  p20  p30 med23obs p20obs p30obs
#> 2     2  0.12 0.10 0.00     0.11   0.10   0.01
#> 4     4  0.21 0.57 0.07     0.18   0.38   0.06
#> 6     6  0.27 0.89 0.31     0.23   0.64   0.18
#> 8     8  0.32 0.98 0.61     0.27   0.81   0.33
#> 10   10  0.35 1.00 0.80     0.29   0.90   0.46
#> 12   12  0.37 1.00 0.91     0.32   0.94   0.58
#> 14   14  0.39 1.00 0.95     0.33   0.96   0.66
#> 16   16  0.40 1.00 0.97     0.34   0.97   0.72
#> 18   18  0.41 1.00 0.98     0.35   0.98   0.77
#> 20   20  0.42 1.00 0.99     0.36   0.98   0.80
#> 30   30  0.46 1.00 1.00     0.38   0.99   0.88

The 10 mg total daily BID dose is likely to provide clinical benefit, but there is considerable risk it will not achieve the commercial target. Increasing the dose to 20 mg produced much higher confidence of achieving the targeted efficacy but it has more risk of long term safety concerns.

To manage this uncertainty, both doses were advanced into Phase 3.
The modeling framework ultimately proved highly predictive of outcomes. The efficacy predictions were confirmed. The 10 mg dose was approved by most regulatory agencies, while the 20 mg dose was not approved due to longer‑term safety concerns. The analyses demonstrate the importance of integrating predictive modeling, missing data considerations, and safety risk into dose selection decisions.