# Exercises for 'Hypothesis Testing' chapter
```{r}
#| include: false
# if it is available, run the setup script that tells quarto to round all df/tibble outputs to three decimal places
if(file.exists("../_setup.R")){source("../_setup.R")}
```
```{r}
#| include: false
# dependencies
library(tidyr)
library(dplyr)
library(readr)
library(furrr)
library(parameters)
library(simhelpers)
library(stringr)
library(forcats)
library(ggplot2)
library(scales)
library(ggstance)
library(patchwork)
library(janitor)
library(effectsize)
library(ggstance)
library(knitr)
library(kableExtra)
# set up parallelization
# use availableCores() to define the number of cores, in case the simulation is being
# run on a HPC that might have more physical cores than are available to the session
plan(multisession, workers = parallelly::availableCores())
```
## Reading
Read [Lakens et al. (2018) "Equivalence Testing for Psychological Research: A Tutorial"](https://doi.org/10.1177/2515245918770963).
## Basis simulation
```{r}
# functions for simulation
generate_data <- function(n_per_condition,
mean_control,
mean_intervention,
sd) {
data_control <-
tibble(condition = "control",
score = rnorm(n = n_per_condition, mean = mean_control, sd = sd))
data_intervention <-
tibble(condition = "intervention",
score = rnorm(n = n_per_condition, mean = mean_intervention, sd = sd))
data_combined <-
bind_rows(data_control,
data_intervention) |>
mutate(condition = factor(condition, level = c("intervention", "control")))
return(data_combined)
}
analyse <- function(data, direction_h1 = "two.sided", alpha = 0.05, mean_difference_h0 = 0){
fit <- t.test(
formula = score ~ condition,
mu = mean_difference_h0,
alternative = direction_h1,
conf.level = 1-alpha,
var.equal = TRUE, # students t test: assumes equal variances
data = data
)
results <- fit %>%
model_parameters() %>%
as_tibble() %>%
janitor::clean_names() %>%
# create estimand
mutate(significant_p = p < .05, # statistically significant in the sense that p < .05
signficiant_ci = !dplyr::between(0, ci_low, ci_high)) |> # statistically significant in the sense that 95% CI excludes a mean difference of 0
# select columns of interest
select(mean_diff = difference,
ci_low,
ci_high,
signficiant_ci,
p,
significant_p)
return(results)
}
# simulation parameters
experiment_parameters <- expand_grid(
n_per_condition = seq(from = 25, to = 100, by = 25),
mean_control = 0,
mean_intervention = c(0, 0.5),
sd = 1,
direction_h1 = "two.sided",
alpha = 0.05,
mean_difference_h0 = 0,
iteration = 1:1000
) |>
mutate(pop_mean_diff = mean_intervention - mean_control)
# run simulation
# results are cached to disk: the simulation is only run the first time, and read back
# in on subsequent renders. delete the .rds file if you change the simulation code.
dir.create("results", showWarnings = FALSE)
if(file.exists("results/simulation_hypothesis_testing_exercises_basis.rds")){
simulation <- read_rds("results/simulation_hypothesis_testing_exercises_basis.rds")
} else {
# set seed
set.seed(42)
simulation <- experiment_parameters |>
# shuffle the rows of the grid to load balance the computationally expensive
# iterations evenly across the workers
slice_sample(prop = 1) |>
mutate(generated_data = future_pmap(.l = list(n_per_condition = n_per_condition,
mean_control = mean_control,
mean_intervention = mean_intervention,
sd = sd),
.f = generate_data,
.progress = TRUE,
.options = furrr_options(seed = TRUE))) |>
mutate(results = future_pmap(.l = list(data = generated_data),
.f = analyse,
.progress = TRUE,
.options = furrr_options(seed = TRUE))) |>
unnest(results) |>
# undo the shuffle so the rows are back in the order of the experiment grid
arrange(n_per_condition, pop_mean_diff, iteration) |>
# the generated data sets are not needed after they have been analyzed, and keeping
# them would make the saved file enormous
select(-generated_data)
write_rds(x = simulation, file = "results/simulation_hypothesis_testing_exercises_basis.rds", compress = "gz")
}
# calculate summary
results_summary <- simulation |>
group_by(pop_mean_diff,
n_per_condition) |>
reframe(mean_difference_in_means = mean(mean_diff),
mean_ci_low = mean(ci_low),
mean_ci_high = mean(ci_high),
# rather than calculating the proportion of significant results by hand, use
# {simhelpers}'s dedicated function for the empirical rejection rate.
# calc_rejection() also returns the number of iterations and the Monte Carlo
# Standard Error of the rate; select() keeps and renames just the rate itself
calc_rejection(data = pick(everything()),
p_values = p,
alpha = 0.05) |>
select(proportion_significant_p = rej_rate))
results_summary |>
mutate_if(is.numeric, round, digits = 2)
```
## Exercises
For each of the exercises below, create a copy of the basis simulation above. Modify it appropriately to answer the exercise.
### Statistical power for a classic Null Hypothesis Significance Test (Two-Sided)
When using a classic Null Hypothesis Significance Test (independent t-test testing a population difference in means (mu) of 0), what is the statistical power to detect a population mean difference of 1, 2, or 3 across different sample sizes (n per condition = 25, 50, 75, 100)?
Notes:
- We're only interested in statistical power here, so set the population difference in means to zero. Check you understand why.
```{r}
# functions for simulation
generate_data <- function(n_per_condition,
mean_control,
mean_intervention,
sd) {
data_control <-
tibble(condition = "control",
score = rnorm(n = n_per_condition, mean = mean_control, sd = sd))
data_intervention <-
tibble(condition = "intervention",
score = rnorm(n = n_per_condition, mean = mean_intervention, sd = sd))
data_combined <-
bind_rows(data_control,
data_intervention) |>
mutate(condition = factor(condition, level = c("intervention", "control")))
return(data_combined)
}
analyse <- function(data, direction_h1 = "two.sided", alpha = 0.05, mean_difference_h0 = 0){
fit <- t.test(
formula = score ~ condition,
mu = mean_difference_h0,
alternative = direction_h1,
conf.level = 1-alpha,
var.equal = TRUE, # students t test: assumes equal variances
data = data
)
results <- fit %>%
model_parameters() %>%
as_tibble() %>%
janitor::clean_names() %>%
# create estimand
mutate(significant_p = p < .05, # statistically significant in the sense that p < .05
signficiant_ci = !dplyr::between(0, ci_low, ci_high)) |> # statistically significant in the sense that 95% CI excludes a mean difference of 0
# select columns of interest
select(mean_diff = difference,
ci_low,
ci_high,
signficiant_ci,
p,
significant_p)
return(results)
}
# simulation parameters
experiment_parameters <- expand_grid(
n_per_condition = seq(from = 25, to = 100, by = 25),
mean_control = 0,
mean_intervention = c(1, 2, 3),
sd = 1,
direction_h1 = "two.sided",
alpha = 0.05,
mean_difference_h0 = 0,
iteration = 1:1000
) |>
mutate(pop_mean_diff = mean_intervention - mean_control)
# run simulation
# results are cached to disk: the simulation is only run the first time, and read back
# in on subsequent renders. delete the .rds file if you change the simulation code.
dir.create("results", showWarnings = FALSE)
if(file.exists("results/simulation_hypothesis_testing_exercises_nhst.rds")){
simulation <- read_rds("results/simulation_hypothesis_testing_exercises_nhst.rds")
} else {
# set seed
set.seed(42)
simulation <- experiment_parameters |>
# shuffle the rows of the grid to load balance the computationally expensive
# iterations evenly across the workers
slice_sample(prop = 1) |>
mutate(generated_data = future_pmap(.l = list(n_per_condition = n_per_condition,
mean_control = mean_control,
mean_intervention = mean_intervention,
sd = sd),
.f = generate_data,
.progress = TRUE,
.options = furrr_options(seed = TRUE))) |>
mutate(results = future_pmap(.l = list(data = generated_data),
.f = analyse,
.progress = TRUE,
.options = furrr_options(seed = TRUE))) |>
unnest(results) |>
# undo the shuffle so the rows are back in the order of the experiment grid
arrange(n_per_condition, pop_mean_diff, iteration) |>
# the generated data sets are not needed after they have been analyzed, and keeping
# them would make the saved file enormous
select(-generated_data)
write_rds(x = simulation, file = "results/simulation_hypothesis_testing_exercises_nhst.rds", compress = "gz")
}
# calculate summary
results_summary <- simulation |>
group_by(pop_mean_diff,
n_per_condition) |>
reframe(mean_difference_in_means = mean(mean_diff),
mean_ci_low = mean(ci_low),
mean_ci_high = mean(ci_high),
# rather than calculating the proportion of significant results by hand, use
# {simhelpers}'s dedicated function for the empirical rejection rate.
# calc_rejection() also returns the number of iterations and the Monte Carlo
# Standard Error of the rate; select() keeps and renames just the rate itself
calc_rejection(data = pick(everything()),
p_values = p,
alpha = 0.05) |>
select(proportion_significant_p = rej_rate))
results_summary |>
mutate_if(is.numeric, round, digits = 2)
```
### Statistical power for a Minimal-Effects Test (One-Sided)
When using a Minimal-Effects Test (One-Sided) (independent t-test testing a population difference in means (mu) of 1), what is the statistical power to detect a population mean difference of 1, 2, or 3 across different sample sizes (n per condition = 25, 50, 75, 100)? What is the false positive rate across these sample sizes?
```{r}
# functions for simulation
generate_data <- function(n_per_condition,
mean_control,
mean_intervention,
sd) {
data_control <-
tibble(condition = "control",
score = rnorm(n = n_per_condition, mean = mean_control, sd = sd))
data_intervention <-
tibble(condition = "intervention",
score = rnorm(n = n_per_condition, mean = mean_intervention, sd = sd))
data_combined <-
bind_rows(data_control,
data_intervention) |>
mutate(condition = factor(condition, level = c("intervention", "control")))
return(data_combined)
}
analyse <- function(data, direction_h1 = "two.sided", alpha = 0.05, mean_difference_h0 = 0){
fit <- t.test(
formula = score ~ condition,
mu = mean_difference_h0,
alternative = direction_h1,
conf.level = 1-alpha,
var.equal = TRUE, # students t test: assumes equal variances
data = data
)
results <- fit %>%
model_parameters() %>%
as_tibble() %>%
janitor::clean_names() %>%
# create estimand
mutate(significant_p = p < .05, # statistically significant in the sense that p < .05
signficiant_ci = !dplyr::between(0, ci_low, ci_high)) |> # statistically significant in the sense that 95% CI excludes a mean difference of 0
# select columns of interest
select(mean_diff = difference,
ci_low,
ci_high,
signficiant_ci,
p,
significant_p)
return(results)
}
# simulation parameters
experiment_parameters <- expand_grid(
n_per_condition = seq(from = 25, to = 100, by = 25),
mean_control = 0,
mean_intervention = c(1, 2, 3),
sd = 1,
direction_h1 = "greater",
alpha = 0.05,
mean_difference_h0 = 1,
iteration = 1:1000
) |>
mutate(pop_mean_diff = mean_intervention - mean_control)
# run simulation
# results are cached to disk: the simulation is only run the first time, and read back
# in on subsequent renders. delete the .rds file if you change the simulation code.
dir.create("results", showWarnings = FALSE)
if(file.exists("results/simulation_hypothesis_testing_exercises_minimal_effects.rds")){
simulation <- read_rds("results/simulation_hypothesis_testing_exercises_minimal_effects.rds")
} else {
# set seed
set.seed(42)
simulation <- experiment_parameters |>
# shuffle the rows of the grid to load balance the computationally expensive
# iterations evenly across the workers
slice_sample(prop = 1) |>
mutate(generated_data = future_pmap(.l = list(n_per_condition = n_per_condition,
mean_control = mean_control,
mean_intervention = mean_intervention,
sd = sd),
.f = generate_data,
.progress = TRUE,
.options = furrr_options(seed = TRUE))) |>
mutate(results = future_pmap(.l = list(data = generated_data),
.f = analyse,
.progress = TRUE,
.options = furrr_options(seed = TRUE))) |>
unnest(results) |>
# undo the shuffle so the rows are back in the order of the experiment grid
arrange(n_per_condition, pop_mean_diff, iteration) |>
# the generated data sets are not needed after they have been analyzed, and keeping
# them would make the saved file enormous
select(-generated_data)
write_rds(x = simulation, file = "results/simulation_hypothesis_testing_exercises_minimal_effects.rds", compress = "gz")
}
# calculate summary
results_summary <- simulation |>
group_by(pop_mean_diff,
n_per_condition) |>
reframe(mean_difference_in_means = mean(mean_diff),
mean_ci_low = mean(ci_low),
mean_ci_high = mean(ci_high),
# rather than calculating the proportion of significant results by hand, use
# {simhelpers}'s dedicated function for the empirical rejection rate.
# calc_rejection() also returns the number of iterations and the Monte Carlo
# Standard Error of the rate; select() keeps and renames just the rate itself
calc_rejection(data = pick(everything()),
p_values = p,
alpha = 0.05) |>
select(proportion_significant_p = rej_rate))
results_summary |>
mutate_if(is.numeric, round, digits = 2)
```
### Statistical power for equivalence test
When using a Two One-Sided Test (TOST), with a Smallest Effect Size of Interest (SESOI) of difference-in-means = ±0.2, what is the statistical power across different sample sizes (n per condition = 25, 50, 75, 100)? What N is needed for 80% power to detect a true null effect size as equivalent to zero?
Notes:
- We're only interested in statistical power here, so set the population difference in means to zero. Why?
- Because reasons (see Lakens et al., 2018), use the 90% Confidence Interval instead of 95% CI.