Calculation of the probability of success with the package probSuccess
Paul Jordan
Last update: 12 October, 2021
- 1 probSuccess
- 2 Installation
- 3 Types of functions contained in probSuccess
- 3.1 Distribution functions
- 3.2 Normal conjugate analysis
- 3.3 Posterior distributions with arbitrary priors
- 3.4 Probability of success - external information
- 3.5 Probability of success after an interim analysis
- 3.6 Probability of success for survival analyses
- 3.7 Conditional and predictive power
- 3.8 Power of a logrank test with one futility analysis
- 3.9 Miscellaneous
- 4 Examples for probability of success (DDCP)
- 5 Examples for conditional and predictive power
- 6 Examples for Bayesian methods
- 7 Theoretical background
- 7.1 Probability of Success after Interim Analysis
- 7.2 Conditional Normal Distributions at an Interim Analysis
- 7.3 Conjugate Bayesian Estimation of the Full Normal Model
- 7.4 Gamma-related densities
- 7.4.1 \(\Gamma\) Distribution
- 7.4.2 \(\text{Inv}\!-\!\Gamma\) Distribution
dInvGamma,pInvGamma,rInvGamma - 7.4.3 \(\text{Scaled}\!-\!\chi^2\) distribution
dScaledChisq,pScaledChisq,rScaledChisq - 7.4.4 \(\text{Scaled-inv}\!-\!\chi^2\) distribution
dScaledIChisq,pScaledIChisq,rScaledIChisq - 7.4.5 Relation between Distributions
- References
1 probSuccess
This R package provides functions for calculating the probability of success (or DDCP). The precursor was called PTS for Probability of Technical Success. As this terminology can be misleading and ambiguous, we changed it to the simpler Probability of Success. Nevertheless there are still function names that kept this wording, such as PTSafterStudy.
Furthermore, the package provides auxiliary functions for Bayesian calculations, such as for the full conjugate analysis of th normal model with unknown mean and variance.
For theoretical background, we refer to the last section ‘Theoretical background’, the references cited therein or to (Rufibach, Jordan, and Abt 2014).
2 Installation
2.1 Stable Release
3 Types of functions contained in probSuccess
3.1 Distribution functions
Specific distributions (PDF, CDF, random numbers).
Gamma_Related_DensitiesScaledBetaBetaBinom
3.2 Normal conjugate analysis
Posterior distributions and predictive posterior distribution for the normal model.
Full_Conjugate_Normal_ModelpostMeanConjPredNormNoninf
3.3 Posterior distributions with arbitrary priors
Posterior distributions for the univariate normal model allowing for arbitrary prior distriburions.
updatePriorMeanPredDist
3.4 Probability of success - external information
Probability of succes for univariate or multivariate normal endpoints allowing for arbitrary prior distributions (“…MV” in the function name stands for ‘multivariate)’.
PTSafterStudyPTSafterStudyMVProbGTcut
3.5 Probability of success after an interim analysis
Probability of succes after an interim analysis for univariate or multivariate normal endpoints allowing for arbitrary prior distributions. ProbSuccess_IA and PTSafterIA deal with an univariate normal endpoint. In case of a noniformative prior, ProbSuccess_IA is more efficient! PTSafterIAMV covers the case of a multivariate normal endpoint with arbitrary prior (“…MV” in the function name stands for ‘multivariate)’.
ProbSuccess_IAPTSafterIAPTSafterIAMV
3.5.1 Differences between PTSafterIA and ProbSuccess_IA
In PTSafterIA, the variance \(\sigma^2\) of the endpoint is assumed to be known and the observed standard deviation(s) is/are imputed. In ProbSuccess_IA, the variance of the endpoint is unknown (see table below).
| Function | Variance | Prior |
|---|---|---|
PTSafterIA |
known | arbitrary |
ProbSuccess_IA |
unknown | no prior |
3.6 Probability of success for survival analyses
Calculation of the probability of success based on the logarithm of the hazard ratio log(HR).
PSsurvival
3.7 Conditional and predictive power
Conditional and Bayesian predictive power for normal endpoints, binary endpoints and time-to-event endpoints. * CondPower * PredPower
The functions are based on the unified approach for adaptive designs using the concept of Brownian motion. PredPower.defaultuses the notations used in this concept; for technical details see the references (Proschan, Gordon Lan, K.K., and Turk Wittes 2006). The functions PredPowerMeans, PredPowerRates and PredPowerTE are wrappers which transform the design-specific parameters to the parameters used in PredPower.default.
By setting \(\sigma_{0} = 0\), the predictive power reduces to the conditional power under the originally hypothesized treatment effect \(\theta_{0}\). By choosing a very large variability, the influence of the prior is negligible and the Bayesian Predictive power becomes an integration of the conditional power over the posterior density of the study effect at the interim analysis. By determinating the prior variance, one can choose how much weight one gives the prior treatment estimate relative to the estimate from the trial at the interim analysis. If the weight of the prior treatment is \(w, \quad(0 \le w \le 1)\), it follows that \(\sigma_0^2= (1-w)/w\).
3.8 Power of a logrank test with one futility analysis
Computes the power of a logrank test accounting for one futility analysis.
powerfun
3.9 Miscellaneous
Miscellaneous utility function.
CovCondMkCovsearchPriorVarPredCondIAmv
4 Examples for probability of success (DDCP)
4.1 External information, probability of succsess for a continuous endpoint
4.1.1 Basic example (fictitious)
Suppose that there are two comparable studies on the same continuous endpoint with the following results:
- Means (difference in treatment arms): \(\bar{x}_1 =12.3,\bar{x}_2= 11.4\).
- Pooled standard deviations: \(\text{sd}_1=3.7,\text{sd}_2=4.4\).
- Sample sizes (per arm): \(n_1=48, n_2=62\).
A new study is planned based on:
- Sample size: 110 per arm.
- Expected standard deviation (in either arm): 4.0 .
- Criterion for succes : Treatment difference > 12 .
For the true treatment difference \(\theta\) the following almost uninformative prior is chosen: \(\theta \sim \mathcal{N}(10,10^6)\).
library(probSuccess)
library(tidyverse)
pSucc <- PTSafterStudy(TPP=12,Mean=c(12.3,11.4),SD=c(3.7,4.4),N=c(48,62)
,nNew=110,sdNew=4.0,diffRef=FALSE,nStudy=1,nSucc=1
,prior=dnorm,mean=10,sd=1000)
print(pSucc)## [1] 0.40560754.1.2 A more complicated example
The table below summarises the results from three previous studies:
| Study | Nact | Nplac | TrtDiff | SE |
|---|---|---|---|---|
| 1 | 52.00 | 49.00 | -0.83 | 0.75 |
| 2 | 65.00 | 81.00 | -2.89 | 0.71 |
| 3 | 53.00 | 67.00 | -0.10 | 0.92 |
with:
Nact: Number of Patients in the active armNplac: Number of Patients in the placebo armTrtDiff: Point estimator of treatment differenceSE: Standard Error ofTrtDif
One or more new (Phase III) studies with 100 patients per arm are planned. The success criterion is: Effect size > 0.35
An almost uninformative normal prior for the the effect size was used: \(\mathcal{N}(0,1000)\).
In order to fit this into the parametrization, the results of the three previous studies are transformed to effect sizes:
sd1 <- 0.75/sqrt(1/49+1/52)
es1 <- 0.83/sd1
sd2 <- 0.71/sqrt(1/81+1/65)
es2 <- 2.89/sd2
sd3 <- 0.92/sqrt(1/67+1/53)
es3 <- 0.1/sd3As diffRef = TRUE, the means are differences to a reference group. Therefore, the (pooled) standard deviation per arm in the tree sudies is set to:
SD1 <- sqrt(0.5*(1/49+1/52))
SD2 <- sqrt(0.5*(1/81+1/65))
SD3 <- sqrt(0.5*(1/67+1/53))The sample sizes are then set to 1 !
PoS for one study out of one:
PTSafterStudy(TPP=0.35,Mean=c(es1,es2,es3),SD=c(SD1,SD2,SD3),N=c(1,1,1)
,nNew=100,sdNew=1,diffRef=TRUE,nStudy=1,nSucc=1)## [1] 0.4697344PoS for one study out of two new studies:
PTSafterStudy(TPP=0.35,Mean=c(es1,es2,es3),SD=c(SD1,SD2,SD3),N=c(1,1,1)
,nNew=100,sdNew=1,diffRef=TRUE,nStudy=2,nSucc=1)## [1] 0.661358PoS for two studies out of two new studies:
PTSafterStudy(TPP=0.35,Mean=c(es1,es2,es3),SD=c(SD1,SD2,SD3),N=c(1,1,1)
,nNew=100,sdNew=1,diffRef=TRUE,nStudy=2,nSucc=2)## [1] 0.27811074.2 Internal information: probability of success after an interim analysis
Let us assume that in a interim analysis (IA) for a two-arm study (active drug vs. placebo), the following values are known or reported:
- Sample-size at interim analysis: \(n_{IA}=20\), for both arms.
- Final sample-size: \(n_{fin}=40\), for both arms.
- Continuous endpoint: \(\theta = \theta_{act} - \theta_{plb}\).
- Reported pooled standard deviation for the mean difference to placebo: \(sd_{IA}=4\) (
diffRef=TRUE). - Critical value (success criterion, cutpoint, TPP, …): \(10\).
nIA <- 20
nFinal <- 40
Crit <- 104.2.1 Mean difference to placebo reported
Let us assume that addditionally:
- Observed at IA:
- Mean difference of 10.4 (= estimator for \(\theta\)).
- sd of 4.1 (standard deviation, not standard error of the mean!).
meanIA <- 10.4
sdIA <- 4.14.2.1.1 Using PTSafterIA
Moreover, for \(\theta\) the following prior distribution has been chosen: \(\mathcal{N}(\mu=11,\sigma^2=10^2)\) ( Ifprior is missing, an uninformative (improper) prior is assumed).
PTSafterIA(TPP=Crit,type="value",MeanIA=meanIA,IntIA=NULL,sdIA=sdIA,nIA=nIA,nFinal=nFinal
,diffRef=TRUE,prior=dnorm,mean=11,sd=10)## [1] 0.69415334.2.1.2 Using ProbSuccess_IA
ProbSuccess_IA(crit=Crit,meanIA=meanIA,sdIA=sdIA,nIA=nIA,nFinal=nFinal)## [1] 0.72367244.2.2 Only interval containing the mean difference is known at IA
Let us assume that at the IA it is only known, that the study will be continued. I.e. it was stopped either for futility or early efficacy.
4.2.2.1 Two-sided interval
- Futility boundary: 8.5
- Early efficacy boundary: 12
- Prior on \(\theta: p_0(\theta) =\mathcal{N}(\mu=10,\sigma^2=20^2)\).
PTSafterIA(TPP=10,type="interval",IntIA=c(8.5,12),sdIA=sdIA,nIA=nIA,nFinal=nFinal
,diffRef=TRUE,prior=dnorm,mean=10,sd=20)## [1] 0.56700424.2.2.2 One-sided interval
- No stopping for futility.
- Early efficacy boundary: 12
- Prior on \(\theta: p_0(\theta) =\mathcal{N}(\mu=10,\sigma^2=5^2)\).
PTSafterIA(TPP=10,type="interval",IntIA=c(-Inf,12),sdIA=sdIA,nIA=nIA,nFinal=nFinal
,diffRef=TRUE,prior=dnorm,mean=10,sd=5)## [1] 0.2321571Note: If the data are given the form of an interval containing the mean, -Inf and Inf are not allowed when no proper prior is given! We recomend that the variance of the prior should not be greater than ca. 25 times the sd at IA. The conditional distribution of the final result given the interim data for standardized variables can be derived via unified distribution theory (Jennison and Turnbull 2000) (Rufibach, Jordan, and Abt 2014).
4.3 Multivariate continuous endpoints
4.3.1 External information, 2-dimensional endpoint, without prior
The endpoint is represented as a 2-dim vctor: \(\vec{y}=(y_1,y_2)\).
- Critical value: \(\vec{y}_c = (12, 11)\)
- Endpoint: Difference to placebo.
- Sample-size of first study: \(N_{old}=30\)
- Sample-size of planned new study: \(N_{new}=80\)
- Observed mean in a previous stud : \((\bar{y_1}, \bar{y_2}) =(11.6,11.3)\).
- Observed covariance matrix: \(\mathbf{S}= \begin{bmatrix} \mathbf{9.3} & \mathbf{4.2}\\ \mathbf{4.2} & \mathbf{15.8}\\ \end{bmatrix}\)
COV <- matrix(c(9.3,4.2,4.2,15.8),ncol=2)
M <- c(11.6,11.3)
PTSafterStudyMV(TPP=c(12,11),Mean=M,N=30,nNew=80,Cov=COV,diffRef=TRUE)## [1] 0.2473748Adding a prior:
Let us assume that our previous knowledge on the 2-dimesional endpoint is represented by:
\[\mathcal{N}\left( \mu= \begin{bmatrix} 10\\ 10 \end{bmatrix}, \Sigma= \begin{bmatrix} \mathbf{20} & \mathbf{0}\\ \mathbf{0} & \mathbf{20}\\ \end{bmatrix} \right)\]
Sigma <- 20*diag(2)
PTSafterStudyMV(TPP=c(12,11),Mean=M,N=30,nNew=80,Cov=COV,diffRef=TRUE,
,prior=dmvnorm,mean=c(10,10),sigma=Sigma)## [1] 0.22082534.3.2 Internal information
Let us consider the following example:
- Bivariate endpoint at the interim analysis (20 subjects).
- The means are differences to a reference group (
diffRef=TRUE). - the mean for endpoint I is \(11.2\) and the mean for endpoint II \(11.9\).
- The covariance matrix of the observations is:\(\mathbf{S}= \begin{bmatrix} \mathbf{9.3} & \mathbf{4.2}\\ \mathbf{4.2} & \mathbf{15.8}\\ \end{bmatrix}\)
- The cutpoint for success is \(11\) for either endpoint.
- Final analysis after 55 subjects.
- Prior on the endpoint: \(\mathcal{N}\left( \mu= \begin{bmatrix} 10\\ 10 \end{bmatrix}, \Sigma= \begin{bmatrix} \mathbf{50} & \mathbf{0}\\ \mathbf{0} & \mathbf{50}\\ \end{bmatrix} \right)\)
MeanIA <- c(11.2,11.9)
nIA <- 20
nFinal <- 55
CovIA <- matrix(c(9.2,3.1,3.1,15.8),ncol=2)
TPP <- c(11,11)
PTSafterIAMV(TPP,type="value",MeanIA=MeanIA,IntIA=NULL,CovIA=CovIA,nIA=nIA
,nFinal=nFinal,diffRef=TRUE
,prior=dmvnorm,mean=c(10,10),sigma=50*diag(2))## [1] 0.5061063Now let us assume, that at the IA, there is only known (other specification equal):
- Mean I lies in the interval [10,12]
- Mean II lies in the interval [9,13]
IntIA <- IntIA <- rbind(c(10,12),c(9,13))
PTSafterIAMV(TPP=c(11,11),type="interval",MeanIA=NULL,IntIA=IntIA
,CovIA=CovIA,nIA=nIA,nFinal=nFinal,diffRef=TRUE,prior=dmvnorm
,mean=c(10,10),sigma=50*diag(2))## [1] 0.25744584.4 Survival analyses
The probability of success is defined as the probability to observe a hazard ratio less than a critical value in a trial with n events. The previous knowledge about the true hazard ratio \(\lambda\) is expressed in the original scale by a lognormal distribution with mean \(\mu\) and variance \(\sigma^2\).
PSsurvival(Mlambda=0.9,Slambda=0.04,Crit=0.814,n=630)## [1] 0.14023065 Examples for conditional and predictive power
The functions are based on the unified approach using the appropriate z-statistic (Proschan, Gordon Lan, K.K., and Turk Wittes 2006). More details are given in the documentationto probSuccess.
5.1 Conditional power: Comparing two means
The example is taken from (Dmitrienko and Wang 2006) (page 2187). In the paper, alpha is one-sided.
meanT <- 9.2
sdT <- 7.3
n1T <- 11
n2T <- 41
meanC <- 8.4
sdC <- 6.4
n1C <- 10
n2C <- 41
CondPowerMeans(meanT=meanT, meanC=meanC, sdT=sdT, sdC=sdC, n1T=n1T, n1C=n1C, n2T=n2T, n2C=n2C,
delta = 4, sigma=8.5, alpha = 0.2)## [1] 0.69484435.2 Conditional power: Comparing two proportions
The example is taken from (Proschan, Gordon Lan, K.K., and Turk Wittes 2006), Example 3.3. (page 50ff).
CondPowerRates(kT=24,kC=22,n1T=48,n1C=44,n2T=200,n2C=200,PT=0.45,PC=0.6)## [1] 0.65673765.3 Bayesian predictive power: Comparing two means
The example is taken from (Proschan, Gordon Lan, K.K., and Turk Wittes 2006), Example 3.2.
delta <- 2 ### expected difference between treatment and control
sigma <- 5 ### expected SD of difference
n1T <- 52
n1C <- 50
n2T <- 132
n2C <- 132
diff <- 1 ### observed treatment difference at interim analysis
s_pool <- 6.1 ### observed pooled standard deviationThe prior information is given a weight of ca. 10% compared to data at the end of the trial; therefore: \[\frac{1/\sigma_0^2} {1/\sigma_0^2 + 1} = 0.1 \rightarrow \sigma_0^2 = 9\]
sigma.0 <- 3
## Type I error one-sided: 0.025
PredPowerMeans(meanT=diff, meanC=0, sdT=s_pool, sdC=s_pool, n1T=n1T, n1C=n1C, n2T=n2T, n2C=n2C
,delta.0 = delta, sigma=sigma ,sigma.0=sigma.0,alpha = 0.025)## [1] 0.3775832For \(\sigma_0^2 = 0\), the predictive power coincides with the conditional power:
PP <- PredPowerMeans(meanT=diff, meanC=0, sdT=s_pool, sdC=s_pool, n1T=n1T, n1C=n1C, n2T=n2T, n2C=n2C
,delta.0 = delta, sigma=sigma, sigma.0=0, alpha = 0.025)
print(PP)## [1] 0.7582574CP <- CondPowerMeans(meanT=diff, meanC=0, sdT=s_pool, sdC=s_pool, n1T=n1T, n1C=n1C, n2T=n2T, n2C=n2C
,delta = delta, sigma=sigma, alpha = 0.05) ### two-sided alpha!
print(CP)## [1] 0.75825746 Examples for Bayesian methods
6.1 Updating a prior distribution by the result of a previous study
Data from a previous pilot study (continuous and normally distributed endpoint):
- Observed mean: \(\bar{x}=11\)
- Observed standard deviation: \(\text{sd}(x)=4\)
- Number of observations: 8
As prior for the population mean \(\mu\), a generalized t-distribution with location \(\mu=8\), scale \(\sigma=2\) and \(\nu=1\) degrees of freedom is chosen:
Mean <- 11
SD <- 4
N <- 8
m_p <- 8
sd_p <- 2
nu_p <- 1
mu <- seq(4,15,0.1)
prior <- dT(mu,nu_p,m_p,sd_p)
post <- updatePriorMean(mu,type="value",Mean=Mean,SD=SD,Int=NULL,N=N
,diffRef=FALSE,prior=dT,nu=nu_p,m=m_p,s=sd_p)
Density <- c(prior,post)
Type <- rep(c("prior","posterior"),each=length(mu))
mu <- rep(mu,2)
PostDens <- data.frame(mu,Density,Type)
plotPost <- ggplot(PostDens,aes(x=mu,y=Density,col=Type)) + geom_line() +
geom_vline(xintercept=Mean,linetype="dashed",col="blue") +
annotate("text",x=Mean,y=0.025,label = expression(bar(x)==11),
hjust=0,col="blue")
print(plotPost)
If instead of the mean, only an interval containing the mean is known: \(\bar{x} \in [0,9]\):
mu <- seq(-5,15,0.1)
prior <- dT(mu,nu_p,m_p,sd_p)
post<- updatePriorMean(mu,type="interval",SD=2,Int=c(0,9),N=16
,diffRef=TRUE,prior=dT,nu=2,m=11,s=10)
Density <- c(prior,post)
Type <- rep(c("prior","posterior"),each=length(mu))
mu <- rep(mu,2)
PostDens <- data.frame(mu,Density,Type)
plotPost <- ggplot(PostDens,aes(x=mu,y=Density,col=Type)) + geom_line() +
geom_vline(xintercept=Mean,linetype="dashed",col="blue")
print(plotPost)
6.2 Posterior predictive distribution
6.2.1 Continuous (normal) endpoint
In a pilot study with \(n_1=30\) observations \(x_1,\ldots,x_{n_1}\) a mean \(\bar{x_1}=10.3\) and a standard deviation \(s^2_1=3.8\) were observed.
The function dPred provides the density of the predictive distributions for the observed mean \(\bar{x_2}\) of \(n_2=45\) new observations, using a noniformative prior.
A success criterion is given by: \(\bar{x_2} \ge x_{crit}\), with \(x_{crit}=9.5\).
Then, the probabilty of success is given by the integral over the predictive distribution: \[\int_{-\infty}^{crit} p(\overline{x}_2|\overline{x}_1,s^2_1)dx\] This interval can be calculated by means of the function evalPred.
#### Old observations:
m1 <- 10.3
SD <- 3.8
n1 <- 30
#### Posterior predictive distribution for the mean of 20 new observations
Nnew <- 45
crit <- 9.5
mx <- seq(5,15,0.1)
predDist <- dPred(mx,m1=m1,sd1=SD,n1=n1,n2=Nnew,VarKnown=TRUE)
Dist <- data.frame(mx,predDist)
### probability of succes
pSucc <- evalPred(m1=m1,sd1=SD,n1=n1,n2=Nnew,crit=crit,VarKnown=TRUE)
Dist2 <- Dist %>% filter(mx >= crit)
PlotPred <- ggplot(Dist,aes(x=mx,y=predDist)) + geom_line(col="blue") +
xlab(paste("Observed mean of",Nnew,"new observations")) + ylab("Posterior predictive density")+
geom_vline(xintercept=crit,linetype="dashed",col="red") +
geom_area(data=Dist2,aes(x=mx,y=predDist),fill="red",inherit.aes = FALSE,alpha=0.3)+
annotate("text",x=9.8,y=0.05,label = paste("P(success)=",round(pSucc,3)),
hjust=0,col="red")
print(PlotPred)
cat("\n prob. Success: ",round(pSucc,3),"\n")##
## prob. Success: 0.814Note: Using the function PredDist,an arbitrary prior for the mean can be chosen. However, it much slower tdue to iterated numeric integrations.
6.2.2 Binary endpoint (conjugate analysis)
The posterior predictive distribution for the binomial model with a Beta-prior \(\mathbf{Beta}(\alpha,\beta)\) on the parameter \(\theta\) is the beta-binomial distribution \(\text{Beta-bin} (k,n,\alpha,\beta)\) implemented int the function BetaBinom.
Example:
- Prior on mortality rate: \(\mathbf{Beta} (\alpha=3,\beta=27)\).
- Suppose we now operate on \(n = 10\) patients and observe \(k = 0\) deaths.
- What is the probability that:
- The next patient will survive the operation?: \(\text{Beta-bin} (0,1,\alpha=3+0,\beta=27+10)\)
- That there are 2 or more deaths in the next 20 operations? \(1- \sum_{i=2}^{20}\text{Beta-bin}(i,20,\alpha=3+0,\beta=27+10)\)
## Probability that the next patient will survive the operation
dbetabin(0,10,a=3+0,b=27+10)## [1] 0.4960378## Probability that here are 2 or more deaths in the next 20 operations
1-pbetabin(1,20,a=3+0,b=27+10)## [1] 0.4176755## Compare
1-(dbetabin(0,20,a=3+0,b=27+10)+dbetabin(1,20,a=3+0,b=27+10))## [1] 0.41767556.3 Finding parameters for the priors in the full normal model with searchPriorVar
A convenient conjugate prior for the full normal model (\(\mu\) and \(\sigma^2\) unknown), is given by (Gelman et al. 1995): \[ \mu|\sigma^2 \sim \mathcal{N}(\mu_0,\sigma^2/\kappa_0)\] \[ \sigma^2 \sim \text{Inv}\!-\!\chi^2(\nu_0,\sigma_0^2).\] Here, $ ^2 ^2(_0,_0^2)$ denotes the scaled inverse-\(\chi^2\) density as implemented in dScaledIChisq and pScaledIChisq (see Gamma_Related_Densities). The function searchPriorVar helps finding parameters for the prior of the variance in the full normal model (scaled inverse-chisq distribution).
First, one has to choose an interval reflecting the prior belief on \(\sigma^2\), i.e. an interval that contains the variance with a predefined coverage probability;
Example:
With 90% probability, the variance is between 9 and 49: \[Prob(9 \le sigma^2 \le 49) = 0.9\]
parm <- searchPriorVar(9,49,0.9)
parm ## df scale
## 8.357316 17.246252## Checking the result
x <- rScaledIChisq(100000,parm[1],parm[2]) ### Random numbers
quantile(x,c(0.05,0.95))## 5% 95%
## 9.013292 48.876680The degrees of freedom df corresponds to \(\nu_0\) and scale to \(\sigma_0^2\).
6.4 Sampling from the Sufficient Statistics of a Normal Distribution
For simulations it is far more efficient to sample directly the sufficient statistics for a given distribution. In the case of the normal model \[x \sim \mathcal{N}(\mu,\sigma^2)\] it is the sample mean and the sample variance \((\overline{x},s^2)\).
The sample mean is distributed as: \[\overline{x} \sim \mathcal{N}(\mu,\frac{\sigma^2}{n})\] where \(n\) denotes the sample size, and the ample variance:
\[s^2 \sim \chi^2_{n-1}\times \frac{\sigma^2}{n-1}\]
The parametrization was chosen in order to be consistent with Bayesian applications (Gelman et al. 1995).
Moments: \(\text{E}(X)=scale; \quad \text{Var}(X)=2\times scale^2/df\).
With this notation, one gets:
\[s^2 \sim \text{Scaled}\!-\!\chi^2(df=n-1,scale=\sigma^2).\] \(\text{Scaled}\!-\!\chi^2(\nu_0,\sigma_0^2)\) scaled inverse-\(\chi^2\) density as implemented in Gamma_Related_Densities.
mu <- 0
sigma <- 2
N <- 5
Nsim <- 100000 ### Number of simulations
xVar <- rchisq(Nsim,N-1)/(N-1)*sigma^2
hist(xVar,prob=TRUE,nclas=50,xlab="Sampled variance",main="")
curve(dchisq((N-1)*x/sigma^2,N-1)*(N-1)/sigma^2,ad=TRUE,col=3,lwd=3)
7 Theoretical background
7.1 Probability of Success after Interim Analysis
7.1.1 Mean at IA known
The probability of success (or DDCP) is defined as the probability that the mean of a random variable is greater than a critical value \(x_{crit}\), given the result of an interim analysis \((\sigma^2|\overline{x}_{ia};s_{ia}^2)\). It can be formally written as: \[\begin{equation} \label{PoS} \text{prob}(\overline{x}_{fin} \geq x_{crit}) = \int(1-\Phi(x_{crit}; \theta = \theta_c, \sigma^2 = \sigma_c^2)) \times p(\theta,\sigma^2|\overline{x}_{ia};s_{ia}^2)d\theta d\sigma^2 \quad (1) \end{equation}\]
In \((1)\), \(\Phi\) denotes the normal CDF and \(\theta_c\), \(\sigma_c^2\) (which depend on \(\theta\) and \(\sigma^2\)) are as defined in “Conditional Normal Distributions at an Interim Analysis”.
\(p(\theta,\sigma^2|\overline{x}_{ia};s_{ia}^2)\) is the posterior density of the parameters \(\theta\) and \(\sigma^2\), given \((\overline{x}_{ia};s_{ia}^2)\).
7.1.2 Interval for mean known at IA
If at the interim analysis only an interval for the mean is known \[\theta_0 \leq \overline{x}_{ia} \leq \theta_1,\] the conditional probability for \(\overline{x}_{fin}\) is given as
\[\begin{equation} \text{Prob}(\overline{x}_{fin} \geq x_{crit}| \overline{x}_{ia} \in [\theta_0,\theta_1], \theta, \sigma^2) = \frac{p(\overline{x}_{fin} \geq x_{crit}, \overline{x}_{ia} \in [\theta_0,\theta_1]|\theta,\sigma^2 )}{p(\overline{x}_{ia} \in [\theta_0,\theta_1]| \theta, \sigma^2)} \quad(2) \end{equation}\]
The joint distribution of \((\overline{x}_{ia},\overline{x}_{fin})\) is given in “Conditional Normal Distributions at an Interim Analysis”.
In order to obtain the PoS, \((2)\) has again to be integrated over the posterior distribution of the parameters as in \((1)\).
7.1.3 Multivariate normal mean
The multivariate case is analogous and can be easily derived; here we omit it for simplicity.
7.1.4 Relation to conditional and predictive power
The variance \(\sigma^2\) is often assumed to be known and the integration in\((1)\) is carried out over a posterior density for \(\theta\) alone. For \(\sigma^2\), \(s_{ia}^2\) is imputed.
With known \(\sigma^2\), the \(1-\alpha\) quantile of \(\mathcal{N}(\theta, \sigma^2)\) can be be chosen as \(x_{crit}\). In this case \((1)\) becomes the predictive power. If the values such as used for a ‘classical’ power calculation are chosen for \(\theta\) and \(\sigma^2\) in \(x_{crit}, (1)\) reduces to: \[\begin{equation} \label{CondPow1} \text{prob}(\overline{x}_{fin} \geq x_{crit}) = 1-\Phi(x_{crit}; \theta = \theta_c, \sigma^2 = \sigma_c^2) \end{equation}\]
7.2 Conditional Normal Distributions at an Interim Analysis
Let \(x \sim \mathcal{N}(\theta, \sigma^2)\) and \(\overline{x}_{ia}=\hat{\theta}_{ia}\) be the (observed) mean of \(n_{ia}\) observations at the time point of the interim analysis (IA) and \(s^2_{ia}\) the corresponding sample variance. Then the conditional density of the observed mean \(\overline{x}_{fin}\) at the final analysis with \(n_{fin}\) observations is given by: \[\overline{x}_{fin}|\overline{x}_{ia} \sim \mathcal{N}(\theta_c, \sigma_c^2)\] where
\[\theta_c=\frac{n_{ia}}{n_{fin}}\overline{x}_{ia}+(1-\frac{n_{ia}}{n_{fin}})\theta\] and
\[\sigma_c^2=\frac{\sigma^2}{n_{fin}} (1 - \frac{n_{ia}}{n_{fin}})\]
7.2.1 Multivariate case
The multivariate case is analogous; the conditional distribution of the observed multivariate mean \(\overline{\mathbf{x}}_{fin}\) given \(\overline{\mathbf{x}}_{ia}\) and the observed covariance matrix \(\Sigma_{ia}\) is:
\[ \overline{\mathbf{x}}_{fin}|\overline{\mathbf{x}}_{ia} \sim \mathcal{N}(\mathbf{\theta}_c, \Sigma_c)\]
where \[\mathbf{\theta}_c=\frac{n_{ia}}{n_{fin}} \overline{\mathbf{x}}_{ia}+(1-\frac{n_{ia}}{n_{fin}})\mathbf{\theta}\] and \[\Sigma_c=\frac{\Sigma}{n_{fin}} (1 - \frac{n_{ia}}{n_{fin}})\]
7.2.2 Outline of proof (univariate case)
In a case of \(k=1\ldots n_{fin}\) analyses (\(n_{fin}-1\) interim analyses), the estimator
\[\mathbf{\hat{\theta}}=(\mathbf{\bar{x}}_1, \ldots,\mathbf{\bar{x}}_n)\]
is (multivariate) normally distributed with expected value \((\theta,\ldots,\theta)\) and covariance matrix: \[\begin{equation} \label{distIA} \sigma^2 \times \begin{pmatrix} n_1^{-1} & n_2^{-1} & n_3^{-1} & \cdots & n_{fin}^{-1}\\ n_2^{-1} & n_2^{-1} & n_3^{-1} & \cdots & n_{fin}^{-1}\\ n_3^{-1} & n_3^{-1} & n_3^{-1} & \cdots & \vdots\\ \vdots &\vdots & & \ddots & \vdots\\ n_{fin}^{-1} & n_{fin}^{-1} & \ldots & \ldots & n_{fin}^{-1}\\ \end{pmatrix} \end{equation}\]
The well-known theorem on conditional normal distributions, where \(\mathbf{x}\) is partitioned as follows \[\mathbf{x} = \begin{bmatrix} \mathbf{x_1}\\ \mathbf{x_2} \end{bmatrix}, \mathbf{\theta}= \begin{bmatrix} \mathbf{\theta_1}\\ \mathbf{\theta_2} \end{bmatrix}\] \[\mathbf{\Sigma}= \begin{bmatrix} \mathbf{\Sigma_{11}} & \mathbf{\Sigma_{12}}\\ \mathbf{\Sigma_{12}} & \mathbf{\Sigma_{22}}\\ \end{bmatrix}\]
states that the distribution of \(\mathbf{x_2}\) conditional on \(\mathbf{x_1}=\mathbf{\hat{\theta_1}}\) is multivariate normal with expected value \[\mathbf{\theta_2}+\mathbf{\Sigma}_{12}\mathbf{\Sigma}_{11}^{-1} \times (\mathbf{x_1} - \mathbf{\theta_1})\] and covariance matrix \[\mathbf{\Sigma}_{22} - \mathbf{\Sigma}_{12}\mathbf{\Sigma}_{11}^{-1}\mathbf{\Sigma}_{12}\]
From that, the result follows immediately.
7.3 Conjugate Bayesian Estimation of the Full Normal Model
7.3.1 Noninformative Prior
Be \({x_1,x_2, \dots, x_n}\) i.i.d. with \(x_i \sim \mathcal{N}(\mu,\sigma^2)\) with unknown \(\mu\) and \(\sigma^2\). Denoting
\[ s^2=\frac{1}{n-1}\sum_{i=1}^{n}(x-\bar{x})^2\] and \[\bar{x} = \frac{1}{n}\sum_{i=1}^n x_i\] the posterior density (noninformative prior) for \(\sigma^2\) is: \[ \sigma^2|\{x_i\} \sim \text{Scaled-inv}\!-\!\chi^2(n-1,s^2).\] The posterior density for the precision \(\tau\) is: \[ \tau|\{x_i\} \sim \text{Scaled}\!-\!\chi^2(n-1,s^{-2})\] and the posterior for \(\mu\) is: \[\mu|\{x_i\} \sim \text{t}_{n-1}(\bar{x},s^2/n),\] where \(\text{t}_{\nu}(\mu,\sigma^2)\) is the Student-t distribution with \(\nu\) degrees of freedom, location \(\mu\) and scale \(\sigma\).
7.3.2 Conjugate prior for variance
A convenient conjugate prior for the full normal model (\(\mu\) and \(\sigma^2\) unknown), is given by: \[ \mu|\sigma^2 \sim \mathcal{N}(\mu_0,\sigma^2/\kappa_0)\] \[ \sigma^2 \sim \text{Inv}\!-\!\chi^2(\nu_0,\sigma_0^2).\]
7.3.3 Marginal posterior distributions
7.3.3.1 Posterior distribution for \(\sigma^2\)
\[\sigma^2|\{x_i\} \sim \text{Scaled-inv}\!-\!\chi^2(\nu_n,\sigma^2_n),\]
where:
\[\nu_n = \nu_0 + n\]
\[\sigma^2_n = \frac{1}{\nu_n} \left[\nu_0\sigma_0^2+(n-1)s^2+\frac{\kappa_0 n}{\kappa_0+n}(\mu_0-\bar{x})^2 \right].\]
7.3.3.2 Posterior distribution for \(\mu\)
\[\mu|\{x_i\} \sim \text{t}_{\nu_n}(\mu_n,\sigma_n^2/\kappa_n),\] with \[\kappa_n = \kappa_0 +n\] and \[\mu_n = \frac{\kappa_0}{\kappa_0 + n} \mu_0 + \frac{n}{\kappa_0 + n}\bar{x}.\]
7.3.4 Posterior predictive distribution (noninformative prior)
The mean \(\bar{x}_{new}\) of new (independent) future observations \(x_{n+1}, \ldots,x_{n+m}\) having observed \(\{\bar{x}_{old}.s^2_{old}\}\), the mean and variance of \(x_1,\ldots,x_n\) is distributed as: \[\bar{x}_{new}|\bar{x}_{old},s^2_{old} \sim t_{n-1} \left( \bar{x}_{old},s^2_{old}(\frac{1}{n}+\frac{1}{m})\right)\]
7.3.5 The Beta-Binomial distribution
The Beta distribution is a conjugate distribution of the binomial distribution. This fact leads to an analytically tractable compound distribution where one can think of the parameter p in the binomial distribution as being randomly drawn from a beta distribution. Namely, if \[X \sim \text{Bin}(n,\pi)\] \[Prob(X=k|n,\pi)=\binom{n}{k}\pi^k(1-\pi)^{n-k}\] where \(\text{Bin}(n,\pi)\) stands for the binomial distribution, and where \(\pi\) is a random variable with a Beta distribution. \[p(\pi|\alpha,\beta) = \text{Beta}(\pi|\alpha,\beta)=\frac{\pi^{\alpha-1}(1-\pi)^{\beta-1}}{\text{Beta}(\alpha,\beta)}\] then the compound distribution is given by \[f(k|n,\alpha,\beta)=\int_0^1 P(k|s,n)\times p(s|\alpha,\beta)ds\]
\[=\binom{n}{k} \int_0^1 s^{k+\alpha-1} \times (1-s)^{n-k+\beta-1}ds=\binom{n}{k} \frac{\text{B}(k+\alpha,n-k+\beta)}{\text{B}(\alpha,\beta)}\]
References
Dmitrienko, Alexei, and Ming-Dau Wang. 2006. “Bayesian Predictive Approach to Interim Monitoring in Clinicl Trials.” Statisics in Medicine 25: 2178–95.
Gelman, Andrew, John B. Carlin, Hal S. Stern, and Donald B. Rubin. 1995. Bayesian Data Analysis. Chapman & Hall.
Jennison, Christopher, and Bruce W. Turnbull. 2000. Group Sequential Methods with Applications to Clinical Trials. Chapman & Hall/CRC.
Proschan, Michel A., Gordon Lan, K.K., and Janet Turk Wittes. 2006. Statistical Monitoring of Clinical Trials. Springer.
Rufibach, Kaspar, Paul Jordan, and Markus Abt. 2014. “Sequentially Updating the Likelihood of Success of a Phase III Pivotal Time-to-Event Trial Based on Interim Analyses or External Information.” J Biopharm Stat, 1–13.