Getting Started

Install StudyPrior package

devtools::install_github('igrave/StudyPrior/StudyPrior')

or install manually the dependencies and the package.

for(dep in c('cubature','parallel','optimr','optimx','VGAM','BB','mvtnorm')) { if (!require(dep, character.only=TRUE)) install.packages(dep)}

https://igrave.github.io/StudyPrior/packages/StudyPrior_0.0.15.zip (for Windows)

https://igrave.github.io/StudyPrior/packages/StudyPrior_0.0.15.tar.gz (source package)

Install INLA

install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)

A separate file with the R code from this Workshop can be found at https://igrave.github.io/StudyPrior/Workshop/code.R.

To see vignette on using StudyPrior package:

??StudyPrior

or online here.

Practical 1: Analysis of clinical trials using single historical control dataset

Data Description: Taken from Neuenschwander 2010 clinical trials paper, designing a Phase IV trial comparing a control treatment and a new intensified treatment of the same compound. The primary outcome is treatment failure (binary). There are 11 internal historical trials of essentially the same trial design.

library(StudyPrior)
library(INLA)


x.hist <- c(6,8,17,28,26,8,22,8,6,16,53)
n.hist <- c(33,45,74,103,140,49,83,59,22,109,213)
p.hist <- x.hist/n.hist

knitr::kable(data.frame(Study = 1:11,
                        Failures = paste0(x.hist,'/',n.hist),
                        Percent = paste0(round(p.hist*100),'%')),
             col.names = c("Study","Failures/Controls","Failure %"),
             align='c')

A new study takes place and there are 12 failures out of 80 controls (15%) and 8 failures out of 80 people on the new intensified treatment (10%).

x.c <- 12
n.c <- 80

x.t <- 8
n.t <- 80
  1. Take one historical study (different people should choose different studies, so we don’t all do the same one!) and calculate the probability that the new treatment is better than the “control” less intensive treatment (i.e. new treatment has lower failure rate) using each of the methods below. Complete the table provided in the Word document with the results from each analysis. You may also find it useful to produce plots of the posteriors for the control arm rate under each method.
my.choice <- 3    # Enter a study number between 1 and 11

x.h <- x.hist[[my.choice]]
n.h <- n.hist[[my.choice]]
  1. Only using the new study data (i.e. no historical data) (Hint: use simulation or numerical integration to compute Pr(treatment failure < control failure); see Lecture 2 slide 5)
#approximate with a Riemann sum
p <- seq(0,1,len=2000)
mean(dbeta(p, 1+x.c, 1+n.c-x.c) *
        pbeta(p, 1+x.t, 1+n.t-x.t, lower.tail = TRUE))

# numerical integration
integrate(function(p) dbeta(p, 1+x.c, 1+n.c-x.c) *
            pbeta(p, 1+x.t, 1+n.t-x.t, lower.tail = TRUE),
          lower=0, upper=1)


#Using the package StudyPrior
q1.nohist <- create.mixture.prior("beta",pars=matrix(c(1,1),ncol=2))
q1.nohist.post <- posterior.mixture.prior(x.c, n.c, mix=q1.nohist)

1-prob.control.smaller(xt=x.t, nt=n.t, posterior=q1.nohist.post)

#Plot posteriors for control and active rate
plot(q1.nohist.post)
q1.active <- create.mixture.prior("beta",pars=matrix(c(1,1),ncol=2))
q1.active.post <- posterior.mixture.prior(x.t, n.t, mix=q1.active)
plot(q1.active.post, col=2, add=TRUE)
  1. Pooling historical and new trial data (Hint: use simulation or numerical integration; see Lecture 2 slide 11-12)
q1.pool <- binom.PP.FIX(x=x.h, n=n.h, d=1, mix = TRUE)
q1.pool.post <- posterior.mixture.prior(x.c, n.c, mix=q1.pool)
1-prob.control.smaller(xt=x.t,nt=n.t,posterior=q1.pool.post)

plot(q1.pool.post)
  1. Using a power prior with fixed power parameters (Lecture 2 slide 14-15)
    • 0.8
    • 0.1
q1.pp08 <- binom.PP.FIX(x=x.h, n=n.h, d=0.8, mix = TRUE)
q1.pp08.post <- posterior.mixture.prior(x.c, n.c, mix=q1.pp08)

1-prob.control.smaller(xt=x.t,nt=n.t,posterior=q1.pp08.post)

q1.pp01 <- binom.PP.FIX(x=x.h, n=n.h, d=0.1, mix = TRUE)
q1.pp01.post <- posterior.mixture.prior(x.c, n.c, mix=q1.pp01)

1-prob.control.smaller(xt=x.t,nt=n.t,posterior=q1.pp01.post)

plot(q1.pp01.post)
  1. Estimating the posterior mean of the power prior and using this as a plug-in value (Hint: to compute the posterior mean of a0 you will need to select an appropriate prior for a0. You will then need to use either simulation to sample directly from the posterior for a0 (see Lecture 2 slide 20; use MCMC) or numerical integration to compute the integral on slide 24 of Lecture 2)
f.a0 <- function(a0,Nh,Rh,N,R,a,b) {
  logf <- {lgamma(a0*Nh+a+b) + lgamma(a0*Rh+R+a) + lgamma(a0*(Nh-Rh)+N-R+b)}-
  {lgamma(a0*Rh+a) + lgamma(a0*(Nh-Rh)+b) + lgamma(a0*Nh+N+a+b)}
  exp(logf)
}

k <- integrate(function(a0) f.a0(a0, n.h, x.h, n.c, x.c, a=1, b=1), lower=0,upper=1)

p <- seq(0,1,len=2000)
plot(p, f.a0(p, n.h, x.h, n.c, x.c, a=1, b=1)/k$value, type='l',
     xlab="a0",ylab="posterior density")

mean.a0 <- integrate(function(a0) a0 * f.a0(a0, n.h, x.h, n.c, x.c, a=1, b=1)/k$value,
                     lower=0,upper=1) 
mean.a0$value

q1.pp.mean <- binom.PP.FIX(x.h, n.h, d=mean.a0$value, mix=TRUE)

q1.pp.mean.post <- posterior.mixture.prior(x.c,n.c, mixture.prior = q1.pp.mean)
plot(q1.pp.mean.post)
1-prob.control.smaller(xt=x.t,nt=n.t,posterior=q1.pp.mean.post)
  1. Estimating the power prior jointly with the control failure rate (see Lecture 2 slide 18)
q1.joint <- binom.PP.FB(x.h, n.h,  d.prior.a=1, d.prior.b=1, mix = TRUE)
q1.joint.01 <- binom.PP.FB(x.h, n.h, d.prior.a=0.02, d.prior.b=0.02, mix = TRUE)

q1.joint.post <- posterior.mixture.prior(x.c,n.c, mixture.prior = q1.joint)
q1.joint.01.post <- posterior.mixture.prior(x.c,n.c, mixture.prior = q1.joint.01)

plot(q1.joint.post)
plot(q1.joint.01.post, add=TRUE, col=2)

1-prob.control.smaller(xt=x.t,nt=n.t,posterior=q1.joint.post)
1-prob.control.smaller(xt=x.t,nt=n.t,posterior=q1.joint.01.post)
plot(q1.nohist.post, ylim=c(0,14), xlim=c(0,0.4))
plot(q1.pool.post, add=TRUE, col=2)
plot(q1.pp08.post, add=TRUE, col=3)
plot(q1.pp01.post, add=TRUE, col=4)
plot(q1.pp.mean.post, add=TRUE, col=5)
plot(q1.joint.post, add=TRUE, col=6)
plot(q1.joint.01.post, add=TRUE, col=7)
legend("topright", legend=c("No hist", "Pooled", "PP FIX 0.8", "PP FIX 0.1", "PP MEAN", "PP FULL 1", "PP FULL 0.2"), lty=rep(1,7), col=1:7)
  1. If you have time, repeat some of the above with either a different historical trial or a different dataset for the new study.
my.choice <- 7

#Repeat as before

Practical 2: Analysis of clinical trial using multiple historical control datasets

Before you analyse your new trial using multiple historical control datasets, there are two additional priors that were covered in Lectures 3 (robust mixture prior) and 4 (empirical power prior) that we will first look at with a single historical control dataset.

Practical 1 extension:

Use the same historical study that you chose for Practical 1, and do the following:

  1. Add a robust component to the beta prior for the historical controls that you used for question 1b. (See Lecture 3 slide 22).
q1.pool.robust <- make.robust(q1.pool, 0.5)
plot(q1.pool)
plot(q1.pool.robust, add=TRUE, col=2)

q1.pool.robust.post <- posterior.mixture.prior(x.c,n.t, mix=q1.pool.robust)

plot(q1.pool.post, lty=2, add=TRUE)
plot(q1.pool.robust.post, add=TRUE, col=2, lty=2)

1-prob.control.smaller(xt=x.t,nt=n.t,posterior=q1.pool.post)
1-prob.control.smaller(xt=x.t,nt=n.t,posterior=q1.pool.robust.post)
  1. Fit a power prior to the single historical control dataset using the empirical Bayes approach to choose the parameters. (See Lecture 4 slides 7-8).
q1.EB <- binom.PP.EB(x.h, n.h,X = x.c, n.c, mix = TRUE)
attr(q1.EB,"powers")

q1.EB.post <- posterior.mixture.prior(x.c, n.c, mix=q1.EB)
1-prob.control.smaller(xt=x.t,nt=n.t,posterior=q1.EB.post)

Practical 2: Multiple historical studies

We will now look at how to extend the methods used in practical 1 to include multiple historical control datasets. For this exercise, you will use all 11 historical trials from Neuenschwander 2010 introduced in Practical 1.

  1. Do a meta analysis of all the studies to obtain a meta-analytic predictive (MAP) prior distribution.
q2.MAP <- binom.MAP.FB(x.hist, n.hist)
  1. Complete the following steps to fit a mixture of beta distributions to approximate the MAP prior:
  1. Fit a mixture of beta distributions to the MAP prior distribution with 1,2 and 3 components. Report the values of the fitted Beta parameters and the mixture weights for each fitted distribution.
MAP4 <- conj.approx(q2.MAP, type="beta", max.degree = 4)
MAP3 <- conj.approx(q2.MAP, type="beta", max.degree = 3)
MAP2 <- conj.approx(q2.MAP, type="beta", max.degree = 2)
MAP1 <- conj.approx(q2.MAP, type="beta", max.degree = 1)

p <- seq(0,1,len=500)
plot(p, q2.MAP(p), type='l')

plot(MAP1, col=2, add=TRUE)
plot(MAP2, col=3, add=TRUE)
plot(MAP3, col=4, add=TRUE)
plot(MAP4, col=5, add=TRUE)
  1. Compare the components to the fitted density using the Kullback-Leibler divergence criterion to decide upon the best fit
KLD <- function(prior){
  integrate(function(p) q2.MAP(p) *(log(q2.MAP(p)) - log(eval.mixture.prior(p,mixture=prior))),
            lower=0.1, upper=.9)
}

KLD(MAP1)
KLD(MAP2)
KLD(MAP3)
KLD(MAP4)
  1. With the same new study data as Practical 1 (i.e. 12 failures out of 80 controls and 8 failures out of 80 people on the new intensified treatment), do an analysis of the new trial data and all 11 historical control datasets using each of the following priors.
  1. the best MAP prior approximation
q2.MAP4.post <- posterior.mixture.prior(x.c,n.c, mix=MAP4)
  1. a robust version of the MAP prior built by adding a fourth component Beta(1,1) with
  • weight 0.1
  • weight 0.5
q2.MAP.robust.01 <- make.robust(MAP4, 0.1)
q2.MAP.robust.01.post <- posterior.mixture.prior(x.c,n.c, mix=q2.MAP.robust.01)

q2.MAP.robust.05 <- make.robust(MAP4, 0.5)
q2.MAP.robust.05.post <- posterior.mixture.prior(x.c,n.c, mix=q2.MAP.robust.05)
  1. Fit a power prior with fixed parameter for all studies (you can choose to fix the power parameter to the same value for all studies, or specify different values for each study)
q2.fix <- binom.PP.FIX(x.hist, n.hist, d=0.4 , mix = TRUE)
q2.fix.diff <- binom.PP.FIX(x.hist, n.hist, d=rep(c(0.8,0.4), times=c(5,6)), mix = TRUE)

q2.fix.post <- posterior.mixture.prior(x.c, n.c, mix=q2.fix)
q2.fix.diff.post <- posterior.mixture.prior(x.c, n.c, mix=q2.fix.diff)
  1. Fit a power prior using the joint approach (as in 1e in practical 1). Choose the parameters for the prior on the the weights.
q2.FB <- binom.PP.FB(x.hist, n.hist,
                     d.prior.a = 1, d.prior.b=1,
                     p.prior.a = 1, p.prior.b = 1,
                     mix = TRUE)

q2.FB.post <- posterior.mixture.prior(x.c, n.c, mix=q2.FB)
  1. Fit a power prior using the empirical Bayes approach to choose the parameters
q2.EB <- binom.PP.EB(x.hist, n.hist, x.c, n.c, mix = TRUE)
attr(q2.EB,"powers")

q2.EB.post <- posterior.mixture.prior(x.c, n.c, mix=q2.EB)
  • For each analysis, compute the posterior probability that the new treatment is better than control, and also record the prior and posterior weights for each mixture component (MAP and robust MAP priors), and the values of the power parameter for each study (fixed and empirical power priors) in the table provided in the Word document.
1-prob.control.smaller(x.t,n.t, posterior=q2.MAP4.post)
1-prob.control.smaller(x.t,n.t, posterior=q2.MAP.robust.01.post)
1-prob.control.smaller(x.t,n.t, posterior=q2.MAP.robust.05.post)

1-prob.control.smaller(x.t,n.t, posterior=q2.fix.post)
1-prob.control.smaller(x.t,n.t, posterior=q2.fix.diff.post)

1-prob.control.smaller(x.t,n.t, posterior=q2.FB.post)
1-prob.control.smaller(x.t,n.t, posterior=q2.EB.post)
  • Produce a plot to overlay the posterior distribution of the control rate for each prior.
plot(q2.MAP4.post, ylim=c(0,25), xlim=c(0,0.5), xlab=expression(theta))
plot(q2.MAP.robust.01.post, col=2, add=TRUE)
plot(q2.MAP.robust.05.post, col=3, add=TRUE)

plot(q2.fix.post, add=TRUE, col=4)
plot(q2.fix.diff.post, add=TRUE, col=5)

plot(q2.FB.post, add=TRUE, col=6)
plot(q2.EB.post, add=TRUE, col=7)

legend("topright", legend=c("MAP", "Robust MAP 0.1", "Robust MAP 0.5", "Fixed PP (1)",
                            "Fixed PP (2)", "Joint PP", "Empirical PP", 
                            "Historical control rates", "Current control rate"), lty=c(rep(1,7), NA, NA), col=c(1:7, 1,2), pch=c(rep(NA, 7), 16,16))
points(y=rep(0,length(x.hist)), x=x.hist/n.hist,  pch=16)
points(y=0, x=x.c/n.c, pch=16, col=2)
  1. If you have time, repeat some of the above with a different current study (e.g. 36 failures out of 320 controls and 32 failures out of 320 people on the new intensified treatment, or make up your own new study data)
# x.c <- 36
# n.c <- 320
# x.t <- 32
# n.t <- 320
#Continue as before

Practical 3: Designing at trial that includes historical controls

Use all the historical trials from Neuenschwander 2010 introduced in Practical 1. In this question, you will look at how the historical controls can be used to either reduce the number of current controls recruited in the new study, or increase the overall power.

  1. Use the MAP prior that you fitted to all the historical control trials in Practical 2 to help design a new study.
# q2.MAP
q3.map <- MAP4
  1. Use a standard sample size calculation (Lecture 5 slide 13) to calculate the sample size for a new study ignoring the historical data to have 80% power and 5% type 1 error rate, and assuming that the control treatment has a 22% failure chance and the experimental treatment has a 15% failure chance
power <- 0.8
type1e <- 0.05
p1 <- 0.22 #control failure
p2 <-  0.15 # experimental treatment failure

(qnorm(1-power)+qnorm(type1e/2))^2 * ((p1*(1-p1))+p2*(1-p2))/(p1-p2)^2

power.prop.test(n=NULL, p1,p2, type1e, power)
  1. Using the historical data, by how much could you reduce the sample size by? (Hint calculate the ESS of the MAP prior; see Lecture 5 slides 7-10)
q3.map.ess <- ess.morita.mixture.prior(q3.map,LEN = 500) #takes about 2 minutes

q3.map.ess
  1. Alternatively you could use the historical data to increase the power of the new study by recruiting the full sample size (as calculated in part a) and including the historical data via the MAP prior. What is the new power of this study? (Hint: This can be done via simulation - see Lecture 5 slide 15. In order to compute power for a Bayesian analysis, you should define the success criterion for a trial to be Pr(Treat < Control)>0.975).
N.new <- 482
#Calculate the posterior since we need it for the next questions
q3.map.post <- lapply(0:N.new, posterior.mixture.prior, ns=N.new, mixture.prior=q3.map)

# sig.matrix calculates the test [Pr(Control < Treat)>level] for all values of Xt, Xc
# so we need  1-0.975=level for Treat<Control and then the TRUEs <-> FALSEs
SIG <- !sig.matrix(n.control=N.new, n.treatment=N.new, level=0.025, posterior = q3.map.post)

# Calculate power to detect a treatment difference of -0.05
q3c.power <- calc.power(sig.mat = SIG, treatment.difference = 0.15-0.22, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)
# Plot power as a function of assumed control failure rate for new trial
plot(y=q3c.power, x=seq(0.1,0.5,len=41), xlab="Control Failure Rate (new trial)", ylab="Power", ty='b')
abline(v=0.22, lty=2)  # mark value of control rate assumed for new trial 
x=seq(0.1,0.5,len=41)
q3c.power[x==0.22]   # Power of study using MAP prior plus 482 current control subjects

# Calculate power for standard design with no historical data
q3.nohist.post <- lapply(0:N.new, posterior.mixture.prior, ns=N.new, mixture.prior=q1.nohist)

SIG.nohist <- !sig.matrix(n.control=N.new, n.treatment=N.new, level=0.025, posterior = q3.nohist.post)

q3c.power.nohist <- calc.power(sig.mat = SIG.nohist, treatment.difference = 0.15-0.22, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)

# Add power curve to plot
lines(y=q3c.power.nohist, x=seq(0.1,0.5,len=41), col=2)
  1. For the scenario in (c) above, compute the type 1 error for the trial design (Hint: use simulation as before, but generate replicate datasets under the null rather than alternative hypothesis).
# We can reuse the result of sig.matrix() since it doesn't depend on the treatment parameter or null/alternative hypothesis

#we use calc.power with treatment.difference=0 to find the type 1 error
q3c.t1e <- calc.power(sig.mat = SIG, treatment.difference = 0,
                        prob.range = c(0.1,0.5),
                        length=41, n.binom.control = N.new)
plot(y=q3c.t1e, x=seq(0.1,0.5,len=41), xlab="Control Failure Rate (new trial)", ylab="Type I Error", ty='b', ylim=c(0, 0.1))
abline(v=0.22, lty=2)  # mark value of control rate assumed for new trial 
x=seq(0.1,0.5,len=41)
q3c.t1e[x==0.22]   # Type 1 error of study using MAP prior plus 482 current control subjects

#Type 1 error for standard design without historical data
q3c.t1e.nohist <- calc.power(sig.mat = SIG.nohist, treatment.difference = 0,
                        prob.range = c(0.1,0.5),
                        length=41, n.binom.control = N.new)

#Add to plot
lines(y=q3c.t1e.nohist, x=seq(0.1,0.5,len=41), col=2)
  1. Repeat (c) - (d) in the previous question using
q3.rob.05 <- make.robust(q3.map, 0.5)
q3.rob.01 <- make.robust(q3.map, 0.1)
q3.EB <- lapply(0:N.new, function(X) binom.PP.EB(x.hist, n.hist, X, N.new, mix = TRUE))
q3.FB <- q2.FB
q3.fix <- q2.fix

q3.rob.05.post <- lapply(0:N.new, posterior.mixture.prior, ns=N.new, mixture.prior=q3.rob.05)
q3.rob.01.post <- lapply(0:N.new, posterior.mixture.prior, ns=N.new, mixture.prior=q3.rob.01)
q3.fix.post <- lapply(0:N.new, posterior.mixture.prior, ns=N.new, mixture.prior=q3.fix)
q3.FB.post <- lapply(0:N.new, posterior.mixture.prior, ns=N.new, mixture.prior=q3.FB)
q3.EB.post <- lapply(0:N.new, function(i) {posterior.mixture.prior( xs=i, ns=N.new, mixture.prior=q3.EB[[i+1]])})

SIG.05 <- !sig.matrix(n.control=N.new, n.treatment=N.new, level=0.025, posterior = q3.rob.05.post) 

q3c.power.05 <- calc.power(sig.mat = SIG.05, treatment.difference = 0.15-0.22, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)
q3c.t1e.05 <- calc.power(sig.mat = SIG.05, treatment.difference = 0, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)


SIG.01 <- !sig.matrix(n.control=N.new, n.treatment=N.new, level=0.025, posterior = q3.rob.01.post) 

q3c.power.01 <- calc.power(sig.mat = SIG.01, treatment.difference = 0.15-0.22, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)
q3c.t1e.01 <- calc.power(sig.mat = SIG.01, treatment.difference = 0, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)

SIG.fix <- !sig.matrix(n.control=N.new, n.treatment=N.new, level=0.025, posterior = q3.fix.post) 

q3c.power.fix <- calc.power(sig.mat = SIG.fix, treatment.difference = 0.15-0.22, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)
q3c.t1e.fix <- calc.power(sig.mat = SIG.fix, treatment.difference = 0, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)

SIG.FB <- !sig.matrix(n.control=N.new, n.treatment=N.new, level=0.025, posterior = q3.FB.post) 

q3c.power.FB <- calc.power(sig.mat = SIG.FB, treatment.difference = 0.15-0.22, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)
q3c.t1e.FB <- calc.power(sig.mat = SIG.FB, treatment.difference = 0, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)


SIG.EB <- !sig.matrix(n.control=N.new, n.treatment=N.new, level=0.025, posterior = q3.EB.post) 

q3c.power.EB <- calc.power(sig.mat = SIG.EB, treatment.difference = 0.15-0.22, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)
q3c.t1e.EB <- calc.power(sig.mat = SIG.EB, treatment.difference = 0, prob.range = c(0.1,0.5), length=41, n.binom.control = N.new)


# Plot power as a function of assumed control failure rate for new trial, and overlay power curves for each prior
plot(y=q3c.power, x=seq(0.1,0.5,len=41), xlab="Control Failure Rate (new trial)", ylab="Power", ty='l', ylim=c(0,1))
lines(y=q3c.power.05, x=seq(0.1,0.5,len=41), col=2)
lines(y=q3c.power.01, x=seq(0.1,0.5,len=41), col=3)
lines(y=q3c.power.fix, x=seq(0.1,0.5,len=41), col=4)
lines(y=q3c.power.FB, x=seq(0.1,0.5,len=41), col=5)
lines(y=q3c.power.EB, x=seq(0.1,0.5,len=41), col=6)
lines(y=q3c.power.nohist, x=seq(0.1,0.5,len=41), col=7)
abline(v=0.22, lty=2)  # mark value of control rate assumed for new trial 
legend("bottomleft", legend=c("MAP", "Robust MAP 0.5", "Robust MAP 0.1", "Fixed PP", "FB PP", "EB PP", "No historical"), lty=rep(1,7), col=1:7, cex=0.6)

# Plot type 1 error as a function of assumed control failure rate for new trial, and overlay type 1 error curves for each prior
plot(y=q3c.t1e, x=seq(0.1,0.5,len=41), xlab="Control Failure Rate (new trial)", ylab="Type 1 Error", ty='l', ylim=c(0,0.1))
lines(y=q3c.t1e.05, x=seq(0.1,0.5,len=41), col=2)
lines(y=q3c.t1e.01, x=seq(0.1,0.5,len=41), col=3)
lines(y=q3c.t1e.fix, x=seq(0.1,0.5,len=41), col=4)
lines(y=q3c.t1e.FB, x=seq(0.1,0.5,len=41), col=5)
lines(y=q3c.t1e.EB, x=seq(0.1,0.5,len=41), col=6)
lines(y=q3c.t1e.nohist, x=seq(0.1,0.5,len=41), col=7)
abline(v=0.22, lty=2)  # mark value of control rate assumed for new trial 
legend("topright", legend=c("MAP", "Robust MAP 0.5", "Robust MAP 0.1", "Fixed PP", "FB PP", "EB PP", "No historical"), lty=rep(1,7), col=1:7, cex=0.6)