10.2 - power curves and simulations

library(tidyverse)

Warning: package 'tidyverse' was built under R version 4.3.3

Warning: package 'tibble' was built under R version 4.3.3

Warning: package 'tidyr' was built under R version 4.3.3

Warning: package 'readr' was built under R version 4.3.3

Warning: package 'purrr' was built under R version 4.3.3

Warning: package 'dplyr' was built under R version 4.3.3

Warning: package 'stringr' was built under R version 4.3.3

Warning: package 'forcats' was built under R version 4.3.3

Warning: package 'lubridate' was built under R version 4.3.3

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   4.0.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(purrrfect)


Attaching package: 'purrrfect'

The following objects are masked from 'package:base':

    replicate, tabulate

W1

Suppose \(Y_1,...,Y_n \stackrel{i.i.d.}{\sim} N(\mu,\sigma^2 = 36)\). We want to test:

\[H_0: \mu = 3\] \[H_a: \mu > 3\]

We will reject \(H_0\) using the test statistic:

\[T = \frac{\bar Y -3}{\sigma/\sqrt{n}}\]

and will reject when \(T > c\).

A. Suppose we want to set \(c\) such that \(\alpha = 0.05\). Find \(c\).

B. Find an expression for the power as a function of \(n\) and the “effect size”, \(d=(\mu_a - \mu_0)/\sigma\).

C. Consider an effect size of \(d=1\). Determine the sample size necessary to achieve a power of at least 0.8.

D. Plot the analytic power function as a function of the effect size \(d\) for \(n = 20\).

ggplot() +
  geom_function(fun = \(d) 1 - pnorm(1.645 - sqrt(20)*d)) +
  labs(x = "Effect Size", y = "Rejection Rate") +
  geom_hline(aes(yintercept = .05),linetype = 2) +
  theme_classic()

E. Plot the analytic power function as a function of \(n\) if \(d = 1\).

ggplot() +
  geom_function(fun = \(n) 1 - pnorm(1.645 - sqrt(n)*1)) +
  labs(x = "Sample Size", y = "Power") +
  xlim(c(1,20)) +
  theme_classic()

F. Although not necessary, use a simulation study to approximate the values of \(c\) for each \(n\in \{5,10,...,50\}\). Compare the simulated values of \(c\) to their analytic counterparts.

w1df <- (parameters(~n,seq(5,50,5))
         %>%add_trials(10000)
  %>%mutate(Ysamp = map(n,\(n) rnorm(n,3,6)),T = pmap_dbl(list(Ysamp,n),\(y,n) sqrt(n)*(mean(y) - 3)/6))
  %>%summarize("C" = quantile(T,.95),.by = n)
  
  
)
w1df

# A tibble: 10 × 2
       n     C
   <dbl> <dbl>
 1     5  1.63
 2    10  1.66
 3    15  1.63
 4    20  1.62
 5    25  1.63
 6    30  1.63
 7    35  1.63
 8    40  1.65
 9    45  1.65
10    50  1.62

G. Simulate the rejection rates of this test using the analytic rejection regions for \(\mu\in\{3,3.5,4,4.5\}\). Plot the rejection rates as a function of \(n\) faceted by the effect size \(d\).

w1dg <- (parameters(~n,~mu,seq(5,50,5),c(3,3.5,4,4.5))
         %>%add_trials(10000)
  %>%mutate(Ysamp = pmap(list(n,mu),\(n,m) rnorm(n,m,6)),T = pmap_dbl(list(Ysamp,n),\(y,n) sqrt(n)*(mean(y) - 3)/6))
  %>%inner_join(.,w1df, by = "n")
  %>%mutate(Reject = ifelse(T > C,1,0))
  %>%summarize("Rejection" = mean(Reject), .by = c(n,mu))
)
w1dg

# A tibble: 40 × 3
       n    mu Rejection
   <dbl> <dbl>     <dbl>
 1     5   3      0.0508
 2     5   3.5    0.0697
 3     5   4      0.105 
 4     5   4.5    0.14  
 5    10   3      0.0501
 6    10   3.5    0.0792
 7    10   4      0.132 
 8    10   4.5    0.201 
 9    15   3      0.0484
10    15   3.5    0.0962
# ℹ 30 more rows

ggplot(data = w1dg)+
  geom_line(aes(x = n,y = Rejection)) +
  geom_hline(aes(yintercept = .05), linetype = 2) +
  facet_grid(~mu,labeller = label_both) +
  theme_classic()

Problem 1

Suppose \(Y_1,...,Y_n\) are an i.i.d. \(UNIF(0,\theta)\) sample. Consider two different tests of \(H_0: \theta = 1\) versus \(H_a: \theta > 1\).

Test 1 rejects \(H_0\) if \(Y_{(n)}> c_1\).
Test 2 rejects \(H_0\) if \(\bar Y> c_2\).

We want both tests to have 5% type-I error rate.

A. Show that \(c_1\) can be found analytically to be \(0.95^{1/n}\).

\(f_y(y) = \frac{1}{\theta}\) and \(F_y(y) = \frac{y}{\theta}\) thus \(F_n(y) = F_y(y)^n = \frac{y^n}{\theta^n}.\)
\(\alpha = .05 = P_0(Y_{(n)}>c_1)=1-P_0(Y_{(n)}<c_1) \Rightarrow .05 = 1 - c_1^n \Rightarrow .95^{1/n} = c_1\)

Find a general expression for the power function of Test 1 as a function of \(n\) and \(\theta_a\). Plot curve as a function of \(n\) for \(\theta = 1.2.\)

\(Power = P_a(Y_{(n)}>.95^{1/n}) = 1- \left(\frac{.95^{1/n}}{\theta_a}\right)^n = 1-.95(\theta_1)^{-n}\)

ggplot() +
  geom_function(fun = \(n) 1-.95*(1/1.2^n)) +
  xlim(1,100) +
  labs(y = "Power",x = "n") +
  theme_classic()

C. Explain why an exact value of \(c_2\) cannot be found analytically and must be simulated.
Unsure about distribution of \(\bar{Y}\) , as normally could find through MGF method, but \(M_Y(t) = \frac{e^{t\theta}-1}{t\theta} \Rightarrow M_{\bar{Y}}(t) = M_Y(\frac{t}{n})^n = (\frac{e^{\frac{t}{n}\theta}-1}{\frac{t}{n}\theta})^n\)

So we do not have a well known form for distribution of \(\bar{Y}\), and thus cannot find \(c_2\).

D. Simulate values of \(c_2\) that yield 5% type-I error rates for Test 2, for \(n \in \{5,10,...,50\}\).

RRsim <- (parameters(~n,seq(5,50,5))
          %>%add_trials(10000)
          %>%mutate(Yn = map(n,\(n) runif(n,0,1)), Yhat = map_dbl(Yn,mean))
          %>%summarize("c_2sim" = quantile(Yhat,.95),.by = n)
)

E. Use a simulation study to compare the sizes and the power curves of the two tests as a function of \(n\) for \(\theta_{true} \in \{1,1.1,1.2,1.3\}.\) Use analytic rejection regions for \(Y_{(n)}\) and simulated rejection regions for \(\bar Y\), but simulate the rejection probabilities for both. Plot the simulated Type-I error rate/power as a function of \(n\) for each \(\theta_{true}\). Discuss the properties of these tests: does one seem preferable?

Hsim <- (parameters(~n,~theta,seq(5,50,5),c(1,1.1,1.2,1.3))
         %>%add_trials(10000)
         %>%inner_join(., RRsim,by = 'n')
         %>%mutate(Yn = pmap(list(n,theta),\(n,t) runif(n,0,t)),Ymax = map_dbl(Yn,max),Ybar = map_dbl(Yn,mean),YnD = ifelse(Ymax > .95^(1/n),1,0),YmaxD = ifelse((Ybar > c_2sim),1,0))
         %>%summarize("PowerYn" = mean(YnD),
                      "PowerYbar" = mean(YmaxD),
                      .by = c(n,theta))
)

ggplot(data = Hsim) +
  geom_line(aes(x = n, y =PowerYn,col = "Max Test")) +
  geom_line(aes(x= n, y = PowerYbar,col = "Mean Test")) +
  facet_wrap(~theta,labeller = label_bquote(theta == .(theta))) +
  labs(y = "Null Rejection %",col = "") +
  theme_classic()

It seems that the test using \(Y_{(n)}\) is more preferable as for all \(\theta_{true}\) it has a higher Null Rejection %, i.e power with the same data.