nrm_pwr1 <- function(n = 30,
alpha = .05,
H0 = 0,
H1 = 2,
sig = 1) {
# Find the z_critical value
z_alpha <- qnorm(p = alpha,
mean = 0,
sd = 1,
lower.tail = FALSE)
# Find the ybar_critical value
y_bar_cv <- z_alpha*(sig/sqrt(n)) + H0
# Calculate the power under H1.
power <- pnorm(q = y_bar_cv,
mean = H1,
sd = sig/sqrt(n),
lower.tail = FALSE)
# Report the power.
power
}
### TASK 1
##########################################
### Finish the plot by adding curves of different colors
### and a legend based
### on varying alpha: .001, .005, .01, .05, .1
#########################################
pcurv.05 <- lapply(X = 1:200,
FUN = nrm_pwr1,
alpha = .05, H0 = 100, H1 = 105, sig = 15)
plot(x = 1:200, y = pcurv.05, type = "l",
xlab = "Sample Size (n)", ylim = c(0,1),
ylab = expression(paste("Power When ", mu, " = 105")),
col = 3, lwd = 2)
pcurv.01 <- lapply(X = 1:200,
FUN = nrm_pwr1,
alpha = .01, H0 = 100, H1 = 105, sig = 15)
points(x = 1:200, y = pcurv.01,
type = "l", col = 4, lwd = 2)
pcurv.1 <- lapply(X = 1:200,
FUN = nrm_pwr1,
alpha = .1, H0 = 100, H1 = 105, sig = 15)
points(x = 1:200, y = pcurv.1,
type = "l", col = 5, lwd = 2)
pcurv.005 <- lapply(X = 1:200,
FUN = nrm_pwr1,
alpha = .005, H0 = 100, H1 = 105, sig = 15)
points(x = 1:200, y = pcurv.005,
type = "l", col = 6, lwd = 2)
pcurv.001 <- lapply(X = 1:200,
FUN = nrm_pwr1,
alpha = .001, H0 = 100, H1 = 105, sig = 15)
points(x = 1:200, y = pcurv.001,
type = "l", col = 7, lwd = 2)
legend(x = 120, y = .4,
legend = c("alpha = .1", "alpha = .05","alpha = .01","alpha = .005","alpha = .001"),
lwd = 2, col = c(5,4,3,6,7))
grid(col = gray(.3))
### TASK 2
##########################################
### Use the function to create a set of power curves
### based on varying sigma instead of n. Use a fixed
### value of n = 30. Make curves for alpha = .01, .05.
### Let sigma range from 1 to 30.
#########################################
pcurvSig.01 <- lapply(X = 1:30,
FUN = nrm_pwr1,
n=30, alpha = .01, H0 = 100, H1 = 105)
pcurvSig.05 <- lapply(X = 1:30,
FUN = nrm_pwr1,
n=30, alpha = .05, H0 = 100, H1 = 105)
plot(x = 1:30, y = pcurvSig.01, type = "l",
xlab = expression(paste(sig)), ylim = c(0,1),
ylab = expression(paste("Power When n = 30")),
col = 3, lwd = 2)
points(x = 1:30, y = pcurvSig.05,
type = "l", col = 5, lwd = 2)
legend(x = 22, y = 1,
legend = c("alpha = .05", "alpha = .01"),
lwd = 2, col = c(3,5))
grid(col = gray(.3))
### TASK 3
##########################################
### Use the function to create a set of power curves
### based on varying the value of the alternative
### hypothesis from 100 to 120. Do this for
### sample sizes: n = 10, 20, 30 on the same plot.
#########################################
pcurvH1.10 <- lapply(X = 100:120,
FUN = nrm_pwr1,
n=10, alpha = .05, H0 = 100, sig =15)
pcurvH1.20 <- lapply(X = 100:120,
FUN = nrm_pwr1,
n=20, alpha = .05, H0 = 100, sig =15)
pcurvH1.30 <- lapply(X = 100:120,
FUN = nrm_pwr1,
n=30, alpha = .05, H0 = 100, sig =15)
plot(x = 100:120, y = pcurvH1.10, type = "l",
xlab = expression(paste(H1)), ylim = c(0,1),
ylab = expression(paste("Power When", sigma, "= 15")),
col = 3, lwd = 2)
points(x = 100:120, y = pcurvH1.20,
type = "l", col = 4, lwd = 2)
points(x = 100:120, y = pcurvH1.30,
type = "l", col = 5, lwd = 2)
legend(x = 116, y = .2,
legend = c("n=10", "n=20","n=30"),
lwd = 2, col = c(3,4,5))
grid(col = gray(.3))
### TASK 4
##########################################
### Summarize the relationships between power and
### (a) sample size,
### (b) variance,
### (c) effect size (value of H1),
### (d) Type I error rate.
### For example, for part (a), assuming
### a fixed variance, fixed effect
### size, and fixed Type I error rate, what
### happens to power as the sample size is
### increased?
##########################################
# (a) when sample size goes up, power increases
# (b) when variance increase, power decreases
# (c) when value of H1 increase, power also increases
# (d) Type I error rate increases, power also increases
### TASK 5 Bonus
##########################################
### Write a new function called rnr_pwr2 that
### calculates power for the two-sided test.
#########################################
rnr_pwr2 <- function(n = 30,
alpha = .05,
H0 = 0,
H1 = 2,
sig = 1) {
# Find the z_critical value
z_alpha_half <- qnorm(p = alpha/2,
mean = 0,
sd = 1,
lower.tail = FALSE)
# Find the ybar_critical value
y_bar_cv1 <- z_alpha_half*(sig/sqrt(n)) + H0
y_bar_cv2 <- -z_alpha_half*(sig/sqrt(n)) + H0
# Calculate the power under H1.
power <- pnorm(q = y_bar_cv1,
mean = H1,
sd = sig/sqrt(n),
lower.tail = FALSE)+
pnorm(q = y_bar_cv2,
mean = H1,
sd = sig/sqrt(n),
lower.tail = TRUE)
# Report the power.
power
}