10.2 - power curves and simulations

Power curves (aka power functions)

To define power curves, we define the following:

\[\Omega_0 = \mbox{The set of possible values of the parameter(s) under the null}\] \[\Omega_a = \mbox{The set of possible values of the parameter(s) under the alternative}\]

The power curve, or a power function, of a test is one of the following:
- A plot of the power versus the sample size \(n\), for a fixed \(\theta \in \Omega_a\);
- A plot of the power versus \(\Omega_a\), for a fixed \(n\).

Exponential rate revisited

Recall sample of \(n=15\) i.i.d. \(EXP(\lambda)\) observations to test:

\[H_0: \lambda = 1/5\] \[H_a: \lambda < 1/5\]

\(\Omega_0: \{0.2\}\)
\(\Omega_a = (0, 0.2)\)
Size-5% test: reject when \(Y_{(1)} >0.9985774\)

Power curve

The curve below is an analytic power curve; computed using analytic probabilities:

library(tidyverse)
ggplot() + 
  geom_function(fun=\(lam) 1-pexp(0.9985774, rate = 15*lam)) + 
  xlim(c(0.001, 1/5)) + 
  geom_hline(aes(yintercept = 0.05), linetype=2)+ 
  scale_y_continuous(breaks = c(0.05, seq(0.2,1,by=.2)))+
  labs(y='Power', x = expression(lambda))+
  theme_classic(base_size=14)

Simulating type-I error and power

To simulate type-I error and power, we replicate the process of:
- Simulating a sample;
- Forming the test statistic;
- Determining whether the test statistic falls in the rejection region.
Once we have our replications, we summarize to find the percent of replications where we reject \(H_0\).

Simulating exponential example

Recall sample of \(n=15\) i.i.d. \(EXP(\lambda)\) observations to test:

\[H_0: \lambda = 1/5\] \[H_a: \lambda < 1/5\]

Reject when \(Y_{(1)}>0.9985774\).

library(purrrfect)
n <- 15
lambda_0 = 1/5
crit_value <- qexp(0.95, rate = n*lambda_0)
N <- 10000

(parameters(~lambda, seq(0.001, 0.2, l=100))
  %>% add_trials(N)
  %>% mutate(Ysample = map(lambda, \(l) rexp(n, rate = l)),
             Y1 = map_dbl(Ysample, min),
             reject = ifelse(Y1 > crit_value, 1, 0)
             )
  ) -> simulated_decisions

(simulated_decisions
   %>% summarize(percent_reject = mean(reject), .by = lambda)
) -> power_and_type1

head(power_and_type1,3)

# A tibble: 3 × 2
   lambda percent_reject
    <dbl>          <dbl>
1 0.001            0.985
2 0.00301          0.956
3 0.00502          0.931

tail(power_and_type1,3)

# A tibble: 3 × 2
  lambda percent_reject
   <dbl>          <dbl>
1  0.196         0.0494
2  0.198         0.0489
3  0.2           0.05

Plotting simulated power curve

ggplot(data = power_and_type1) + 
  geom_line(aes(x=lambda, y = percent_reject)) + 
  geom_hline(aes(yintercept = 0.05), linetype=2)+ 
  scale_y_continuous(breaks = c(0.05, seq(0.2,1,by=.2)))+
  labs(y='Proportion of replications rejecting null', x = expression(lambda[TRUE]))+
  theme_classic(base_size=14)

Simulating critical values

In the previous example, we were able to define analytic critical value for \(Y_{(1)}\)
- How? We knew the sampling distribution of \(Y_{(1)}\)!
In some cases, we won’t know the sampling distribution of our test statistic \(\Rightarrow\) cannot determine analytic critical value/RR.
Simulating critical values requires simulating the sampling distribution of the test statistic under \(H_0\) and finding the quantile corresponding to the desired significance level (often 0.05).
Process: simulate replications of test statistic(s) under \(H_0\) for various \(n\), summarize.

Finding RRs for exponential min and median

As an example, reconsider testing \(H_0: \lambda = 1/5\) vs \(H_a: \lambda < 1/5\) with i.i.d. \(EXP(\lambda)\) sample.
Will test using sample of \(n\) i.i.d. observations

Test 1: reject when \(T_1 = Y_{(1)} > c_1\) \(\Leftarrow\) can find analytically, since \(Y_{(1)}\sim EXP(n\lambda_0)\) under \(H_0\); will depend on \(n\).
- \(\alpha = P(Y_{(1)}> c_1| \lambda = 1/5) = e^{-n\cdot (1/5)\cdot c} \stackrel{set}{=}0.05 \Rightarrow c_1 = -5\ln(0.05)/n\)

Test 2: reject when \(T_2 = median > c_2\) \(\Leftarrow\) cannot find analytically; don’t know distribution of median! Must find via simulation.

Note that both critical values will depend on the sample size \(n\).

Finding RRs for exponential min and median

lam0 <- 1/5

(exp_crit_vals <- parameters(~n ,
            seq(10,100,by=5)
            )
  %>% add_trials(1000)
  %>% mutate(Ysample = map(n,.f = \(n) rexp(n, rate = lam0)),
             T1 = map_dbl(Ysample, min), # <- unnecessary, but could use to get simulated c1 if needed.
             T2 = map_dbl(Ysample, median)
             )
  %>% summarize(c1_simulated = quantile(T1, 0.95),
                c2_simulated = quantile(T2, 0.95), 
                .by = n
                )
  %>% mutate(c1_analytic = -5*log(0.05)/n)
)  %>% head

# A tibble: 6 × 4
      n c1_simulated c2_simulated c1_analytic
  <dbl>        <dbl>        <dbl>       <dbl>
1    10        1.39          6.50       1.50 
2    15        0.969         6.11       0.999
3    20        0.765         5.51       0.749
4    25        0.532         5.43       0.599
5    30        0.478         5.17       0.499
6    35        0.446         5.10       0.428

Simulating decisions using simulated RRs

To simulate decisions with simulated RRs, we must:
- Replicate test statistic(s) under \(H_a\) parameter grid;
- Join replicates with critical values by \(n\) \(\Leftarrow\) new step;
- Make decisions;
- Summarize.

Exponential power simulation

We will now simulate \(T_1\) and \(T_2\) using \(\lambda \in \Omega_0 \cup\Omega_a\); must use same \(n\) grid as critical values!
- Will use analytic \(c_1\) to make decisions using \(T_1\)
- Will use simulated \(c_2\) to make decisions for \(T_2\)

options(width = 10000)
(exp_decisions <- parameters(~n , ~true_lambda, 
            seq(10,100,by=5), seq(0.05,0.2, by = 0.05)
            )
  %>% add_trials(1000)
  %>% inner_join(., exp_crit_vals, by ='n') # <- important! join critical values to data set 
  %>% mutate(Ysample = pmap(list(n,true_lambda),.f = \(n,l) rexp(n, rate = l)),
             T1 = map_dbl(Ysample, min),
             T2 = map_dbl(Ysample, median),
             T1_decision = ifelse(T1 > c1_analytic, 1,0),
             T2_decision = ifelse(T2 > c2_simulated, 1,0)
             )
)  %>% head()

# A tibble: 6 × 11
      n true_lambda .trial c1_simulated c2_simulated c1_analytic Ysample       T1    T2 T1_decision T2_decision
  <dbl>       <dbl>  <dbl>        <dbl>        <dbl>       <dbl> <list>     <dbl> <dbl>       <dbl>       <dbl>
1    10        0.05      1         1.39         6.50        1.50 <dbl [10]> 2.13   13.7           1           1
2    10        0.05      2         1.39         6.50        1.50 <dbl [10]> 4.41   25.6           1           1
3    10        0.05      3         1.39         6.50        1.50 <dbl [10]> 0.411  14.0           0           1
4    10        0.05      4         1.39         6.50        1.50 <dbl [10]> 0.949  13.9           0           1
5    10        0.05      5         1.39         6.50        1.50 <dbl [10]> 2.50   14.2           1           1
6    10        0.05      6         1.39         6.50        1.50 <dbl [10]> 0.777  11.3           0           1

Exponential power simulation

Finally, summarize to find Power:

(simulated_prob_reject <- exp_decisions
 %>% summarize(power_T1 = mean(T1_decision),
               power_T2 = mean(T2_decision),
               .by=c(n,true_lambda)
               )
) %>% head

# A tibble: 6 × 4
      n true_lambda power_T1 power_T2
  <dbl>       <dbl>    <dbl>    <dbl>
1    10        0.05    0.481    0.947
2    10        0.1     0.21     0.582
3    10        0.15    0.098    0.19 
4    10        0.2     0.045    0.048
5    15        0.05    0.44     0.974
6    15        0.1     0.248    0.63

Plotting power curves vs \(n\)

ggplot(data = simulated_prob_reject) + 
  geom_line(aes(x=n, y = power_T1, col='Test using min')) + 
  geom_line(aes(x=n, y = power_T2, col='Test using median')) + 
  scale_y_continuous(breaks = c(0.05, seq(0.2,1,by=.2)))+
  geom_hline(aes(yintercept = 0.05), linetype=2)+ 
  labs(y='Proportion of replications rejecting null', 
       x='n',
       color='')+
  facet_wrap(~true_lambda, labeller = label_bquote(lambda == .(true_lambda)))+
  theme_classic(base_size=14)