<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/one-arm-bayesian-adaptive-trial-simulation-code/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://biostatmatt.com/archives/2979"> R – BioStatMatt</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>I had an opportunity recently to design a Bayesian adaptive trial with several interim analyses that allow for early stopping due to efficacy or futility. The code below implements the one-arm trial described in the great introductory article by <a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4247348/" target="_blank">Ben Saville et al.</a> Three of the coauthors are current or former colleagues at Vandy. The figures below the code show the sample size distribution under the null hypothesis and an alternative hypothesis, respectively.</p>
<pre><code>
## Simulate Bayesian single-arm adaptive trial
## Allow early termination due to futility or efficacy
## Binary outcome
## Beta-binomial:
##   p   ~ beta(a, b)
##   x_i ~ binomial(p)  i = 1..n
##   p|x ~ beta(a + sum(x), b + n - sum(x))
## Efficacy at interim t if Pr(p > p_0 | x_{(t)}) > \gamma_e
## Futility at interim t if Pr(p > p_0 | x_{(t_max)}) &lt; \gamma_f
## https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4247348/

library('rmutil') ## for betabinom
## Simulate entire trial
## ptru - true probability of outcome (p)
## pref - reference probability of outcome (p_0)
## nint - sample sizes at which to conduct interim analyses
## efft - efficacy threshold
## futt - futility threshold
## apri - prior beta parameter \alpha
## bpri - prior beta parameter \beta
simtrial &lt;- function(
  ptru = 0.15,
  pref = 0.15,
  nint = c(10, 13, 16, 19),
  efft = 0.95,
  futt = 0.05,
  apri = 1,
  bpri = 1) {
  
  ## determine minimum number of 'successes' necessary to 
  ## conclude efficacy if study continues to maximum 
  ## sample size
  nmax &lt;- max(nint)
  post &lt;- sapply(0:nmax, function(nevt)
    1-pbeta(pref, apri + nevt, bpri + nmax - nevt))
  nsuc &lt;- min(which(post > efft)-1)
  
  
  ## simulate samples
  samp &lt;- rbinom(n = nmax, size = 1, prob = ptru)
  
  ## simulate interim analyses
  intr &lt;- lapply(nint, function(ncur) {
    
    ## compute number of current events
    ecur &lt;- sum(samp[1:ncur])
    
    ## compute posterior beta parameters
    abb &lt;- apri + ecur
    bbb &lt;- bpri + ncur - ecur
    sbb &lt;- abb + bbb
    mbb &lt;- abb/(abb+bbb)
    
    ## compute efficacy Pr(p > p_0 | x_{(t)})
    effp &lt;- 1-pbeta(pref, abb, bbb)
    
    ## return for efficacy
    if(effp > efft)
      return(list(action='stop',
                  reason='efficacy',
                  n = ncur))
    
    ## number of events necessary in remainder of 
    ## study to conclude efficacy
    erem &lt;- nsuc-ecur

    ## compute success probability Pr(p > p_0 | x_{(t_max)})
    if(erem > nmax-ncur) { ## not enough possible events
      sucp &lt;- 0
    } else { ## not yet met efficacy threshold
      sucp &lt;- 1-pbetabinom(q = erem-1, 
        size = nmax-ncur, m = mbb, s = sbb)
    }
    
    if(sucp &lt; futt)
      return(list(action='stop',
                  reason='futility',
                  n = ncur))
    
    return(list(action='continue',
                reason='',
                n = ncur))
  })
  
  stpi &lt;- match('stop', sapply(intr, `[[`, 'action'))
  return(intr[[stpi]])
}


## Simulate study with max sample size of 200 where true 
## probability is identical to reference (i.e., the null
## hypothesis is true). This type of simulation helps us 
## determine the overall type-I error rate.
nint &lt;- c(40,80,120,160,200)
nmax &lt;- max(nint)
res &lt;- do.call(rbind, lapply(1:10000, 
  function(t) as.data.frame(simtrial(ptru = 0.72, 
                                     pref = 0.72,
                                     nint = nint,
                                     efft = 0.975,
                                     futt = 0.20))))

## Prob. early termination (PET) due to Futility
mean(res$reason == 'futility' &amp; res$n &lt; nmax)
## PET Efficacy
mean(res$reason == 'efficacy' &amp; res$n &lt; nmax)
## Pr(conclude efficacy) 'type-I error rate'
mean(res$reason == 'efficacy')
## average and sd sample size
mean(res$n); sd(res$n)
barplot(prop.table(table(res$n)),
        xlab='Study Size (N)',
        main=&quot;No Difference&quot;)

## Simulate study where true probability is greater than
## reference (i.e., an alternative hypothesis). This type
## of simulation helps us determine the study power.
res &lt;- do.call(rbind, lapply(1:10000, 
  function(t) as.data.frame(simtrial(ptru = 0.82, 
                                     pref = 0.72,
                                     nint = nint,
                                     efft = 0.975,
                                     futt = 0.20))))

## Prob. early termination (PET) due to Futility
mean(res$reason == 'futility' &amp; res$n &lt; nmax)
## PET Efficacy
mean(res$reason == 'efficacy' &amp; res$n &lt; nmax)
## Pr(conclude efficacy) 'power'
mean(res$reason == 'efficacy')
## average and sd sample size
mean(res$n); sd(res$n)
barplot(prop.table(table(res$n)),
        xlab='Study Size (N)',
        main=&quot;35% Reduction&quot;)

</code></pre>
<p><a rel="nofollow" href="https://i0.wp.com/biostatmatt.com/uploads/null.png" target="_blank"><img alt="" src="https://i0.wp.com/biostatmatt.com/uploads/null.png?w=450" class="wp-image-2983" data-recalc-dims="1"></a></p>
<p> <a rel="nofollow" href="https://i2.wp.com/biostatmatt.com/uploads/alt.png" target="_blank"><img alt="" src="https://i2.wp.com/biostatmatt.com/uploads/alt.png?w=450" class="wp-image-2982" data-recalc-dims="1"></a></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://biostatmatt.com/archives/2979"> R – BioStatMatt</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/HFpAfLaPCPY" height="1"></a></div>