First the theory. Question 2 asks: Why doesn't the true parameter $\theta_a$ under the alternative hypothesis $H_a$ make an appearance in the formula for the power?
Let's assume we have a method to compute an estimator $\widehat{\theta}$ and its standard error $\sqrt{\widehat{\text{Var}}(\widehat{\theta})}$ from a dataset. Given $\widehat{\theta}$ and its SE we can do hypothesis testing: compute a p-value and a confidence interval and decide whether to reject or not reject the null hypothesis.
In a simulation study we generate many datasets and repeat the analysis for each dataset. In the $i$th simulation:
$$
\widehat{\theta}_i \pm z_{\alpha/2}\sqrt{\widehat{\text{Var}}(\widehat{\theta}_i)}
$$
is a $(1-\alpha)$100% confidence interval for the parameter $\theta$.
Furthermore, we reject the null hypothesis $H_0$ at the $\alpha$ significance level if $\theta_0$ is outside the $(1-\alpha)$100% confidence interval. Otherwise we don't reject the null. Note that this is a yes-or-no decision.
We record our reject/don't reject decision for each simulation. The proportion of rejections is an estimate of the power P.
$$
\widehat{P} = \frac{1}{n_{sim}} \sum_{i=1}^{n_{sim}} \mathbb{I}\left\{ \theta_0 \notin \left[ \widehat{\theta}_i \pm z_{\alpha/2}\sqrt{\widehat{\text{Var}}(\widehat{\theta}_i)} \right] \right\}
$$
And since this estimate of power is a proportion (of rejections) in a set of events (simulations), its Monte Carlo standard error has the form:
$$
\text{MCSE}(\widehat{P}) = \sqrt{\widehat{P}\left(1 - \widehat{P}\right) / n_{sim}}
$$
Now let's do a simulation study to estimate the power of the Wald test to test the null hypothesis $H_0: \pi = 0.6$ vs the alternative hypothesis $H_a: \pi \neq 0.6$.
Notes: 1. The Wald test is not the most powerful test for proportions. 2. You asked for the one-sided alternative $\pi < 0.6$ but, if I read the documentation correctly, simsum assumes a two-sided test. 3. So we will compute the one-sided power ourselves!
library("rsimsum")
library("tidyverse")
set.seed(1234)
pi0 <- 0.6 # null hypothesis
pi1 <- 0.5 # alternative hypothesis
delta <- pi1 - pi0
sigma <- sqrt(pi1 * (1 - pi1))
# To be sure we know how to do a simulation study with rsimsum,
# let's find out the sample size necessary to achieve 90% power
alpha <- 0.05
power <- 0.9
za <- qnorm(1 - alpha / 2)
zb <- qnorm(power)
delta <- pi1 - pi0
sigma <- sqrt(pi1 * (1 - pi1))
# (Approximate) Sample size to achieve specified power at significance alpha
n <- ((za + zb) * sigma / delta)^2
n <- ceiling(n)
n
#> [1] 263
# We'll do lots of simulations
reps <- 1000
data <-
crossing(
dataset = seq(reps),
n = seq(n)
) %>%
mutate(
# Simulate data
x = rbinom(n(), size = 1, prob = pi1)
) %>%
group_by(
dataset
) %>%
group_modify(
~ summarise(.x,
n = n(),
# Compute estimator
xbar = mean(x),
# Compute its standard error
se = sqrt(xbar * (1 - xbar) / n)
)
)
data
#> # A tibble: 1,000 × 4
#> # Groups: dataset [1,000]
#> dataset n xbar se
#> <int> <int> <dbl> <dbl>
#> 1 1 263 0.498 0.0308
#> 2 2 263 0.548 0.0307
#> 3 3 263 0.521 0.0308
#> 4 4 263 0.498 0.0308
#> 5 5 263 0.468 0.0308
#> 6 6 263 0.464 0.0308
#> 7 7 263 0.475 0.0308
#> 8 8 263 0.513 0.0308
#> 9 9 263 0.521 0.0308
#> 10 10 263 0.532 0.0308
#> # … with 990 more rows
s <- simsum(
data = data,
estvarname = "xbar", se = "se", true = pi1,
by = c("n")
)
tidy(s, stats = c("thetamean", "bias", "cover", "power"))
#> stat est mcse n
#> 1 thetamean 0.4997110266 NA 263
#> 2 bias -0.0002889734 0.0009306672 263
#> 3 cover 0.9570000000 0.0064149045 263
#> 4 power 1.0000000000 0.0000000000 263
The estimate $\widehat{\theta}$ is correct, the coverage $1-\alpha$ is correct and simsum reports the power is 100%. However, we know that the power is 90%.
The issue is that simsum assumes the null hypothesis is $\theta = 0$. To see this, study the formula for the power in your question. And also note while simsum takes an argument true for the true value of the parameter, there is no argument for the null value.
Let's account for this by subtracting $\pi_0$ both from the estimators $\bar{x}$ and from the true value $\pi_1$.
s <- simsum(
data = data %>% mutate(xbar = xbar - pi0),
estvarname = "xbar", se = "se", true = pi1 - pi0,
by = c("n")
)
tidy(s, stats = c("thetamean", "bias", "cover", "power"))
#> stat est mcse n
#> 1 thetamean -0.1002889734 NA 263
#> 2 bias -0.0002889734 0.0009306672 263
#> 3 cover 0.9570000000 0.0064149045 263
#> 4 power 0.8970000000 0.0096120237 263
Now we get that the power is 90% as expected.
Finally, we can work around the limitations of simsum to show that the power for the one-sided hypothesis $H_a: \pi < 0.6$ is 96%. This also confirms we understand how to estimate the Monte Carlo standard error of the power.
data %>%
mutate(
upper_limit = xbar + qnorm(1 - alpha) * se
) %>%
ungroup() %>%
summarise(
power = mean(pi0 > upper_limit)
)
#> # A tibble: 1 × 1
#> power
#> <dbl>
#> 1 0.96
