samplePower: Determine the power of a test with given sample size

In dkahle/bayesRates: Two-Sample Tests and Sample Size Determination from a Bayesian Perspective

Description

Determine the power of a test with given sample size

Usage

samplePower(n, a1, b1, a2, b2, a = a1, b = b1, pi0 = 0.5, pi1 = 1 -
  pi0, c = 1, family = c("binomial", "poisson"))

Arguments

n	the sample size (this is the t parameter in the poisson case)
a1	alpha, the hyperparameter of the prior distribution
b1	beta, the hyperparameter of the prior distribution
a2	alpha, the hyperparameter of the prior distribution
b2	beta, the hyperparameter of the prior distribution
a	alpha, the hyperparameter of the prior distribution under the null
b	beta, the hyperparameter of the prior distribution under the null
pi0	the prior probability of the null hypothesis
pi1	the prior probability of the alternative hypothesis
c	relative loss constant (loss due to type II error divided by loss due to type I error)
family	"binomial" or "poisson", depending on test

See Also

samplePowerBinomial, samplePowerPoisson, samplePowerEst, findSize

Examples

## Not run: 

a1 <- 3; b1 <- 7
a2 <- 7; b2 <- 3
plotBeta(c(a1, a2), c(b1, b2))
samplePower(n = 10, a1, b1, a2, b2)
samplePower(n = 11, a1, b1, a2, b2)
samplePower(n = 1, a1, b1, a2, b2)

samplePower(100, a1, b1, a2, b2)
samplePowerBinomial(100, a1, b1, a2, b2)





# monte carlo check of samplePower
N <- 10000
n <- 100
th1 <- rbeta(N, a1, b1)
y1s <- rbinom(N, n, th1)
th2 <- rbeta(N, a2, b2)
y2s <- rbinom(N, n, th2)
a <- a1; b <- b1 # belief under null
c <- 1 # relative penalty of type II to type I error
pi0 <- .5; pi1 <- .5 # prior likelihoods of hypotheses

samplePower(n = 100, a1, b1, a2, b2)
# exact power should be 91.95%

test <- function(x) bayes_binom_test(x, n, a1, b1, a2, b2)$reject
mean(apply(cbind(y1s, y2s), 1, test))
# simulated power checks out at ~92%
  


# check the probability of type I error of the test
y2s <- rbinom(N, n, th1) # note : not resampled betas
mean(apply(cbind(y1s, y2s), 1, test))
# prob(type I error) ~= 2.5% 

# compare the power of the bayesian-derived test to 
# that of the standard two-tailed two-sample binomial test for the 
# same probability of type I error
classic_test <- function(x) prop.test(x, c(n, n))$p.value < .025

# confirm it's probability of type I error
mean(apply(cbind(y1s, y2s), 1, classic_test)) # it's low...

# check it's power against the same alternative
th2 <- rbeta(N, a2, b2)
y2s <- rbinom(N, n, th2)
mean(apply(cbind(y1s, y2s), 1, classic_test))
# power ~= 87.5%, lower than the bayesian test of 91.52%








a1 <- 1; b1 <- 1 # E = 1
a2 <- 3; b2 <- 1 # E = 3
plot_gamma(c(a1, a2), c(b1, b2))
samplePower(n = 2, a1, b1, a2, b2, family = "poisson")
samplePower(n = 10, a1, b1, a2, b2, family = "poisson")
samplePower(n = 11, a1, b1, a2, b2, family = "poisson")
samplePower(n = 10:11, a1, b1, a2, b2, family = "poisson")

samplePower(n = 250, a1, b1, a2, b2, family = "poisson") 
samplePowerEst(n = c(250, 500, 1000), a1, b1, a2, b2) # 5s, off by .001



a1 <- 4; b1 <- 2
a2 <- 6; b2 <- 2
plot_gamma(c(a1, a2), c(b1, b2))
samplePower(n = 20, a1, b1, a2, b2, family = "poisson")



## End(Not run)

dkahle/bayesRates documentation built on May 15, 2019, 9:07 a.m.