data.type <- "Bernoulli"
n.t <- 100
n.c <- 100
# Simulate three historical datasets
historical <- matrix(0, ncol=2, nrow=3)
historical[1,] <- c(70, 100)
historical[2,] <- c(60, 100)
historical[3,] <- c(50, 100)
# Generate sampling priors
set.seed(1)
b_st1 <- b_st2 <- 1
b_sc1 <- b_sc2 <- 1
samp.prior.mu.t <- rbeta(50000, b_st1, b_st2)
samp.prior.mu.c <- rbeta(50000, b_st1, b_st2)
# The null hypothesis here is H0: mu_t - mu_c >= 0. To calculate power,
# we can provide samples of mu.t and mu.c such that the mass of mu_t - mu_c < 0.
# To calculate type I error, we can provide samples of mu.t and mu.c such that
# the mass of mu_t - mu_c >= 0.
sub_ind <- which(samp.prior.mu.t < samp.prior.mu.c)
# Here, mass is put on the alternative region, so power is calculated.
samp.prior.mu.t <- samp.prior.mu.t[sub_ind]
samp.prior.mu.c <- samp.prior.mu.c[sub_ind]
N <- 10 # N should be larger in practice
result <- power.two.grp.random.a0(data.type=data.type, n.t=n.t, n.c=n.c, historical=historical,
samp.prior.mu.t=samp.prior.mu.t, samp.prior.mu.c=samp.prior.mu.c,
delta=0, nMC=10000, nBI=250, N=N)
summary(result)
Run the code above in your browser using DataLab