Notes on Dunnett’s Method for Conversion Data
Dunnett’s method is a way to compare any number of treatments against a common control, all while preserving the overall type 1 error rate. This is a far more efficient approach than using Bonferroni corrections or doing a series of A/B tests. Here are my notes on the procedure.
Author’s Note: This post is me pulling together notes and reworking some things, and so things are even less guaranteed to be correct than in other posts (which have no guarantee of correctness, so make of that what you will). If you find errors please let me know.
Dunnett’s method (or Dunnett’s test or Dunnett’s procedure) is a method for analyzing an A/B/n experiment where we can test a number of challengers, each versus a single control, while controling the probability of at least one false alarm in the procedure. This is often called multiple comparisons with control (MCC). In the case where we’re working with conversions we can avoid integration across possible variances and still get good results, which simplifies the expressions involved.
Let there be \(m\) challenger treatments indexed by \(i=1,\ldots,m\), and let there be a single control treatment indexed by \(i=0\). The null hypothesis is always \(H_0: p_0 = p_1 = \ldots = p_m\) but the alternative can be one-sided or two-sided:
- \(H_A: p_i \ne p_0 ~~\text{for some}~~i \in \{1,\ldots,m\}\)
- \(H_A: p_i > p_0 ~~\text{for some}~~i \in \{1,\ldots,m\}\)
- \(H_A: p_i < p_0 ~~\text{for some}~~i \in \{1,\ldots,m\}\)
Dunnett’s test also produces corresponding confidence intervals around the comparison of each non-control treatment with the control, so the intervals will fall around
\[ \mu_i - \mu_0~. \]
These intervals have simultaneous coverage of \(100(1-\alpha)\%\).
Let’s establish some notation for working out some of the details here. Let \(p_i\) represent the proportion converted from treatment \(i\), and \(n_i\) the number of samples in treatment \(i\). Also, let \(\sigma_i^2 = p_i(1 - p_i)/n_i\). For exchangeability reasons, we should have \(n_1 = \ldots = n_m\), but it need not be that \(n_0 = n_1\). Indeed, since the control data is used in comparisons to each of the \(m\) challengers we will find that having a larger sample size for \(n_0\) than for the other treatments is appropriate.
Dunnett’s Test: Alternative \(p_i > p_0\)
I’ll use the following identity in changing some probabilities to conditional expectations:
\[ P(X > Y) = \int_{\Omega_Y} P(X>y)f(y) ~\mathrm{d}y= \mathbb{E}\left[P(X > Y | Y)\right]~. \]
We want to maintain the overall type 1 error rate \(\alpha\) (that is, the probability of at least one false alarm among the tests \(p_i > p_0\) for \(i=1,\ldots,m\)). We first need to find a set of constants \(d_i\) so that the probability (under \(H_0\)) of no type 1 errors occuring is \(1-\alpha\). This is the same as saying our lower confidence bound \(\hat{p}_i - \hat{p}_0 - d \sqrt{\sigma_i^2 + \sigma_0^2}\) should be less than \(p_i - p_0\) for \(i=1,\ldots,m\) with probability \(1-\alpha\): \[\begin{align*} 1-\alpha &= P\left(p_i - p_0 > \hat{p}_i - \hat{p}_0 - d_i \sqrt{\sigma_i^2 + \sigma_0^2} \quad \text{for}~i=1,\ldots,m\right) \\ &= \mathbb{E}\left[ P\left(d_i \sqrt{\sigma_i^2 + \sigma_0^2} > (\hat{p}_i - p_i)-(\hat{p}_0-p_0)~\forall~i > 0~~|~~\hat{p}_0\right) \right] \\ &= \mathbb{E}\left[ P\left(d_i \sqrt{1 + \frac{\sigma_0^2}{\sigma_i^2}} > \underbrace{\frac{\hat{p}_i - p_i}{\sigma_i}}_{Z_i} - \frac{\sigma_0}{\sigma_i} \underbrace{\frac{\hat{p}_0 - p_0}{\sigma_0}}_{Z_0} ~\forall~i > 0~~|~~\hat{p}_0\right) \right] \end{align*}\] where \(Z_0 \overset{\cdot}{\sim} \mathrm{N}(p_0, \sigma_0^2)\) and \(Z_i \overset{\cdot}{\sim} \mathrm{N}(p_i, \sigma_i^2)\).
We can rewrite this as \[\begin{align*} 1-\alpha &= \mathbb{E}\left[ P\left(d_i \sqrt{1 + \frac{\sigma_0^2}{\sigma_i^2}} > Z_i - \frac{\sigma_0}{\sigma_i} Z_0~\forall~i > 0~~|~~Z_0\right) \right] \\ &= \mathbb{E}\left[P\left(Z_i < \frac{\sigma_0}{\sigma_i} Z_0 + d_i \sqrt{1 + \frac{\sigma_0^2}{\sigma_i^2}}~\forall~i > 0 ~~|~~ Z_0\right)\right] \\ &= \int_{-\infty}^\infty \prod_{i=1}^m \Phi\left(\sqrt{\frac{\sigma_0^2}{\sigma_i^2}} z+ d_i \sqrt{1+\frac{\sigma_0^2}{\sigma_i^2}}\right)~\mathrm{d}\Phi(z)~. \end{align*}\] Since we can exchange positive indecies without affecting our analysis we should have \(n_1 = \ldots = n_m\) (no particular index should get more samples than the others), and for that reason we can choose the \(d_i\)’s to be a single constant \(d\), then our expression simplifies to
\[ 1-\alpha=\int_{-\infty}^\infty \left[\Phi(\sqrt{\lambda}z+d \sqrt{1 + \lambda})\right]^m~\mathrm{d}\Phi(z)~. \]
where \(\lambda = \frac{\sigma_0^2}{\sigma_i^2}\). Note that the form of this integral does not rely on \(n_i\) or \(n_0\) anymore. We’ll see later optimal allocation occurs when \(n_0 = \sqrt{m} n_i\).
We can construct \(100(1-\alpha)\%\) confidence intervals in the following way
\[ \left[ \hat{p}_i - \hat{p}_0 - d \sqrt{\frac{\hat{p}_0 (1 - \hat{p}_0)}{n_0} + \frac{\hat{p} _i (1 - \hat{p}_i)}{n_i}}, \infty\right) \quad\text{where}\quad \hat{p}_i = \frac{X_i}{n_i}~. \]
Dunnett’s Test: Alternative \(p_i < p_0\)
Just as above, we can work out \[\begin{align*} 1 - \alpha &= P\left(p_i - p_0 < \hat{p}_i - \hat{p}_0 + d_i \sqrt{\sigma_i^2 + \sigma_0^2}~\text{for}~i=1,\ldots,m\right) \\ &= \int_{-\infty}^\infty\left[1 - \Phi(\sqrt{\lambda}z - d_i \sqrt{1+\lambda})\right]^m~\mathrm{d}\Phi(z)~. \end{align*}\] Hsu (1996) gives a proof that this constant \(d\) is identical to that in the case of the the \(p_i > p_0\) alternative.
Dunnett’s Test: Alternative \(p_i \ne p_0\)
Here, we need to calculate a new two-sided \(d\) using the following:
\[ 1-\alpha = \int_{-\infty}^\infty \left[ \Phi\left(\sqrt\lambda z + d \sqrt{1 + \lambda}\right) - \Phi\left(\sqrt\lambda z - d \sqrt{1 + \lambda}\right) \right]^m~\mathrm{d}\Phi(z)~. \]
We can construct simultaneous \(100(1-\alpha)\%\) confidence intervals for \(i=1,\ldots,m\) as
\[ \hat{p}_i - \hat{p}_0 \pm d \sqrt{\frac{\hat{p}_i(1 - \hat{p}_i)}{n_i} + \frac{\hat{p}_0(1 - \hat{p}_0)}{n_0}}~. \]
Powering the Thing
There’s lots of different alternatives we could use for power analysis. But the situation needing the largest sample sizes is when:
- Every challenger is some fixed improvement from the control treatment: \(p_1 = \ldots = p_m > p_0\), and
- We require that we detect all such effects with probability \(1-\beta\).
Starting with the \(p_i > p_0\) alternative, a little algebra gets us to \[\begin{align*} 1 - \beta &= P\left(\hat{p}_i - \hat{p}_0 - d \sqrt{\sigma_i^2 + \sigma_0^2} > 0~~\text{for}~~i=1,\ldots,m\right) \\ &= \mathbb{E}\left[ \prod_{i=1}^m P\left( \frac{\hat{p}_i - p_i}{\sigma_i} - \frac{\sigma_0}{\sigma_i} \frac{\hat{p}_0 - p_0}{\sigma_0} > \frac{d \sqrt{\sigma_i^2 + \sigma_0^2} - (p_i - p_0)}{\sigma_i} ~~\Bigg|~~\hat{p}_0 \right) \right] \\ &\approx \mathbb{E}\left[ \prod_{i=1}^m P\left( Z_i - \frac{\sigma_0}{\sigma_i} Z_0 > d \sqrt{1 + \frac{\sigma_0^2}{\sigma_i^2}} - \frac{p_i - p_0}{\sigma_i} ~~\Bigg|~~Z_0 \right) \right] \\ &= \int_{-\infty}^\infty \left[ 1 - \Phi\left( \sqrt\lambda z + d \sqrt{1 + \lambda} - \frac{p_i - p_0}{\sqrt{p_i (1 - p_i)}}\sqrt{n_i} \right) \right]^m ~\mathrm{d}\Phi(z)~. \end{align*}\]
When \(p_i < p_0\) then the power calculation becomes \[\begin{align*} 1 - \beta &= P\left(\hat{p}_i - \hat{p}_0 + d \sqrt{\sigma_i^2 + \sigma_0^2} < 0~~\text{for}~~i=1,\ldots,m\right) \\ &= \mathbb{E}\left[ \prod_{i=1}^m P\left( \frac{\hat{p}_i - p_i}{\sigma_i} - \frac{\sigma_0}{\sigma_i} \frac{\hat{p}_0 - p_0}{\sigma_0} < -\frac{p_i - p_0}{\sigma_i} - \frac{d \sqrt{\sigma_i^2 + \sigma_0^2}}{\sigma_i} ~~\Bigg|~~\hat{p}_0 \right) \right] \\ &\approx \mathbb{E}\left[ \prod_{i=1}^m P\left( Z_i - \frac{\sigma_0}{\sigma_i} Z_0 < -\frac{p_i - p_0}{\sigma_i} -d \sqrt{1 + \frac{\sigma_0^2}{\sigma_i^2}} ~~\Bigg|~~Z_0 \right) \right] \\ &= \int_{-\infty}^\infty \left[ \Phi\left( \sqrt\lambda z - d \sqrt{1 + \lambda} - \frac{p_i - p_0}{\sqrt{p_i (1 - p_i)}}\sqrt{n_i} \right) \right]^m ~\mathrm{d}\Phi(z)~. \end{align*}\]
The two-sided alternative can be phrased in terms of the detection event
\[ \left\{\left|\hat{p}_i - \hat{p}_0\right| > d \sqrt{\sigma_i^2 + \sigma_0^2}\right\} \]
which we can write as
\[ Z_i > \sqrt{\lambda}Z_0 + d \sqrt{1 + \lambda} - \frac{p_i-p_0}{\sigma_i} \quad\text{or}\quad Z_i < \sqrt{\lambda}Z_0 - d \sqrt{1 + \lambda} - \frac{p_i-p_0}{\sigma_i} \]
and by similar actions to the above arrive at \[\begin{align*} 1-\beta = \int_{-\infty}^\infty \Biggl[& 1 - \Phi\left( \sqrt\lambda z + d \sqrt{1 + \lambda} - \frac{p_i - p_0}{\sqrt{p_i (1 - p_i)}}\sqrt{n_i} \right) \\&+ \Phi\left( \sqrt\lambda z - d \sqrt{1 + \lambda} - \frac{p_i - p_0}{\sqrt{p_i (1 - p_i)}}\sqrt{n_i} \right) \Biggr]^m ~\mathrm{d}\Phi(z)~. \end{align*}\]
Optimal Tradeoff Between \(n_0\) and \(n_i\)
The quantity we’re estimating in all cases is \(\hat{p}_i - \hat{p}_0\) so one possible way of finding the optimal tradeoff in sample size is by minimizing the sum of such variances (the \(A\)-optimal allocation, which minimizes the trace of the covariance matrix \(\boldsymbol\Sigma\)). Let’s operate under the null hypothesis so \(p = p_0 = p_1 = \ldots = p_m\), and let \(n_0 = \gamma n_1\) assuming \(n_1 = \ldots = n_m = n\). We have the constraint that the overall sample size is fixed at \(N\) where \(N = \gamma n + m n\), so
\[ n = \frac{\gamma + m}{N}~. \]
The sum of the variances under \(H_0\) is then \[\begin{align*} \mathrm{tr}(\boldsymbol\Sigma) &= m\left(\frac{p(1-p)}{n} + \frac{p(1-p)}{\gamma n}\right) \\ &= m\left(1 + \gamma + m + \frac{m}{\gamma}\right) \frac{p(1-p)}{N} \end{align*}\] Differentiating with respect to \(\gamma\) and solving for zero we find
\[ \frac{\partial \mathrm{tr}(\boldsymbol\Sigma)}{\partial \gamma} = \left(m - \frac{m^2}{\gamma^2}\right)\frac{p(1-p)}{N} = 0 \]
so
\[ m - \frac{m^2}{\gamma^2} = 0 \quad\text{and}\quad\gamma=\sqrt{m}~. \]
Thus the \(A\)-optimal solution sets \(n_0 = \sqrt{m} n_i\).
One might ask if the \(D\)-optimal solution, where we minimize the determinant of the covariance matrix \(\boldsymbol\Sigma\), is also a good choice for allocation. Under \(H_0\) \[\begin{align*} \boldsymbol\Sigma &= \mathrm{Var}(\hat{p}_0) \mathbf{J} + \mathrm{Var}(\hat{p}_i) \mathbf{I} = \left(1+\frac{m}{\gamma}\right)\frac{p(1-p)}{N} \mathbf{J} + (\gamma + n) \frac{p(1-p)}{N} \mathbf{I} \\ \boldsymbol\Sigma &\propto \left(1+\frac{m}{\gamma}\right) \mathbf{J} + (\gamma + m) \mathbf{I} \end{align*}\] Using the matrix determinant identity
\[ \mathrm{det}\left(\mathbf{A} + \mathbf{u}\mathbf{v}^\mathrm{T}\right) = \left(1 + \mathbf{v}^\mathrm{T} \mathbf{A}^{-1} \mathbf{u}\right) ~\mathrm{det}(\mathbf{A}) \]
and letting
\[ \mathbf{u} = \mathbf{v} = \sqrt{1+\frac{m}{\gamma}}\mathbf{1} \quad\text{and}\quad \mathbf{A} = (\gamma + m) \mathbf{I} \]
where
\[ \mathbf{A}^{-1} = \frac{1}{\gamma + m}\mathbf{I} \quad\text{and}\quad \mathrm{det}(\mathbf{A}) = (\gamma + m)^{m} \]
we find
\[ \left(1 + \mathbf{v}^\mathrm{T} \mathbf{A}^{-1} \mathbf{u}\right) ~\mathrm{det}(\mathbf{A}) = \left(1 + \left(\frac{\gamma + m}{\gamma}\right)\left(\frac{1}{\gamma + m}\right)\right) (\gamma + m)^m = \frac{1 + \gamma}{\gamma} (\gamma + m)^m~. \]
Taking the derivative and setting it to zero we find the \(D\)-criterion is minimized at
\[ \gamma = \frac{1 - m + \sqrt{5m^2 - 2m + 1}}{2m}~. \]
When \(m=1\) then \(\gamma = 1\) but as \(m\) increases \(\gamma\) decreases quickly to its limit \((\sqrt{5} - 1)/2 \approx 0.618\).
The emphasis of the \(D\)-criterion on reducing correlations causes it to drive the variance of our individual estimates up in a way that I don’t think most people would enjoy.
Appendix A - Test Code
Below is R code I used to implement and test my formulation of Dunnett’s test.
##
## Verify that the Dunnett's constant is matching our table.
##
alpha <- 0.05
book.table <- data.frame(
'k'=2:17,
'd'= c(1.645, 1.916, 2.062, 2.160, 2.234, 2.292, 2.341, 2.381, 2.417, 2.448,
2.476, 2.502, 2.525, 2.546, 2.565, 2.584)
)
constant.dunnetts.onesided <- function(alpha, ns){
g <- function(d){
f <- Vectorize(function(z){
exp(
sum(
pnorm(sqrt(ns[-1]/ns[1]) * z + d * sqrt(1 + ns[-1]/ns[1]), log.p=TRUE)
) + dnorm(z, log=TRUE)
)
})
integrate(f, lower=-Inf, upper=Inf, abs.tol=1e-8)$value
}
rt <- uniroot(
function(q){
g(q) - (1 - alpha)
},
interval=c(0, 100),
tol=1e-8
)
return(rt$root)
}
for(i in 1:dim(book.table)[1]){
ns <- rep(1, book.table[i, 'k'])
d <- constant.dunnetts.onesided(alpha, ns)
if(abs(book.table[i,'d'] - d) > 0.001) stop('Error')
}
print("One-Sided constants matched the book values for infinite degrees of freedom.")
# ------------------------------------------------------------------------------
##
## Two-Sided Dunnett Tests
##
alpha <- 0.05
book.table <- data.frame(
'k'=2:17,
'd'= c(1.960, 2.212, 2.349, 2.442, 2.512, 2.567, 2.613, 2.652, 2.686, 2.716,
2.743, 2.767, 2.790, 2.810, 2.829, 2.846)
)
constant.dunnetts.twosided <- function(alpha, ns){
g <- function(d){
f <- Vectorize(function(z){
exp(
sum(
log(
pnorm(sqrt(ns[-1]/ns[1]) * z + d * sqrt(1 + ns[-1]/ns[1]))
-
pnorm(sqrt(ns[-1]/ns[1]) * z - d * sqrt(1 + ns[-1]/ns[1]))
)
) + dnorm(z, log=TRUE)
)
})
integrate(f, lower=-Inf, upper=Inf, abs.tol=1e-8)$value
}
rt <- uniroot(
function(q){
g(q) - (1 - alpha)
},
interval=c(0, 100),
tol=1e-8
)
return(rt$root)
}
for(i in 1:dim(book.table)[1]){
ns <- rep(1, book.table[i, 'k'])
d <- constant.dunnetts.twosided(alpha, ns)
if(abs(book.table[i,'d'] - d) > 0.001) stop('Error')
}
print("Two-Sided constants matched the book values for infinite degrees of freedom.")
# ------------------------------------------------------------------------------
##
## Define the functions for one-sided and two-sided confidence intervals and
## planning. Note that .d is letting me pass through the constant instead of
## calculating it each run of the simulation study I do later.
##
label.dunnett <- function(v){
names(v) <- sprintf('[%d] - [0]', 1:(length(ns) - 1))
return(v)
}
ci.dunnetts.twosided <- function(xs, ns, alpha, .d=NULL){
if(!is.null(.d)){
d <- .d
} else {
d <- constant.dunnetts.twosided(alpha, ns)
}
ps <- xs / ns
ss <- sqrt(ps[1] * (1 - ps[1]) / ns[1] + ps[-1] * (1 - ps[-1]) / ns[-1])
lcl <- label.dunnett(ps[-1] - ps[1] - d * ss)
ucl <- label.dunnett(ps[-1] - ps[1] + d * ss)
pte <- label.dunnett(ps[-1] - ps[1])
return(
list(
'Lower'=lcl,
'Point Estimate'=pte,
'Upper'=ucl
)
)
}
ci.dunnetts.onesided <- function(xs, ns, alpha, alternative='>', .d=NULL){
if(!is.null(.d)){
d <- .d
} else {
d <- constant.dunnetts.onesided(alpha, ns)
}
ps <- xs / ns
ss <- sqrt(ps[1] * (1 - ps[1]) / ns[1] + ps[-1] * (1 - ps[-1]) / ns[-1])
if(alternative == '>'){
lcl <- label.dunnett(ps[-1] - ps[1] - d * ss)
ucl <- label.dunnett(rep(Inf, length(ns) - 1))
} else if(alternative == '<'){
lcl <- label.dunnett(rep(-Inf, length(ns) - 1))
ucl <- label.dunnett(ps[-1] - ps[1] + d * ss)
} else {
stop('The "alternative" argument must be ">" or "<" only.')
}
pte <- label.dunnett(ps[-1] - ps[1])
return(
list(
'Lower'=lcl,
'Point Estimate'=pte,
'Upper'=ucl
)
)
}
power.dunnetts.kind <- function(m, kind){
ns <- rep(1, m+1)
if(kind == 'A-optimal'){
ns[1] <- sqrt(m)
} else if(kind == 'D-optimal') {
ns[1] <- (1 - m + sqrt(5*m^2 - 2*m + 1))/(2 * m)
} else if(kind == 'balanced') {
# do nothing
} else {
stop('"kind" must be "A-optimal", "D-optimal", or "balanced" only.')
}
ns
}
power.dunnetts.onesided.gt <- function(alpha, beta, m, p0, p1, ns, .d=NULL){
if(!is.null(.d)){
d <- .d
} else {
d <- constant.dunnetts.onesided(alpha, ns)
}
lambda <- (p0 * (1 - p0) / ns[1]) / (p1 * (1 - p1) / ns[-1])
g <- Vectorize(function(n){
f <- Vectorize(function(z){
exp(
sum(
log(
1 - pnorm(sqrt(lambda) * z + d * sqrt(1 + lambda) - sqrt(n)*(p1 - p0)/sqrt(p1*(1-p1)))
)
) + dnorm(z, log=TRUE)
)
})
integrate(f, lower=-Inf, upper=Inf, abs.tol=1e-8)$value
})
ix <- 2^(0:32)
iy <- g(ix) - (1 - beta)
interval <- ix[c(max(which(iy < 0)), min(which(iy > 0)))]
rt <- uniroot(
function(q){
g(q) - (1 - beta)
},
interval=interval,
tol=1e-8
)
return(list('ns'=rt$root * ns, 'd'=d))
}
power.dunnetts.onesided.lt <- function(alpha, beta, m, p0, p1, ns, .d=NULL){
if(!is.null(.d)){
d <- .d
} else {
d <- constant.dunnetts.onesided(alpha, ns)
}
lambda <- (p0 * (1 - p0) / ns[1]) / (p1 * (1 - p1) / ns[-1])
g <- Vectorize(function(n){
f <- Vectorize(function(z){
exp(
sum(
pnorm(sqrt(lambda) * z - d * sqrt(1 + lambda) - sqrt(n)*(p1 - p0)/sqrt(p1*(1-p1)), log.p=TRUE)
) + dnorm(z, log=TRUE)
)
})
integrate(f, lower=-Inf, upper=Inf, abs.tol=1e-8)$value
})
ix <- 2^(0:32)
iy <- g(ix) - (1 - beta)
interval <- ix[c(max(which(iy < 0)), min(which(iy > 0)))]
rt <- uniroot(
function(q){
g(q) - (1 - beta)
},
interval=interval,
tol=1e-8
)
return(list('ns'=rt$root * ns, 'd'=d))
}
power.dunnetts.twosided <- function(alpha, beta, m, p0, p1, ns, .d=NULL){
if(!is.null(.d)){
d <- .d
} else {
d <- constant.dunnetts.twosided(alpha, ns)
}
lambda <- (p0 * (1 - p0) / ns[1]) / (p1 * (1 - p1) / ns[-1])
g <- Vectorize(function(n){
f <- Vectorize(function(z){
exp(
sum(
log(
1 - pnorm(sqrt(lambda) * z + d * sqrt(1 + lambda) - sqrt(n)*(p1 - p0)/sqrt(p1*(1-p1)))
+ pnorm(sqrt(lambda) * z - d * sqrt(1 + lambda) - sqrt(n)*(p1 - p0)/sqrt(p1*(1-p1)))
)
) + dnorm(z, log=TRUE)
)
})
integrate(f, lower=-Inf, upper=Inf, abs.tol=1e-8)$value
})
ix <- 2^(0:32)
iy <- g(ix) - (1 - beta)
interval <- ix[c(max(which(iy < 0)), min(which(iy > 0)))]
rt <- uniroot(
function(q){
g(q) - (1 - beta)
},
interval=interval,
tol=1e-8
)
return(list('ns'=rt$root * ns, 'd'=d))
}
power.dunnetts <- function(alpha, beta, m, p0, p1, alternative='>', kind='A-optimal', .d=NULL){
ns <- power.dunnetts.kind(m, kind)
if(alternative == '>'){
rslt <- power.dunnetts.onesided.gt(alpha, beta, m, p0, p1, ns, .d)
} else if(alternative == '<'){
rslt <- power.dunnetts.onesided.lt(alpha, beta, m, p0, p1, ns, .d)
} else if(alternative == '!='){
rslt <- power.dunnetts.twosided(alpha, beta, m, p0, p1, ns, .d)
} else {
stop('"alternative" must be ">", "!=", or "<" only.')
}
return(
list(
alternative=alternative,
kind=kind,
ns=rslt$ns,
d=rslt$d
)
)
}
# ------------------------------------------------------------------------------
##
## Simulation of Two-Sided CIs. Looking for 95% Coverage and 90% Power
##
library(parallel)
cl <- parallel::makeForkCluster(parallel::detectCores() - 2)
parallel::clusterSetRNGStream(cl)
alpha <- 0.05
beta <- 0.10
p0 <- 0.11
p1 <- 0.115
m <- 4
plan <- power.dunnetts(alpha, beta, m, p0, p1, alternative='!=')
ns <- ceiling(plan$ns)
d <- plan$d
f <- function(){
xs <- rbinom(length(ns), ns, c(p0, rep(p1, m)))
rslt <- ci.dunnetts.twosided(xs, ns, alpha, .d=d)
coverage <- all(rslt$Lower < p1 - p0 & rslt$Upper > p1 - p0)
power <- all(rslt$Lower > 0 | rslt$Upper < 0)
return(c(coverage, power))
}
parallel::clusterExport(
cl,
c('alpha', 'p0', 'p1', 'ns', 'm', 'd', 'f',
"ci.dunnetts.onesided", "ci.dunnetts.twosided",
"constant.dunnetts.onesided", "constant.dunnetts.twosided",
"label.dunnett", "power.dunnetts", "power.dunnetts.kind",
"power.dunnetts.onesided.gt", "power.dunnetts.onesided.lt",
"power.dunnetts.twosided")
)
rslt <- do.call(cbind, parallel::clusterApplyLB(cl, 1:512, function(.){
apply(
vapply(1:10000, function(.){f()}, c(TRUE,TRUE)),
1,
mean
)
}))
rslt <- apply(rslt, 1, mean)
print(sprintf('TwoSided Coverage: %0.3f%%, Power: %0.3f%%', 100*rslt[1], 100*rslt[2]))
plan <- power.dunnetts(alpha, beta, m, p0, p1, alternative='>')
ns <- ceiling(plan$ns)
d <- plan$d
f <- function(){
xs <- rbinom(length(ns), ns, c(p0, rep(p1, m)))
rslt <- ci.dunnetts.onesided(xs, ns, alpha, alternative='>', .d=d)
coverage <- all(rslt$Lower < p1 - p0 & rslt$Upper > p1 - p0)
power <- all(rslt$Lower > 0 | rslt$Upper < 0)
return(c(coverage, power))
}
parallel::clusterExport(
cl,
c('alpha', 'p0', 'p1', 'ns', 'm', 'd', 'f',
"ci.dunnetts.onesided", "ci.dunnetts.twosided",
"constant.dunnetts.onesided", "constant.dunnetts.twosided",
"label.dunnett", "power.dunnetts", "power.dunnetts.kind",
"power.dunnetts.onesided.gt", "power.dunnetts.onesided.lt",
"power.dunnetts.twosided")
)
rslt <- do.call(cbind, parallel::clusterApplyLB(cl, 1:512, function(.){
apply(
vapply(1:10000, function(.){f()}, c(TRUE,TRUE)),
1,
mean
)
}))
rslt <- apply(rslt, 1, mean)
print(sprintf('OneSided Coverage: %0.3f%%, Power: %0.3f%%', 100*rslt[1], 100*rslt[2]))
parallel::stopCluster(cl)