Assessment of Bayesian Expected Power via Bayesian Bootstrap

TL;DR

This paper introduces a model-free Bayesian bootstrap method to estimate Bayesian expected power for future trials, offering a computationally efficient and robust alternative to traditional model-based approaches, with demonstrated effectiveness through simulations and case studies.

Contribution

The paper proposes a Bayesian bootstrap approach for BEP estimation that avoids model assumptions and simplifies computation, also introducing a frequentist double bootstrap method.

Findings

BBS provides accurate BEP estimates comparable to model-based methods.

BBS is computationally simpler and more robust to model misspecification.

Double bootstrap (BS2) offers a frequentist alternative with similar properties.

Abstract

The Bayesian expected power (BEP) has become increasingly popular in sample size determination and assessment of the probability of success (POS) for a future trial. The BEP takes into consideration the uncertainty around the parameters assumed by a power analysis and is thus more robust compared to the traditional power that assumes a single set of parameters. Current methods for assessing BEP are often based in a parametric framework by imposing a model on the pilot data to derive and sample from the posterior distributions of the parameters. Implementation of the model-based approaches can be analytically challenging and computationally costly especially for multivariate data sets; it also runs the risk of generating misleading BEP if the model is mis-specified. We propose an approach based on the Bayesian bootstrap technique (BBS) to simulate future trials in the presence of individual-level pilot data, based on which the empirical BEP can be calculated. The BBS approach is model-free with no assumptions about the distribution of the prior data and circumvents the analytical and computational complexity associated with obtaining the posterior distribution of the parameters. Information from multiple pilot studies is also straightforward to combine. We also propose the double bootstrap (BS2), a frequentist counterpart to the BBS, that shares similar properties and achieves the same goal as the BBS for BEP assessment. Simulation studies and case studies are presented to demonstrate the implementation of the BBS and BS2 techniques and to compare the BEP results with model-based approaches.

Tables

Table 1. Table 1: BEP assessment via the BBS

DO $j=1,\dots ,m$

1) draw ${𝝅}^{\left(j\right)}$ from the Dirichlet distribution in Eqn (3);

2) draw ${\widetilde{𝐲}}^{\left(j\right)}$ from $\text{M}(\tilde{n},{𝝅}^{\left(j\right)})$;

3) test $H_0$ in simulated trial ${\widetilde{𝐲}}^{\left(j\right)}$ and record $I\left(R\left({\widetilde{𝐲}}^{\left(j\right)}\right)\right)$ (1 if $H_0$ is rejected; 0 otherwise);

END DO

OUTPUT: calculate BEP $=m^{-1}{\sum }_{j=1}^{m}I\left(R\left({\widetilde{𝐲}}^{\left(j\right)}\right)\right)$.

Table 2. Table 2: Generation of the posterior distribution of power via the BBS

DO $j=1,\dots ,m$

1) draw ${𝝅}^{\left(j\right)}$ from the Dirichlet distribution in Eqn (3);

2) given the drawn ${𝝅}^{\left(j\right)}$

DO $l=1,\dots ,t$

2.1) draw ${\widetilde{𝐲}}^{(l,j)}$ from $\text{M}(\tilde{n},{𝝅}^{\left(j\right)})$;

2.2) test $H_0$ in simulated trial ${\widetilde{𝐲}}^{(l,j)}$ and record $I\left(R\left({\widetilde{𝐲}}^{(l,j)}\right)\right)$ (1 if $H_0$ is rejected; 0 o.w.);

END DO

3) calculate power $\beta \left({𝝅}^{\left(j\right)}\right)=t^{-1}{\sum }_{l=1}^{t}I\left(R\left({\widetilde{𝐲}}^{(l,j)}\right)\right)$.

END DO

OUTPUT: $m$ samples from the posterior distribution of power: $(\beta \left({𝝅}^{\left(1\right)}\right),\dots ,\beta \left({𝝅}^{\left(m\right)}\right))$

and BEP $=m^{-1}{\sum }_{j=1}^m(\beta \left({𝝅}^{\left(j\right)}\right)$.

Table 3. Table 3: Generation of the conditional distribution of power given pilot data via the BS2

DO $j=1,\dots ,m$

1) bootstrap a sample $S^{\left(j\right)}$ of size $n$ from $𝐲$ with replacement;

DO $l=1,\dots ,t$

2.1) bootstrap ${\widetilde{𝐲}}^{(l,j)}$ of size $\tilde{n}$ from $S^{\left(j\right)}$ with replacement;

2.2) test $H_0$ in simulated trial ${\widetilde{𝐲}}^{(l,j)}$ and record $I\left(R\left({\widetilde{𝐲}}^{(l,j)}\right)\right)$ (1 if $H_0$ is rejected; 0 o.w.);

END DO

3) calculate power $\beta \left({𝝅}^{\left(j\right)}\right)=t^{-1}{\sum }_{l=1}^{t}I\left(R\left({\widetilde{𝐲}}^{(l,j)}\right)\right)$.

END DO

OUTPUT: $m$ samples from the conditional distribution of power given $𝐲$:

$(\beta \left({𝝅}^{\left(1\right)}\right),\dots ,\beta \left({𝝅}^{\left(m\right)}\right))$, and EP $=m^{-1}{\sum }_{j=1}^m(\beta \left({𝝅}^{\left(j\right)}\right)$.

Table 4. Table 4: Monte Carlo power via FTS given pilot data

DO $j=1,\dots ,m$

1) bootstrap ${\widetilde{𝐲}}^{\left(j\right)}$ of size $\tilde{n}$ from $𝐲$ with replacement;

2) test $H_0$ in ${\widetilde{𝐲}}^{\left(j\right)}$ and record $I\left(R\left({\widetilde{𝐲}}^{\left(j\right)}\right)\right)$ (1 if $H_0$ is rejected; 0 o.w.);

END DO

OUTPUT: power $=m^{-1}{\sum }_{j=1}^{m}I\left(R\left({\widetilde{𝐲}}^{\left(j\right)}\right)\right)$.

Table 5. Table 5: POS assessment of each test formulation in the future BE study

	BEP (%)					Power (%)
treatment	model-	model-	BBS	BS2		model-	model-	BS
	based${}_{}{}^{†}{}_1{}^{}$	based${}_{}{}^{†}{}_2{}^{}$				based${}_{}{}^{‡}{}_1{}^{}$	based${}_{}{}^{‡}{}_2{}^{}$
T1	51.0	74.3	86.3	87.6		99.98	65.8	100.0
T2	42.6	68.0	79.4	81.6		99.39	78.6	99.3

Table 6. Table 6: Empirical POS assessment for the future trial to compare two antiretroviral drugs in hazard of death in HIV patients

$\tilde{n}$	$T$ was the same as pilot study							$T=24$ months
allocation	BEP(%)				Power (%)			BEP(%)	Power†(%)
ratio 1:1	model-based	BBS	BS2		model-based†	BS		model-based
200	73.4	75.0	73.6		99.6	90.4		83.0	$>$99.9
400	81.8	83.0	84.5		$>$99.9	99.2		88.3	100.0
600	85.5	87.0	86.5		100.0	99.8		90.5	100.0

Equations

$$\text{mbox}BEP=\text{mbox}E[\beta (\mathbf{\theta })|y]=\int \beta (\mathbf{\theta })p(\mathbf{\theta }|\mathbf{y})d\mathbf{\theta }.$$

$$=$$

$$p({\pi }_1,\dots ,{\pi }_n|\mathbf{y})=\frac{{\prod }_{i=1}^n{\pi }_k^{({\alpha }_k+1)-1}}{B({\alpha }_1+1,\dots ,{\alpha }_n+1)}=\text{mbox}D({\alpha }_1+1,\dots ,{\alpha }_n+1).$$

$$p({\pi }_1,\dots ,{\pi }_n|\mathbf{y})={\left(\Gamma (n)\right)}^{-1}=\text{mbox}D(1,\dots ,1).$$

$$p(\tilde{\mathbf{y}}|\mathbf{y})=\int p(\tilde{\mathbf{y}}|\mathbf{y},\mathbf{\pi })p(\mathbf{\pi }|\mathbf{y})d\mathbf{\pi }=\int p(\tilde{\mathbf{y}}|\mathbf{\pi })p(\mathbf{\pi }|\mathbf{y})d\mathbf{\pi }=\int M(\tilde{n},\mathbf{\pi })p(\mathbf{\pi }|\mathbf{y})d\mathbf{\pi },$$

$$\text{mbox}BEP=\int \int I(R(\tilde{\mathbf{y}}))M(\tilde{n},\mathbf{\pi })p(\mathbf{\pi }|\mathbf{y})d\mathbf{\pi }d\tilde{\mathbf{y}}$$

$$L(\mathbf{\beta },\Sigma ;\mathbf{y},\mathbf{x})={\prod }_{i=1}^n{\left(2\pi |\Sigma |\right)}^{-1/2}\exp \{({\mathbf{y}}_i-{\mathbf{x}}_i^{\prime }\mathbf{\beta }{)}^{\prime }{\Sigma }^{-1}({\mathbf{y}}_i-{\mathbf{x}}_i^{\prime }\mathbf{\beta })\}.$$

$$p(\mathbf{\beta },\Sigma |\mathbf{y},\mathbf{x})\propto |\Sigma {|}^{-({\nu }_0+p+1+n)/2}\exp \{({\mathbf{y}}_i-{\mathbf{x}}_i^{\prime }\mathbf{\beta }{)}^{\prime }{\Sigma }^{-1}({\mathbf{y}}_i-{\mathbf{x}}_i^{\prime }\mathbf{\beta })\}.$$

$$p({\mathbf{y}}_{\text{mbox}mis}|\mathbf{\beta },\Omega ,{\mathbf{y}}_{\text{mbox}obs})=N\left({\mathbf{x}}_{\text{mbox}mis}\mathbf{\beta },M\Omega M^{\prime }\right),$$

$$p(\Sigma |\mathbf{\beta },{\mathbf{y}}_{\text{mbox}mis},{\mathbf{y}}_{\text{mbox}obs})=\text{mbox}Inv-Wishart({\sum }_{i=1}^n{\mathbf{e}}_i{\mathbf{e}}_i^{\prime },{\nu }_0+n)$$

$$=$$

$$p(\mathbf{\beta }|\Omega ,{\mathbf{y}}_{\text{mbox}mis},{\mathbf{y}}_{\text{mbox}obs})=N\left(({\mathbf{x}}^{\prime }{\Omega }^{-1}\mathbf{x}{)}^{-1}{\mathbf{x}}^{\prime }{\Omega }^{-1}\mathbf{y},({\mathbf{x}}^{\prime }{\Omega }^{-1}\mathbf{x}{)}^{-1}\right).$$

$$y_{ij}(t)$$

$$h_i(t)$$

$$\log (h_0(t))$$

$$p(\mathbf{\Sigma }|\mathbf{\beta },\mathbf{y})$$

$$\propto |\Sigma {|}^{-({\nu }_0+p+n+1)/2}\text{mbox}exp\left(-\frac{1}{2}{\sum }_{i=1}^n\text{mbox}tr({\mathbf{e}}_i{\mathbf{e}}_i^T{\Sigma }^{-1})\right)$$

$$\propto |\Sigma {|}^{-({\nu }_0+p+n+1)/2}\text{mbox}exp\left(-\frac{1}{2}\text{mbox}tr\left({\sum }_{i=1}^n{\mathbf{e}}_i{\mathbf{e}}_i^T{\Sigma }^{-1}\right)\right),\text{mbox}where{\mathbf{e}}_i={\mathbf{y}}_i-{\mathbf{x}}_i\mathbf{\beta }.$$

$$\text{mbox}Thusp(\mathbf{\Sigma }|\mathbf{\beta },\mathbf{y})$$

$$p(\mathbf{\beta }|\mathbf{\Sigma },\mathbf{y})$$

$$\propto |\Sigma {|}^{-({\nu }_0+p+n+1)/2}\text{mbox}exp\left(-\frac{1}{2}{\sum }_{i=1}^n{\left({\Sigma }^{-1/2}({\mathbf{y}}_i-{\mathbf{x}}_i\mathbf{\beta })\right)}^T\left({\Sigma }^{-1/2}({\mathbf{y}}_i-{\mathbf{x}}_i\mathbf{\beta })\right)\right)$$

$$\propto |\Sigma {|}^{-({\nu }_0+p+n+1)/2}\text{mbox}exp(-\frac{1}{2}({\mathbf{\beta }}^T({\sum }_{i=1}^n{\mathbf{x}}_i^T{\Sigma }^{-1}{\mathbf{x}}_i)\mathbf{\beta }$$

$$-{\sum }_{i=1}^n{\mathbf{y}}_i^T{\Sigma }^{-1}{\mathbf{x}}_i\mathbf{\beta }-{\mathbf{\beta }}^T{\sum }_{i=1}^n{\mathbf{x}}_i^T{\Sigma }^{-1}{\mathbf{y}}_i+{\sum }_{i=1}^n{\mathbf{y}}_i^T{\Sigma }^{-1}{\mathbf{y}}_i))$$

$$\propto |\Sigma {|}^{-({\nu }_0+p+n+1)/2}$$

$$\text{mbox}exp\left(-\frac{1}{2}{\left(\mathbf{\beta }-({\mathbf{x}}^T{\Omega }^{-1}\mathbf{x}{)}^{-1}({\mathbf{x}}^T{\Omega }^{-1}\mathbf{y})\right)}^T{\left(({\mathbf{x}}^T{\Omega }^{-1}\mathbf{x}{)}^{-1}\right)}^{-1}\left(\mathbf{\beta }-({\mathbf{x}}^T{\Omega }^{-1}\mathbf{x}{)}^{-1}({\mathbf{x}}^T{\Omega }^{-1}\mathbf{y})\right)\right).$$

Full text

Assessment of Bayesian Expected Power via Bayesian Bootstrap

Fang Liu*∗*

Applied and Computational Mathematics and Statistics

University of Notre Dame, Notre Dame, IN 46556

∗email: fang.liu.131@nd.edu

The Bayesian expected power (BEP) has become increasingly popular in sample size determination and assessment of the probability of success (POS) for a future trial. The BEP takes into consideration the uncertainty around the parameters assumed by a power analysis and is thus more robust compared to the traditional power that assumes a single set of parameters. Current methods for assessing BEP are often based in a parametric framework by imposing a model on the pilot data to derive and sample from the posterior distributions of the parameters. Implementation of the model-based approaches can be analytically challenging and computationally costly especially for multivariate data sets; it also runs the risk of generating misleading BEP if the model is mis-specified. We propose an approach based on the Bayesian bootstrap technique (BBS) to simulate future trials in the presence of individual-level pilot data, based on which the empirical BEP can be calculated. The BBS approach is model-free with no assumptions about the distribution of the prior data and circumvents the analytical and computational complexity associated with obtaining the posterior distribution of the parameters. Information from multiple pilot studies is also straightforward to combine. We also propose the double bootstrap (BS2), a frequentist counterpart to the BBS, that shares similar properties and achieves the same goal as the BBS for BEP assessment. Simulation studies and case studies are presented to demonstrate the implementation of the BBS and BS2 techniques and to compare the BEP results with model-based approaches.

KEY WORDS: double bootstrap; trial simulation; probability of success; robust decision making; weighted average power; assurance

1 Introduction

The probability of success (POS) of a clinical trial is affected by many factors such as patient recruitment, ethical considerations, local regulations, resources, study designs and execution, among others. The part where statisticians get involved the most is the study design during the trial planning stage, including the sample size determination and power calculation. Power is defined as the probability of rejecting the null hypothesis $H_{0}$ if the alternative hypothesis $H_{1}$ is true in a future trial. The probabilistic nature of power makes a natural choice as a metric for measuring how the statistical design aspect of the study affects the POS of the trial. On the other hand, the classical power is a conditional probability given a specific effect size $\Delta=\mu/\sigma$, the value of which is an unknown quantity and is often the parameter of primary interest to estimate in the planned trial. As such, the power is sensitive to the assumed effect size and not a robust measurement of POS.

A more robust measurement for POS compared to the classical power is the Bayesian expected power (BEP) or the weighted average power (WAP), which is the expected power over the distribution of the effect size $\Delta$ given existing data. In other words, the BEP is a “marginal” measure of the POS given what’s known, integrating out the unknown $\Delta$. The BEP has been “re-invented” several times under different names in various contexts. Brown et al., (1987) suggested using Bayesian methods to obtain a posterior distribution representing the state of knowledge of the parameters of interest, and to predict the outcome of a specified comparative trial. The approach was implemented and demonstrated in several examples in Spiegelhalter et al., (2004). O’Hagan and Stevens, (2001) and O’Hagan et al., (2005) used the term “assurance” and Chuang-Stein, (2006) used “average success probability” to describe the BEP. In all thedr cases, only the uncertainty around $\mu$ concerned by the hypothesis in the future trial is considered and the uncertainty around the variance $\sigma$, which is an equally important parameter in power calculation, is not taken into consideration. Liu, (2010) extended the BEP by accounting for the uncertainty around both $\mu$ and $\sigma$ (or the effect size $\Delta$), and also added two more versions of the BEP to by removing the “type I error” component from the regular BEP metric. Some recent reviews, discussion, and applications of the BEP in early and late phased clinical trials, and in meta-analysis are given in Kirby et al., (2012); Carroll, (2013); Ibrahim et al., (2015); Du and Wang, (2016); Zierhut et al., (2016). The BEP has been routinely calculated alongside the traditional power in some pharmaceutical companies; and it is also implemented in several sample size and power calculation software (Labes et al.,, 2016; EAST-CYTEL,, 2016).

The current approaches to assessing the BEP often start with constructing the joint posterior distributions of $\mu$ and $\sigma$ given a prior and the likelihood of the parameters given the pilot/exisiting data. In many cases, the posterior distributions of $\mu$ and $\sigma$ may not have closed-form expressions, or even they do, iterative approaches, such as the MCMC algorithms and other sampling techniques, might still be necessary to draw the parameters from the posterior distributions. In the case when there are co-primary hypotheses, $\boldsymbol{\mu}$ is multi-dimensional, analytically and computational it becomes even harder. In summary, the model-based approaches can be analytically challenging and computationally costly especially for multivariate data. On top of all these, there is always the risk of misspecifying the likelihood of the parameters with the pilot data, leading to misleading BEP values subsequently.

We propose an approach based on the Bayesian bootstrap technique (Rubin,, 1981) , referred to as the BBS approach, to assess the BEP when individual-level pilot data $\mathbf{y}$ are available. The Bayesian bootstrap is a Bayesian version of the bootstrap technique (Efron,, 1979). Rubin, (1981) proved that the Bayesian bootstrap is operationally and inferentially similar to the frequentist bootstrap. In the BBS approach, we repeatedly simulate the future trial data given $\mathbf{y}$ via the Bayesian bootstrap technique without imposing any distributional assumptions on $\mathbf{y}$, and test $H_{0}$ according to the planned analysis in each set of the simulated trials; and the overall rejection rate of $H_{0}$ over the repetitions leads to a Monte Carlo (MC) estimate of the BEP. The uncertainty around the underlying true distribution of the pilot data is accounted for by placing a Dirichlet prior on the probability of each individual. The BBS approach is straightforward to implement with minimal analytical work except for the planned analysis on the simulated future data. Computationally, only a few lines of codes are needed to implement the sampling and trial simulation steps, and the whole BBS procedure can be easily parallelized for fast computation. When there are multiple relevant pilot data sets, it is straightforward to combine the pilot data together via the BBS approach.

Directly applying the regular bootstrap to the pilot data to simulate future data is inappropriate for the purpose of BEP assessment since it will not propagate the uncertainty around the distribution of the pilot data. As a matter of fact, the overall rejection rate based on the future data simulated this way is a MC estimate of the classical power assuming the estimated effect size from the pilot data is the true effect size. This motivates us to come up with the “double bootstrap” technique (BS2), a frequentist counterpart to the BBS, that will propagate the uncertainty around the distribution of the pilot data and achieves the same goal as the BBS for BEP assessment. Procedurally, the BS2 first bootstraps the pilot data, and then samples the bootstrapped pilot data to generate repetitions for the future trials based on which a MC estimate of the BEP is obtained. Similar to the BBS, only a few lines of codes are needed to implement the BS2, and it is easy to parallelize computationally.

The rest of the paper is organized as follows. Section 2 introduces the BBS and BS2 approaches for assessing the BEP. Section 3 compares the BBS and BS2 approaches with the parametric approaches in the BEP assessment in two simulation studies. Section 4 implements the BBS and BS2 approaches to assess BEP for an equivalence study and for a HIV survival study. Section 5 concludes the discussion with some final remarks. The R codes for the numerical examples in Section 4 are provided in the online supplementary materials to this paper and are also available for download at TBD.

2 Methods

2.1 assessment of BEP

By definition, the BEP is the probability of rejecting $H_{0}$ in a future trial given existing data $\mathbf{y}$. Denote the classical power by $\beta(\boldsymbol{\theta})=\Pr$(rejecting $H_{0}|\boldsymbol{\theta},\tilde{n}$), where $\boldsymbol{\theta}$ refers to the parameters involved in the power calculation and $\tilde{n}$ is the given sample size of the future trial, then

[TABLE]

To calculate the BEP, we may first draw $\boldsymbol{\theta}$ from its posterior distribution $p(\boldsymbol{\theta}|\mathbf{y})$, and then plug in the drawn $\boldsymbol{\theta}$ in the power function $\beta(\boldsymbol{\theta})$ to obtain a classical power value given the drawn $\boldsymbol{\theta}$. Repeating the two steps many times, say $m$, will lead to the posterior distribution of the power $\beta$ given $\mathbf{y}$, the average $\sum_{j=1}^{m}\beta(\boldsymbol{\theta^{(j)}})$ of the posterior samples is an MC estimate for BEP for the future trial for a given $\tilde{n}$. Besides the BEP, other statistics from the posterior distribution of the power can also be reported, such as the mode and the percentiles.

If a closed-formed function $\beta(\boldsymbol{\theta})$ is not available, the MC method can be applied to numerically approximate the power by simulating the future trial data $\tilde{\mathbf{y}}$ on which $H_{0}$ will be tested. Denote the event of rejecting $H_{0}$ in the future trial by $I(R(\tilde{\mathbf{y}}))$, where $R$ is the rejection rule, a function of $\tilde{\mathbf{y}}$, and $I(R(\tilde{\mathbf{y}}))=1$ if $H_{0}$ is rejected in $\tilde{\mathbf{y}}$ and 0 otherwise, then

[TABLE]

Eqn (2) suggests that we can first draw $\boldsymbol{\theta}$ from the posterior distribution $p(\boldsymbol{\theta}|\mathbf{y})$ given pilot data $\mathbf{y}$, and simulate the future data $\tilde{\mathbf{y}}$ of size $\tilde{n}$ given the drawn $\boldsymbol{\theta}$, then perform the planned statistical analysis and hypothesis testing on the simulated data $\tilde{\mathbf{y}}$, and record 1 if rejected and 0 if not. Repeating the process many times, say $m$, and the rejection rate over the $m$ repetations is an MC estimate of the BEP. We refer to this MC approach for assessing the BEP as the future trial simulation (FTS) approach.

Regardless of whether the model-based approach is employed when closed-form $\beta(\boldsymbol{\theta})$ is available or the FTS approach is applied when it is not, the current approaches for the BEP assessment are mainly model-based by imposing distributional or model assumptions on the pilot data $\mathbf{y}$ to obtain the posterior distribution $p(\boldsymbol{\theta}|\mathbf{y})$, followed by the posterior sampling step. If $\tilde{\mathbf{y}}$ and $\boldsymbol{\theta}$ are multidimensional, the sampling of $\boldsymbol{\theta}$ and simulation of $\tilde{\mathbf{y}}$ can be computationally costly, not to be mention the risk of mis-specification of the parametric model, leading to misleading BEP subsequently. Alternative techniques that avoid making strong parametric assumptions about data $\mathbf{y}$ and at the same time are computationally less complicated are desired.

We propose a FTS technique, referred to as the BBS approach, motivated by the Bayesian bootstrap to evaluate the BEP numerically. The BBS assesses the BEP without imposing a model on the pilot data or sampling from complicated posterior distributions via MCMC approaches. Instead, the BBS approach obtains the posterior distributions of the population distribution underlying $\mathbf{y}$ given $\mathbf{y}$, from which the future data will be simulated. We will also present the frequentist version of the BBS approach – the “double bootstrap” technique (BS2), and discuss why the regular bootstrap will only lead to a MC estimation of the classical power instead of yielding the BEP.

2.2 construction of a posterior distribution of the population distribution via Bayesian bootstrap

Define a finite population of size $N$ with $K$ ($K<\infty$) distinct values over $p$ attributes. Let $\pi_{k}$ denote the probability that value $\mathbf{d}_{k}$ occurs for $k=1,\ldots,K$. The Bayesian bootstrap can be used to obtain a posterior distribution of the population distribution $f$ of the $K$ distinct values given a data set $\mathbf{y}$ of size $n$. Denote by $\boldsymbol{\pi}=(\pi_{1},\ldots,\pi_{K})$ the probabilities associated with the $K$ distinct values ($\sum_{k}\pi_{k}=1$). Since $\boldsymbol{\pi}$ fully characterizes $f$, obtaining the posterior distribution on $f$ is equivalent to obtaining the posterior distribution of $\boldsymbol{\pi}$. A convenient choice on the prior of $\boldsymbol{\pi}$ is the Dirichlet distribution $\mbox{D}(\boldsymbol{\pi}|\alpha_{1},\ldots,\alpha_{K})\propto\prod_{k=1}^{K}\pi_{k}^{\alpha_{k}-1}$, with hyper-parameters $\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{K})>0$. Since Dirichlet priors are conjugate priors for the multinomial likelihood, the posterior distribution of $\boldsymbol{\pi}$ given data $\mathbf{y}$ (of size $n$) also follows a Dirichlet distribution $\mbox{D}(\pi_{1},\ldots,\pi_{K}|\mathbf{n})\propto\prod_{k=1}^{K}\pi_{k}^{n_{k}+\alpha_{k}-1}$, where $\mathbf{n}=(n_{1},\ldots,n_{K})$, $n_{k}$ is the number of observations in the sample that take the value $\mathbf{d}_{k}$, and $\sum_{k}n_{k}=n$.

It is highly likely that not every distinct $\mathbf{d}_{k}$ from the population will occur in the sample; and $n_{k}$’s associated with these “non-appearing” cases will be [math]. Denote by $L$ the number of distinct cases in the sample ($L\leq K$), the posterior distribution of the subset of $\boldsymbol{\pi}$ that is associated with the sample cases is $D(\pi_{1},\ldots,\pi_{L}|\mathbf{n})\propto\prod_{k=1}^{L}\pi_{k}^{n_{k}+\alpha_{k}-1}$. In practical application, we can always set $L=n$ due to two reasons: it does not affect the posterior distribution of $\boldsymbol{\pi}$ associated with the distinct values given the aggregation property of the Dirichlet distribution; and every individual in the sample can be distinct from each other when there are continuous variables among the $p$ attributes or $p$ is large. With $L=n$, $n_{k}=1$ for $k=1,\ldots,n$, and The posterior distribution of $\boldsymbol{\pi}$ can be simplified to

[TABLE]

In terms of the choice for the hyper-parameters $\alpha_{k}$, proper priors should have $\alpha_{k}>0\;\forall\;k=1,\ldots,n$. If there are minimal information about $\boldsymbol{\pi}$ prior to the pilot data $\mathbf{y}$, setting $\alpha_{k}$ at small positive numbers (eg, 0.1, 0.5, 1, etc) leads to weakly informative priors for $\boldsymbol{\pi}$. An improper non-informative but convenient choice is $\alpha_{k}=0\;\forall\;k=1,\ldots,n$, which still generates a proper posterior distribution for $\boldsymbol{\pi}$,

[TABLE]

2.3 BEP assessment with Bayesian bootstrap

The Bayesian bootstrap is employed in the BEP assessment to simulate data $\tilde{\mathbf{y}}$ for the future trial, the posterior predictive distribution of which, given pilot data $\mathbf{y}$, is

[TABLE]

where $p(\boldsymbol{\pi}|\mathbf{y})$ is given in Eqn (3) (Eqn (4) if $\alpha_{k}=0$). Plugging in $p(\tilde{\mathbf{y}}|\mathbf{y})$ from Eqn (5) in Eqn (2), we have

[TABLE]

Eqn (6) suggests that BEP can be calculated via the steps given in Table 1, which produces a single metric of BEP rather than the posterior distribution of power. If it is desired to have the posterior distribution of power, an inner loop will need to be built into the BBS procedure, as given in Table 2. When the number of iterations of the inner loop $t=1$, Table 2 reduces to Table 1.

2.4 the double boostrap

The BBS procedures in Tables 1 and 2 for assessing the BEP are based in the framework of the Bayesian bootstrap. The frequentist counterpart to the BBS approach is the “double bootstrap” (BS2) given in Table 3. The outer-loop bootstrap captures the uncertainty around the population distribution (or the population parameters), corresponding to the sampling of $\boldsymbol{\pi}$ from its posterior distribution in the BBS procedure. The inner-loop bootstrap propagates the sampling variability and error for the future trial, serving the same purposes of drawing a sample data set given $\boldsymbol{\pi}$ in the BBS procedure. The BS2 technique outputs samples from a distribution of power given pilot data $y$, the average of which gives the WAP, taking into account the uncertainty of the sample data. When $m\rightarrow\infty$, the BS2 technique is equivalent to the BBS approach when the hyper-parameters $\alpha_{k}=0\;\forall\;k=1,\ldots,n$. Strictly speaking, the BS2 is not Bayesian conceptually, and the distribution of the power generated from the procedure is thus not a posterior distribution though it is a conditional distribution given the pilot data $\mathbf{y}$, and the mean of which is asymptotically equivalent to the BEP. Setting $t=1$ for the inner loop in Table 3 yields a single WAP estimate without a distribution of power (similar to Table 2 being reduced to Table 1 in the BBS).

If we directly simulate the future data $\tilde{\mathbf{y}}$ via the regular bootstrap, we end up having a procedure (Table 4) that leads to a MC estimate of the classical power assuming what’s observed in the pilot data is the truth.

It is obvious that this procedure does not take account the uncertainty around the unknown population distribution underlying $\mathbf{y}$, thus implicitly assumes the observed pilot data is the whole population and the estimated effect size from $\mathbf{y}$ is the true effect size of the population.

3 Simulation Studies

We run two simulation studies to implement the BBS and BS2 procedures in Tables 2 and 3 to assess the BEP, and compared the results to the model-based BEP with an assumed model on the pilot data. We also computed the classical power we compared the MC power via the bootstrap procedure in Table 4 to the model-based power assuming the observed results from the pilot study is the truth. The type I error rates in the hypothesis testing in future trials were 5% in both simulation studies.

In the first simulation study, 2500 repetitions were run. In each repetition, a pilot of study of size $n$ was simulated from N$(\mu=0.15,\sigma=1)$. We examined 3 cases of $n$: $n=$10, 30 and 100, respectively. The future trial under planning had $\tilde{n}=500$ and the hypotheses were $H_{0}\!:\!\mu=0$ vs $H_{1}\!:\!\mu\neq 0$ . The power assuming the effect size from the pilot study was true and the BEP were calculated analytically and via the bootstrap procedures. Figure 1 depicts the distributions for the differences over the 2500 repetitions between the various metrics, and the results are summarized as follows. First, the BEP calculated analytically and via the BBS and BS2 was similar; the small discrepancy between the two decreased with $n$ (the first and second rows). Second, the power calculated via the bootstrap procedure was similar to the model-based power; the small discrepancy between the two decreased with the pilot study size $n$ (the first row). Third, the BEP could be larger or smaller than power, depending on the pilot data. The mean difference between the two increased with $n$ and the dispersion of the difference decreased (the fourth row). From the histograms and the empirical CDFs (cumulative distribution functions) of the power and BEP presented in the online supplementary materials, we observe that the CDFs of the power and the BEP intersected when power and BEP around $0.5$ regardless of $n$. As $n$ increased, the distribution of the power became less ”bi-polar” while that of the BEP became less uniform, and the two distributions became more similar.

In the second simulation study, 1000 repetitions were run. In each repetition, a pilot of study of size $n$ was simulated from a bivariate lognormal distribution LN$(\boldsymbol{\mu}=(3,5)^{\prime},\Sigma=(^{1}_{0.65}\;{}^{0.65}_{1}))$. We examined 3 cases of $n$: 10, 30, and 50. The future study had $\tilde{n}=80$ and the hypotheses were $H_{0}\!:\!(\mu_{1}\!\leq\!2.7)\cup(\mu_{2}\!\leq\!4.5)$ vs. $\!\!H_{1}\!:\!(\mu_{1}\!>\!2.7)\cap(\mu_{2}\!>\!4.5)$. Power assuming the effect size from the pilot study was true was calculated both analytically and via the bootstrap procedure, and the model-based BEP with a correctly specified model on the pilot data (a lognormal distribution) and with a wrongly specified model (a normal distribution), and the BEP via the BBS and BS2 procedures were calculated. Figure 2 depicts the distributions of the differences over the 1000 repetitions between the metrics. The histograms of the metrics (2 powers, and 4 BEPs) over the 1000 simulations and their empirical CDFs are available in the online supplementary materials. The results on the comparisons between the model-based (correctly specified) vs bootstrap power, between the model-based (correctly specified), BBS, and BS2 BEP, and between the power and BEP are similar to those in Simulation 1. In terms of the model-based BEP based on the correctly specified model vs the wrongly specified model on the pilot data, Figure 2 (the first row) suggests the BEP based on the wrong model was smaller than the correct BEP by 2.4% to 5.6% on average, with larger $n$ blunting the effect of the model mis-specification on the BEP assessment.

4 Cast Studies

We implemented the BBS and the BS2 procedures in two case studies. The future trial in the first example is a 2-period crossover study with an equivalence hypothesis and the pilot study is a 3-period crossover study. In the second example, the future trial is a two-arm study comparing the survival time of AIDS patients on two different treatments via the joint modelling of the survival data and the longitudinal data on CD4 counts; and the pilot study is of the same design but is smaller in size. The type I error rates in both examples are 5%.

4.1 example 1: crossover study with equivalence hypothesis

A 2-period crossover study (sample size $\tilde{n}=200$) is under planning with the goal of developing a new formulation for a fixed-dose combination (FDC) drug for treating dislipidemia that has similar pharmacokinetic (PK) profile as a reference formulation R. The $H_{1}$ states that the geometric mean ratios (GMRs) between the FDC and R are within the interval of (0.80, 1.25) on 4 PK endpoints: 2 measures on the area under the PK curve (AUC) and 2 measures of the maximum concentration (Cmax) on two chemical entities. There exists a small pilot 3-period crossover study ($n=36$) with two candidate FDC formulations FDC1 and FDC2, and the reference formulation R. The team aims to use the test formulation with the higher POS in the planned trial based on the information collected in the pilot study.

When modelling a single normally distributed endpoint $\mathbf{y}_{i}=(y_{i1},\ldots,y_{ip})^{T}$ from a crossover study with $p$ repeated measures, the linear mixed-effects (lme) model $\mathbf{y}_{i}=\mathbf{x}_{i}\boldsymbol{\beta}+\mathbf{z}_{i}\boldsymbol{\gamma}_{i}+\epsilon_{i}$ is often used. The fixed-effects term $\boldsymbol{\beta}$ often includes an intercept, treatment effects, and period effects. The random effects term $\boldsymbol{\gamma}_{i}\sim\mbox{N}(0,\mathbf{G})$ and the error term $\epsilon_{i}\sim\mbox{N}(0,R)$ together define the variance/covariance V$(\mathbf{y}_{i})=\Sigma=\mathbf{z}_{i}G\mathbf{z}_{i}+R$; and $\Sigma$ can be as simple as “compound symmetry” (CS) structure (constant marginal variance on the diagonal and constant covariance on all the off-diagonals) or as complex as “unstructured” (UN) (fully parameterized with $p(p+1)/2$ parameters).

In the model-based assessment of BEP, we first obtained the posterior distribution of the treatment differences (on the log-scale) between T1 and R, and between T2 and R given the pilot data (AUC and Cmax in BE studies are often analyzed on the log scale). The hypothesis to be established in the future trial is a union of 4 BE hypotheses, one per primary endpoint. Therefore, the power is the joint probability that the 95% CIs for the 4 GMRs fell within the interval of (0.80, 1.25) simultaneously. To model the dependency structure among the 4 endpoints, we analyzed all 4 endpoints in one lme model with a fully parameterized covariance matrix $\Sigma_{12\times 12}$. The likelihood function of $(\boldsymbol{\beta},\Sigma)$ was

[TABLE]

We imposed the prior $p(\boldsymbol{\beta},\Sigma)\!\propto\!|\Sigma|^{-(\nu_{0}+p+1)/2}\!$ on $\Sigma$, which is an inverse Wishart distribution with the a priori degrees of freedom $\nu_{0}$, scale matrix $|\Sigma_{0}|\rightarrow 0$, and $p=12$. We tried two different $\nu_{0}$ at [math] and $p+1$, respectively, to examine the impact of different parametric assumptions on the model-based BEP. $\nu_{0}=0$ corresponds to the (improper) Jeffreys prior. When $\nu_{0}=p+1$, the conditional posterior mean $\mbox{E}(\Sigma|\boldsymbol{\beta},\mathbf{y})$ is the same as the MLE for $\Sigma$ if $\boldsymbol{\beta}$ is known. The joint posterior distribution of $(\boldsymbol{\beta},\Sigma)$ was

[TABLE]

In the pilot study, 5 individuals out of 36 had missing values (missing at random) from at least one periods. To draw $(\boldsymbol{\beta},\Sigma)$ from their posterior distribution, we used the imputation-posterior (IP) algorithm and the Gibbs sampler by imputing $\mathbf{y}_{\mbox{\small{mis}}}$ given $\Sigma$ and $\boldsymbol{\beta}$ (the imputation step in Eqn (9)) and drawing $\Sigma$ and $\boldsymbol{\beta}$ respectively from their full conditional posterior distributions given $\mathbf{y}_{\mbox{\small{obs}}}$ and the imputed $\mathbf{y}_{\mbox{\small{mis}}}$ (the posterior step in Eqns (10) and (11)). The technical details on the derivation of the equations are provided in the Appendix.

[TABLE]

$M_{np\times np}$ in Eqn (9) was a diagonal indicator matrix of missingness on each observation with $m_{jj}=1$ if the $j\textsuperscript{th}$ response was missing and 0 otherwise for $j=1,\ldots,np$; $\sum_{i=1}^{n}\mathbf{e}_{i}\mathbf{e}^{\prime}_{i}=\sum_{i=1}^{n}(\mathbf{y}_{i}-\mathbf{x}_{i}\boldsymbol{\beta})(\mathbf{y}_{i}-\mathbf{x}_{i}\boldsymbol{\beta})^{\prime}$ in Eqn (10) was the scale matrix of the inverse Wishart distribution and $2\nu_{0}+n-p-1$ was the degree of freedom; $\mathbf{y}=(\mathbf{y}_{\mbox{\small{obs}}},\mathbf{y}_{\mbox{\small{mis}}})$ and $\Omega$ was the block diagonal matrix with $n$ blocks of $\Sigma$ in Eqn (11). Upon the convergence of the IP algorithm, after the burn-in and thinning periods, 1000 sets of posterior samples of $(\Sigma,\boldsymbol{\beta})$ were kept, based on which, 1000 power values were calculated analytically for the future trial via a multivariate $t$ distribution that incorporated the dependency structure among the 4 endpoints (refer to the online supplementary materials for the R function on the power calculation). The average of the 1000 power values yielded the model-based BEP.

The BEP assessment via the BBS and the BS2 techniques was much more straightforward than the model-based approach. Not only they involved minimal analytical work, coding of the procedures was also easy. After the steps in Table 2 and Table 3 ($m=1000$ and $t=1$), we tested the 4 equivalence hypotheses via a lme model by endpoint in each of the simulated trials via the BBS and BS2 procedures. If the 95% Cls of the log treatment differences for all 4 endpoints fell within $(\log(0.8),\log(1.25))$, the future trial was claimed a success. The success rate out of the 1000 simulated trials was the BEP.

We also computed the classical power analytically and via the bootstrap approach (Table 4). In the model-based approach, the log-scaled treatment differences in the 4 endpoints and the variance/covariance structure $\Sigma_{12\times 12}$ were estimated from the pilot study by fitting a lme model. We modelled $\Sigma$ with two different structures to examine the sensitivity of the model-based power to the assumed model: a UN $\Sigma$ and a Kronecker product $\Sigma=\Sigma_{1}\bigotimes\Sigma_{2}$ between a UN $\Sigma_{1}$ across the 4 endpoints and a CS $\Sigma_{2}$ across the 3 periods on the same endpoint. The estimated treatment differences and $\Sigma$ were then fed into the multivariate $t$-based power function to obtain the power. In the bootstrap approach, 1000 sets were generated and the hypothesis testing in each simulated $\tilde{\mathbf{y}}$ was the same as in the BBS and BS2 procedures. The success rate out of the 1000 trials was the bootstrapped power.

Though the pilot study is a 3-period crossover study while the future study is 2-period, the model-based power and BEP values calculated above were legitimate since the treatment differences were adjusted for the period effects through the lme models fitted to the pilot data (model-based) or to the future data (bootstrap-based); and it was not necessary to produce the “exact” 2-period crossover data in the FTS approach. In other words, the design inconsistency between the pilot and future trials was not a concern from the perspective of the treatment comparison as along as the period effects were properly taken care of.

The results are presented in Table 5 and summarized as follows. First, the BEP values via the BBS and BS2 procedures were similar. Second, the model-based BEP and power values were notably sensitive to the underlying parametric assumptions. The assumptions on the dependency structure $\Sigma$ among the endpoints influenced the estimates on the treatment differences in the pilot data, leading different power and BEP values under each assumption. Specifically, for BEP, different priors on $\Sigma$ led to a $\sim\!25\%$ difference; in the case of power, different structural assumptions on $\Sigma$ led to a $25\%\!\sim\!35\%$ difference. Third, the parametric and bootstrapped POS estimates were different in most cases. Fourth, T1 had consistently higher BEP than T2 though the absolute values varied by the assessment approach. The relative POS of T1 and T2 measured by the classical power were slightly inconsistent depending on the assessment approach. All taken together, we would recommend T1 for the future trial and cite the BEP obtained via the BBS and BS2 procedures given their robustness (model-free incorporation of the uncertainty around the parameters in the power calculation).

4.2 example 2: joint modeling of survival time and longitudinal data

The pilot study is a randomized trial in 40 AIDS patients who had failed or were intolerant of zidovudine (AZT) therapy (simulated from a published data set in Goldman et al., (1996)). The patients were allocated to two antiretroviral drug treatment groups (A and B) in a 1:1 ratio. Both the time to death due to AIDS and the CD4 count at 4 times points (0, 2, 6 and 12 months) were collected in the pilot study. Each patient was followed up for a minimum of 12 months and with an average of 15.6 months. Based on the results from the pilot study, a larger study of size $\tilde{n}$ is under planning to compare the two treatments in a more confirmatory manner and to test the hazard ratio of death between the two against 1. The team is interested in the assessment of the POS of the future trial.

In the model-based assessment of the BEP, the Bayesian joint modelling of the survival time (with right censoring) and the CD4 counts was applied. Specifically, the square root of the CD4 count $y_{ij}$ in subject $i$ at time $t_{j}$ was modeled with a lme model after the square-root transformation; and the Cox model was applied to analyze the survival time with the conditional expected mean of $y_{ij}$ as a covariate, plus the effect from the treatments (Eqn 12).

[TABLE]

$\mathbf{x}_{ij}(t)$ is the covariates for patient $i$ at time $t_{j}$ that includes an intercept, the time when $y_{ij}$ was measured ($x_{ij1}(t)$), the dummy variables for treatment $x_{i2}$, whether there was previous AIDS diagnosis at the study entry ($x_{i3}$), and whether the patient had AZT intolerance or AZT failure ($x_{i4}$); $\boldsymbol{\beta}$ contains the fixed effects associated with $\mathbf{x}_{ij}(t)$; the random effect $\mathbf{b}_{i}=(b_{i0},b_{i1})^{\prime}\sim N(\mathbf{0},\mathbf{G}_{2\times 2})$ corresponds a random intercept and a random slope for $x_{ij1}(t)$; and error term $\epsilon_{ij}\sim N(0,\sigma^{2})$. $\lambda$ in Eqn (12) quantifies the association between the CD4 count up to time $t$ and the hazard for death at $t$, and $\gamma$ is the log hazard ratio for death between the two treatments. The baseline hazard $h_{0}(t)$ was modelled with the penalized-splines approach, where $B_{l}(t,\mathbf{k})$ is the $l\textsuperscript{th}$ basis function of the splines at knots $\mathbf{k}=(k_{1},\ldots,k_{L})$ (we set $L=8$) and $\gamma_{l}$’s are the spline coefficients. The R package JMbayes was employed to fit the joint model in Eqn (12). After the convergence of the MCMC algorithm, 2000 posterior samples on the model parameters, including $\mathbf{b}_{i}$, were obtained (after 1000 burn-in and 40 thinning), from which the posterior probability of survival $p_{i}=\Pr(t_{i}\leq T)$ by some time $T$ in subject $i$ was obtained for $i=1,\ldots,40$. We examined two types of $T$: the last recorded time in the pilot study, which was 20.87 and 20.27 months in Treatments A and B, respectively; and 24 months, which was an extrapolation of the pilot data in study duration. The model-based power for testing $H_{0}:$ HR $=1$ is $\Phi\!\left(\!\sqrt{\tilde{n}(1-p)\gamma^{2}}/2-\Phi^{-1}(1-\alpha/2)\right),$ where $\tilde{n}(1-p)$ is the expected number of events in the future trial, where $p=n^{-1}\sum_{i=1}^{n}p_{i}$, and $\alpha$ is the type-I error rate. Power was calculated at each of the 2000 posterior samples of ($p,\gamma)$, the average of which gave the model-based BEP. The model-based power given the posterior means of $\gamma$ and $p$ at two types of $T$ was also calculated.

In the bootstrap-based approaches, $m=500$ sets of future data were generated according to the steps in Table 1 and Table 3 ($t=1$). In each simulated trial, the joint model in Eqn (12) was applied, and the 95% posterior interval on the HR was obtained. If the lower bound of the interval was $>1$ or the upper bound was $<1$, then the future trial was claimed a success; the success rate out of 500 was the BEP. The bootstrap-based power based on Table 3 was also calculated. Different from the model-based approach, extrapolation to a longer study duration ($T=24$) was not possible without further parametric assumptions. In other words, the BEP and power evaluated via the bootstrap-based approaches implicitly assumed the study duration in the future trial was the same as that in the pilot study.

The results of POS at different $\tilde{n}$ are provided in Table 6. The model-based BEP and the BEP via the BBS and BS2 procedures were similar in this example, so were between the model-based power and the bootstrap power. In summary, the POS for the future trial reached $80\sim 85\%$ at $\tilde{n}=400$, which might be an acceptable level of POS. As expected, a longer study duration (24 months) led to an increase in the POS (by $5\%\sim 10\%$) compared to the same sized future trial with the same study duration as the pilot study. By contrast, all the POS values assessed by the power were close to 100% at all $\tilde{n}$ and for both study durations, thus not a differentiable nor a valuable metric in determining $\tilde{n}$ .

5 Discussion

We developed the BBS procedure based on the Bayesian bootstrap to calculate the BEP given individual-level pilot data. We also presented a non-Bayesian counterpart to the BBS procedure, named the double bootstrap (BS2), that achieves the same goal as the BBS for the BEP assessment. Neither procedures make assumptions about the distribution of the pilot data and can handle multidimensional data sets without imposing a dependence structure among the variables. The implementation of both procedures are straightforward: only a few lines of codes are needed and the whole procedures can be easily parallelized for fast computation. By contrast, the parametric approaches specify a likelihood function on the pilot data set and priors for the model parameters, followed by posterior sampling of the parameters involved in the analytical power calculation. The derivation of the posterior distributions and the posterior sampling can be complicated and computational costly, not to mention the possibility of model mis-specification on the pilot data.

In the case of existence of multiple sets of pilot data, the data sets can be easily combined to simulate the future trial. The weight associated with each pilot study, in terms of their contribution to the future data, by default would be the size of the pilot studies: the larger the sample size of a pilot study, the more likely the subjects from that study will be bootstrapped. If other weights are desired and specified (e.g., the design agreement between the pilot and future trials, the quality of the pilot data), they can be conveniently incorporated using a 2-stage sampling procedure: first sampling the studies with probabilities proportional to the corresponding weights and followed by the regular BBS and BS2 procedure within each study.

The bootstrap-based approaches require the pilot information to be available in the form of individual-level data. When there exist only historical summary/aggregate statistics, the bootstrap-based approaches are not directly applicable. Though future data can be simulated from the aggregate statistics, it would require additional distributional assumptions, defeating the purposes of developing the nonparametric bootstrap methods in the first place. In practice, the pilot study might not be perfectly matched up to simulate the future trial such as the design discrepancy between the two. Some design differences may not have direct impacts on the planned analysis on the future trial or can be adjusted for in the BEP calculation, and thus will not be a concern (e.g., the first case study). If a discrepancy that relates to the effect size in the power calculation is not easy to adjust for without making further parametric assumptions (such as the second example where the study duration could be different), the model-based approaches might be the only choice. However, if unreasonable extrapolation has to made, the parametric approaches would not be appropriate either.

In summary, the BBS and BS2 approaches provide an alternative to the parametric approaches to assess the BEP. The bootstrap-based procedures will appeal to non-Bayesian practitioners given their analytical and computational simplicity and easiness in implementation. We provide the sample R functions on the BBS and BS2 procedures in the online supplementary materials to facilitate their practical implementation.

Supplementary materials

The supplementary materials can be found at http://www3.nd.edu/~fliu2/BBS-supp.pdf

Appendix A: Conditional posterior distributions of $\Sigma$ and $\boldsymbol{\beta}$ in case study 1

This appendix presents the derivation of the posterior distribution of $\Sigma$ and $\boldsymbol{\beta}$ in linear mixed models. Denote the number of subjects by $n$ and the number of measurements per subject by $p$. The likelihood of $\Sigma$ and $\boldsymbol{\beta}$ is $L(\Sigma,\boldsymbol{\beta}|\mathbf{y})\!\propto\!\textstyle\prod_{i=1}^{n}\!\left\{\!|\Sigma|^{-1/2}\mbox{exp}\left(-\frac{1}{2}(\mathbf{y}_{i}-\mathbf{x}_{i}\boldsymbol{\beta})^{T}\Sigma^{-1}(\mathbf{y}_{i}-\mathbf{x}_{i}\boldsymbol{\beta})\right)\!\right\}$. Let the prior on $\boldsymbol{\beta}$ and $\Sigma$ be $p(\boldsymbol{\beta},\Sigma)\!=\!|\Sigma|^{-(\nu_{0}+p+1)/2}$. The full conditional posterior distribution of $\Sigma$ given $(\boldsymbol{\beta},\mathbf{y})$ is

[TABLE]

which is an inverse Wishart distribution $W^{-1}(\sum_{i=1}^{n}\mathbf{e}_{i}\mathbf{e}_{i}^{t},n)$ with scale matrix $\sum_{i=1}^{n}\mathbf{e}_{i}\mathbf{e}_{i}^{t}$ and the degrees of freedom $n$. The full conditional distribution of $\boldsymbol{\beta}$ given $(\Sigma,\mathbf{y})$ is

[TABLE]

[TABLE]

Thus $p(\boldsymbol{\beta}|\boldsymbol{\Sigma},\mathbf{y})=\mbox{N}\left((\mathbf{x}^{T}\Omega^{-1}\mathbf{x})^{-1}(\mathbf{x}^{T}\Omega^{-1}\mathbf{y}),(\mathbf{x}^{T}\Omega^{-1}\mathbf{x})^{-1}\right)$, where $\Omega$ is a block diagonal matrix with $n$ blocks of $\Sigma$ on the diagonal, $\mathbf{x}=\{\mathbf{x}_{i}\}_{i=1,\ldots,n}$, and $\mathbf{y}=\{\mathbf{y}_{i}\}_{i=1,\ldots,n}$