Description Usage Arguments See Also Examples
Determine the power of a test with given sample size
1 2 | samplePower(n, a1, b1, a2, b2, a = a1, b = b1, pi0 = 0.5, pi1 = 1 -
pi0, c = 1, family = c("binomial", "poisson"))
|
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 |
samplePowerBinomial
, samplePowerPoisson
, samplePowerEst
, findSize
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 | ## 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.