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\).
E. Plot the analytic power function as a function of \(n\) if \(d = 1\).
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.
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\).
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\).
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}\).
n_vals <- seq(5, 50, by = 5)
c1_vals <- 0.95^(1 / n_vals)
data.frame(n = n_vals, c1 = c1_vals) n c1
1 5 0.9897938
2 10 0.9948838
3 15 0.9965863
4 20 0.9974386
5 25 0.9979504
6 30 0.9982917
7 35 0.9985356
8 40 0.9987185
9 45 0.9988608
10 50 0.9989747n_vals <- seq(5, 50, by = 5)
power_test1_theta12 <- 1 - 0.95 / (1.2^n_vals)
plot( n_vals, power_test1_theta12, type = "b", pch = 19, xlab = "n", ylab = "Power", main = expression("Power of Test 1 for " * theta == 1.2) )
C. Explain why an exact value of \(c_2\) cannot be found analytically and must be simulated.
c2 is harder to get exactly because it depends on the distribution of YˉYYˉ, not the maximum. Under H0H_0H0, YˉYYˉ is the average of uniform random variables, and there is not a simple formula to directly solve for the 95th percentile for each nnn. Because of that, simulation is used to estimate the cutoff value that gives a 5% type-I error rate
D. Simulate values of \(c_2\) that yield 5% type-I error rates for Test 2, for \(n \in \{5,10,...,50\}\).
n_vals <- seq(5, 50, by = 5)
B <- 50000
c2_vals <- numeric(length(n_vals))
for (j in seq_along(n_vals)) {
n <- n_vals[j]
ybar_sim <- replicate(B, mean(runif(n, 0, 1)))
c2_vals[j] <- quantile(ybar_sim, 0.95)
}
data.frame(n = n_vals, c2 = c2_vals) n c2
1 5 0.7119813
2 10 0.6492318
3 15 0.6213370
4 20 0.6071735
5 25 0.5952204
6 30 0.5860735
7 35 0.5801063
8 40 0.5752111
9 45 0.5704200
10 50 0.5671412E. 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?
theta_vals <- c(1, 1.1, 1.2, 1.3)
n_vals <- seq(5, 50, by = 5)
c1_vals <- 0.95^(1 / n_vals)
B <- 20000
results <- data.frame()
for (theta_true in theta_vals) {
for (j in seq_along(n_vals)) {
n <- n_vals[j]
c1 <- c1_vals[j]
c2 <- c2_vals[j]
reject1 <- replicate(B, {
y <- runif(n, 0, theta_true)
max(y) > c1
})
reject2 <- replicate(B, {
y <- runif(n, 0, theta_true)
mean(y) > c2
})
results <- rbind(
results,
data.frame(
theta_true = theta_true,
n = n,
test = "Test1",
rejection_probability = mean(reject1)
),
data.frame(
theta_true = theta_true,
n = n,
test = "Test2",
rejection_probability = mean(reject2)
)
)
}
}
results theta_true n test rejection_probability
1 1.0 5 Test1 0.04915
2 1.0 5 Test2 0.04955
3 1.0 10 Test1 0.05165
4 1.0 10 Test2 0.05225
5 1.0 15 Test1 0.05035
6 1.0 15 Test2 0.05185
7 1.0 20 Test1 0.05120
8 1.0 20 Test2 0.04810
9 1.0 25 Test1 0.04950
10 1.0 25 Test2 0.04765
11 1.0 30 Test1 0.04710
12 1.0 30 Test2 0.05110
13 1.0 35 Test1 0.05075
14 1.0 35 Test2 0.04860
15 1.0 40 Test1 0.04960
16 1.0 40 Test2 0.04750
17 1.0 45 Test1 0.05165
18 1.0 45 Test2 0.05085
19 1.0 50 Test1 0.05030
20 1.0 50 Test2 0.04780
21 1.1 5 Test1 0.40960
22 1.1 5 Test2 0.13475
23 1.1 10 Test1 0.62510
24 1.1 10 Test2 0.16160
25 1.1 15 Test1 0.77195
26 1.1 15 Test2 0.19930
27 1.1 20 Test1 0.85735
28 1.1 20 Test2 0.21105
29 1.1 25 Test1 0.91320
30 1.1 25 Test2 0.24320
31 1.1 30 Test1 0.94420
32 1.1 30 Test2 0.26930
33 1.1 35 Test1 0.96395
34 1.1 35 Test2 0.28440
35 1.1 40 Test1 0.97850
36 1.1 40 Test2 0.31025
37 1.1 45 Test1 0.98670
38 1.1 45 Test2 0.33220
39 1.1 50 Test1 0.99255
40 1.1 50 Test2 0.35215
41 1.2 5 Test1 0.61665
42 1.2 5 Test2 0.24180
43 1.2 10 Test1 0.84800
44 1.2 10 Test2 0.32875
45 1.2 15 Test1 0.94130
46 1.2 15 Test2 0.40410
47 1.2 20 Test1 0.97465
48 1.2 20 Test2 0.46320
49 1.2 25 Test1 0.98985
50 1.2 25 Test2 0.52180
51 1.2 30 Test1 0.99550
52 1.2 30 Test2 0.58215
53 1.2 35 Test1 0.99840
54 1.2 35 Test2 0.63635
55 1.2 40 Test1 0.99910
56 1.2 40 Test2 0.67200
57 1.2 45 Test1 0.99990
58 1.2 45 Test2 0.71280
59 1.2 50 Test1 0.99990
60 1.2 50 Test2 0.74685
61 1.3 5 Test1 0.74615
62 1.3 5 Test2 0.35965
63 1.3 10 Test1 0.92890
64 1.3 10 Test2 0.50630
65 1.3 15 Test1 0.98160
66 1.3 15 Test2 0.61885
67 1.3 20 Test1 0.99495
68 1.3 20 Test2 0.69810
69 1.3 25 Test1 0.99860
70 1.3 25 Test2 0.76695
71 1.3 30 Test1 0.99990
72 1.3 30 Test2 0.82520
73 1.3 35 Test1 0.99995
74 1.3 35 Test2 0.86820
75 1.3 40 Test1 0.99995
76 1.3 40 Test2 0.89300
77 1.3 45 Test1 1.00000
78 1.3 45 Test2 0.92285
79 1.3 50 Test1 1.00000
80 1.3 50 Test2 0.93930par(mfrow = c(2, 2))
for (th in theta_vals) {
subdat <- subset(results, theta_true == th)
y1 <- subset(subdat, test == "Test1")$rejection_probability
y2 <- subset(subdat, test == "Test2")$rejection_probability
plot(
n_vals, y1,
type = "b", pch = 19, ylim = c(0, 1),
xlab = "n",
ylab = "Rejection Probability",
main = bquote(theta == .(th))
)
lines(n_vals, y2, type = "b", pch = 17, lty = 2)
legend(
"bottomright",
legend = c("Test1", "Test2"),
pch = c(19, 17),
lty = c(1, 2),
bty = "n"
)
}
At θ=1=1θ=1, both tests should stay close to 0.05 since that is the target size. As θincreases, both rejection probabilities increase, but Test 1 usually increases faster because the maximum reacts more directly when the upper endpoint gets larger.