# R/bpower.s In harrelfe/Hmisc: Harrell Miscellaneous

#### Documented in ballocationbpowerbpower.simbsamsize

```bpower <- function(p1, p2, odds.ratio, percent.reduction, n, n1, n2,
alpha=.05)
{
if(!missing(odds.ratio))
p2 <- p1*odds.ratio/(1-p1+p1*odds.ratio)
else if(!missing(percent.reduction))
p2 <- p1*(1-percent.reduction/100)

if(!missing(n)) {
n1 <- n2 <- n/2
}
z <- qnorm(1-alpha/2)
q1 <- 1-p1
q2 <- 1-p2
pm <- (n1*p1+n2*p2)/(n1+n2)
ds <- z*sqrt((1/n1 + 1/n2)*pm*(1-pm))
ex <- abs(p1-p2)
sd <- sqrt(p1*q1/n1+p2*q2/n2)
c(Power = 1-pnorm((ds-ex)/sd)+pnorm((-ds-ex)/sd) )
}

bsamsize <- function(p1, p2, fraction=.5, alpha=.05, power=.8)
{
z.alpha <- qnorm(1-alpha/2)
z.beta  <- qnorm(power)

ratio <- (1-fraction)/fraction
p <- fraction*p1+(1-fraction)*p2

n1 <- (z.alpha * sqrt((ratio+1) * p * (1-p)) +
z.beta * sqrt(ratio * p1 * (1-p1) + p2 * (1 - p2))
)^2/ratio/((p1-p2)^2)

n2 <- ratio*n1
c(n1=n1, n2=n2)
}

ballocation <- function(p1, p2, n, alpha=.05)
{
q1 <- 1-p1
q2 <- 1-p2

f.minvar.diff <- 1/(1+sqrt(p2*q2/(p1*q1)))
f.minvar.ratio <- 1/(1+sqrt(p1*q2/p2/q1))

z <- c(fraction.group1.min.var.diff=f.minvar.diff,
fraction.group1.min.var.ratio=f.minvar.ratio,
fraction.group1.min.var.logodds=1-f.minvar.diff)

if(!missing(n)) {
possf <- seq(.001,.999,length=1000)
pow <- bpower(p1, p2, n1=n*possf, n2=n*(1-possf), alpha=alpha)
## fun <- function(f, n, p1, p2, alpha) bpower(p1, p2, n1=f*n, n2=(1-f)*n, alpha=alpha)
## f.maxpow <- optimize(fun, lower=.01, upper=.99, maximum=T,
##                      n=n, p1=p1, p2=p2, alpha=alpha)\$maximum
f <- possf[pow==max(pow)]
f <- f[abs(f-.5)==min(abs(f-.5))]
z <- c(z, fraction.group1.max.power=f[1])
}
z
}

bpower.sim <- function(p1, p2, odds.ratio, percent.reduction, n, n1, n2,
alpha=.05, nsim=10000)
{
if(!missing(odds.ratio))
p2 <- p1*odds.ratio/(1-p1+p1*odds.ratio)
else if(!missing(percent.reduction))
p2 <- p1*(1-percent.reduction/100)

if(!missing(n)) {
n1 <- n2 <- round(n/2)
}
n <- n1+n2

if(length(p1)+length(p2)+length(n1)+length(n2)+length(alpha)+length(nsim)!=6)
stop('all arguments must have length 1')

chi2 <- qchisq(1-alpha, 1)

d1 <- rbinom(nsim, n1, p1)
d2 <- rbinom(nsim, n2, p2)
chisq <- n*(d1*(n2-d2)-(n1-d1)*d2)^2/(d1+d2)/(n-d1-d2)/n1/n2
power <- mean(chisq>chi2)
se <- sqrt(power*(1-power)/nsim)
c(Power=power,Lower=power-1.96*se,Upper=power+1.96*se)
}
```
harrelfe/Hmisc documentation built on May 19, 2024, 4:13 a.m.