sampleAlpha: Determine the probability of type I error of a test with...

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

Description

Determine the probability of type I error of a test with given sample size

Usage

sampleAlpha(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

sampleAlphaBinomial, sampleAlphaPoisson, findSize

Examples

## Not run: 


 

# generate samples of size 30
N <- 10000#00
n <- 30

a1 <- 3; b1 <- 7
a2 <- 7; b2 <- 3
plotBeta(c(a1, a2), c(b1, b2))

pi <- .3
pi <- rbeta(N, a1, b1)
data <- data.frame(
  x1 = rbinom(N, n, pi),
  x2 = rbinom(N, n, pi)  
)
test <- function(x) bayesBinomTest(x, n, a1, b1, a2, b2)$reject
mean( apply(data, 1, test) ) # .061832 at 1E6
sampleAlpha_ryan("binomial", n, 1, a1, b1) # .0204
sampleAlpha(n, a1, b1, a2, b2) 
bayesRates:::sampleAlphaPoissonCpp(n, a1, b1, a2, b2, a1, b1, .5, .5, 1)




# note that the significance depends on all parameters
a1 <- 3; b1 <- 7
a2 <- 90; b2 <- 10
plotBeta(c(a1, a2), c(b1, b2))
f <- function(x) bayesBinomTest(x, n, a1, b1, a2, b2)$reject
f(c(11, 8))
mean( apply(data, 1, f) )
sampleAlpha(n, a1, b1, a2, b2) 











if(FALSE){


sim_alpha <- function(pi, N = 5E4, n = 30){

  data <- data.frame(
    x1 = rbinom(N, n, pi),
    x2 = rbinom(N, n, pi)  
  )


  a1 <- 3; b1 <- 7
  a2 <- 7; b2 <- 3

  f <- function(x) bayes_binom_test(x, n, a1, b1, a2, b2)$reject

  mean( apply(data, 1, f) )
}

sim_alpha(.3)

s <- seq(.01, .99, length.out = 250)
t1 <- sapply(as.list(s), function(x){
  message(".", appendLF = F)
  sim_alpha(x)
})

qplot(s, t1, geom = "line") +
  labs(x = expression(paste(pi, " = ", pi[1], " = ", pi[2])), y = expression(alpha)) +
  coord_equal()
  
  
sim_alpha_classic <- function(pi, N = 5E3, n = 30){

  data <- data.frame(
    x1 = rbinom(N, n, pi),
    x2 = rbinom(N, n, pi)  
  )
  
  f <- function(x) prop.test(x, c(n, n))$p.value < .05

  mean( apply(data, 1, f) )
}
sim_alpha_classic(.3)  
  
s <- seq(.01, .99, length.out = 25)
t1 <- sapply(as.list(s), function(x){
  message(".", appendLF = F)
  sim_alpha_classic(x)
})

qplot(s, t1, geom = "line") +
  labs(x = expression(paste(pi, " = ", pi[1], " = ", pi[2])), y = expression(alpha))





}




## End(Not run)

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