$\$
# install.packages("latex2exp") library(latex2exp) options(scipen=999) knitr::opts_chunk$set(echo = TRUE) set.seed(123)
# get some images and data that are used in this document SDS230::download_image("student_t.png") SDS230::download_data("gingko_RCT.rda") SDS230::download_data("alcohol.rda")
$\$
$\$
Below is the code for running a randomization test for the experiment by Solomon et al (2002) to see if there is evidence that taking a gingko pills affects cognition.
We have the code here so that we can compare the results of running this randomization test to the results from running a parametric t-test.
$\$
$H_0: \mu_{gingko} - \mu_{placebo} = 0$ $H_A: \mu_{gingko} - \mu_{placebo} \ne 0$
$\alpha = 0.05$
$\$
load("gingko_RCT.rda") # plot the data boxplot(gingko, placebo, names = c("Gingko", "Placebo"), ylab = "Memory score") # create a stripchart data_list <- list(gingko, placebo) stripchart(data_list, group.names = c("Gingko", "Placebo"), method = "jitter", xlab = "Memory score", col = c("red", "blue"))
$\$
The formula for a t-statistic is:
$$t = \frac{\bar{x}_t - \bar{x}_c}{\sqrt{\frac{s^2_t}{n_t} + \frac{s^2_c}{n_c}}}$$
numerator <- mean(gingko) - mean(placebo) demoninator <- sqrt( var(gingko)/length(gingko) + var(placebo)/length(placebo) ) (obs_stat <- numerator/demoninator)
$\$
# combine the data from the treatment and control groups together combined_data <- c(gingko, placebo) n_gingko <- length(gingko) total <- length(combined_data) # use a for loop to create shuffled treatment and control groups and shuffled statistics null_distribution <- NULL for (i in 1:10000) { # shuffle data shuff_data <- sample(combined_data) # create fake treatment and control groups shuff_gingko <- shuff_data[1:n_gingko] shuff_placebo <- shuff_data[(n_gingko + 1):total] numerator_shuff <- mean(shuff_gingko) - mean(shuff_placebo) demoninator_shuff <- sqrt(var(shuff_gingko)/length(shuff_gingko) + var(shuff_placebo)/length(shuff_placebo)) # save the statistic of interest null_distribution[i] <- numerator_shuff/demoninator_shuff } # plot the null distribution as a histogram hist(null_distribution, breaks = 100, main = "Null distribution", xlab = "t-statistic")
$\$
# plot the null distribution again with a red line a the value of the observed statistic hist(null_distribution, breaks = 100, main = "Null distribution", xlab = "t-statistic") abline(v = obs_stat, col = "red") abline(v = abs(obs_stat), col = "red") # calculate the p-value (p_value_left_tail <- sum(null_distribution <= obs_stat)/10000) (p_value_right_tail <- sum(null_distribution >= abs(obs_stat))/10000) (p_value <- p_value_left_tail + p_value_right_tail)
$\$
Since r p_value
is greater than $\alpha = 0.05$ we can not reject the null hypothesis. Thus if we are using the Neyman-Pearson paradigm, we do not have sufficient evidence to say that the pill is effective.
When we used a statistic of $\bar{x}_t - \bar{x}_c$ in our randomization test in class 7 we got a p-value of 0.127. How do these results compare?
$\$
Let's redo our analysis of Solomon et al (2002) using a parametric probability distribution, which in this case is the t-distribution. The same 5 steps of hypothesis testing apply here as well!
$\$
Same as before...
$H_0: \mu_{gingko} - \mu_{placebo} = 0$
$H_A: \mu_{gingko} - \mu_{placebo} > 0$
$\alpha = 0.05$
$\$
Same as before:
$$t = \frac{\bar{x}_t - \bar{x}_c}{\sqrt{\frac{s^2_t}{n_t} + \frac{s^2_c}{n_c}}}$$
numerator <- mean(gingko) - mean(placebo) demoninator <- sqrt( var(gingko)/length(gingko) + var(placebo)/length(placebo) ) (obs_stat <- numerator/demoninator)
$\$
We will now use a parametric t-distribution (i.e., density function) as a null distribution. The t-distribution has one parameter called "degrees of freedom". We will set this parameter as the minimum of $n_t - 1$ or $n_c - 1$.
What are the degrees of freedom for this study?
Let's visualize the t-distribution
# visualize the t-distribution density curve using the dt() function (degree_free <- min(length(gingko), length(placebo)) - 1) x_vals <- seq(-5, 5, by = .01) y_vals <- dt(x_vals, degree_free) plot(x_vals, y_vals, type = "l", col = "blue") # how does this compare to our t-distribution created by shuffling? hist(null_distribution, freq = FALSE, breaks = 50, main = "Null distribution", xlab = "t-statistic") points(x_vals, y_vals, type = "l", col = "blue")
$\$
We can get $Pr(X < stat)$ for a t-distribution using the pt()
function.
(p_value_left_tail <- pt(obs_stat, degree_free, lower.tail = TRUE)) (p_value_right_tail <- pt(abs(obs_stat), degree_free, lower.tail = FALSE)) (p_value <- p_value_left_tail + p_value_right_tail)
$\$
How does our p-value and decision compare to the p-value decision we got from the permutation test?
$\$
We can use the built in t.test()
function to run a t-test as well.
Note: If you want to run one-tailed tests you can use the extra argument alternative
argument.
t.test(gingko, placebo)
Why is the p-value slightly different than what we got when we used the pt()
function?
$\$
$\$
Humans sleep for 8 hours a day (at least according to a data set that is included with ggplot package). Is this different from mean amount of sleep other mammals get? Let's investigate using permutation tests and t-tests!
$H_0: \mu = 8$ $H_A: \mu_{treat} \ne 8$
$\alpha = 0.05$
$\$
# install.packages("ggplot2") library(ggplot2) msleep <- msleep # make the data frame visible amount_of_sleep <- msleep$sleep_total # how much sleep humans get (human_sleep <- amount_of_sleep[msleep$name == "Human"]) # remove humans from the data amount_of_sleep <- amount_of_sleep[msleep$name != "Human"] boxplot(amount_of_sleep) abline(h = human_sleep, col = "red")
$\$
The formula for a t-statistic is:
$$t = \frac{\bar{x} - \mu_0}{\frac{s}{\sqrt{n}}}$$
numerator <- mean(amount_of_sleep) - 8 demoninator <- sd(amount_of_sleep)/sqrt(length(amount_of_sleep)) # the SE for the mean (obs_stat <- numerator/demoninator)
$\$
A computational method we can use the create the null distribution consists of:
Transforming our sample of data to have a mean value equal to the parameter value specified in the null hypothesis.
Resampling with replacement to get estimate the amount of variability we expect to see in the null distribution (and recomputing out statistic from this resampled data).
Repeating step 2 10,000 times to get a full null distribution.
# create a data set consistent with the null distribution that the mean amount of sleep is 8 hours null_sleep_data <- amount_of_sleep - mean(amount_of_sleep) + 8 # use a for loop to create shuffled treatment and control groups and shuffled statistics null_distribution <- NULL for (i in 1:10000) { # resample the data resample_data <- sample(null_sleep_data, replace = TRUE) # save the statistic of interest numerator_resample <- mean(resample_data) - 8 demoninator_resample <- sqrt( var(resample_data)/length(resample_data) ) null_distribution[i] <- numerator_resample /demoninator_resample }
$\$
# plot the null distribution again with a red line a the value of the observed statistic hist(null_distribution, breaks = 100, freq = FALSE, main = "Null distribution", xlab = "t-statistic", xlim = c(-6, 6)) # create a t-distribution x_vals <- seq(-6, 6, by = .01) y_vals <- dt(x_vals, df = length(amount_of_sleep) - 1) points(x_vals, y_vals, type = "l", col = "blue") abline(v = obs_stat, col = "red") # calculate the randomization method p-value (p_value_resample <- (sum(null_distribution >= obs_stat) + sum(null_distribution <= -obs_stat))/length(null_distribution)) # calculate the t-statistic p-value (p_value_t <- 2 * pt(-obs_stat, df = (length(amount_of_sleep) - 1))) # using R's t-test function to get a p-value t.test(amount_of_sleep, mu = 8)
$\$
Since the p-value is small we have strong evidence that humans get less sleep than other mammals. zzzzz
$\$
$\$
For confidence intervals we want to keep data on the original scale; i.e., we want to create a confidence interval for a range of plausible values for the actual mean length of sleep time for mammals.
Thus in our bootstrap (resampling) procedure we want to use the mean statistic applied to our resampled data (rather than calculating t-statistics from the resampled data, since creating a confidence interval for a t-statistics not that informative about how much mammals speed on average).
# use a for loop to create shuffled treatment and control groups and shuffled statistics boot_distribution <- NULL for (i in 1:10000) { # resample the data resample_data <- sample(amount_of_sleep, replace = TRUE) boot_distribution[i] <- mean(resample_data) }
$\$
# plot the bootstrap distribution hist(boot_distribution, breaks = 100, freq = FALSE, main = "Bootstrap distribution", xlab = "x-bar") # create a normal distribution as a parametric approximation if using a noncentral t-distribution x_vals <- seq(8, 13, by = .01) SE <- sd(amount_of_sleep)/sqrt(length(amount_of_sleep)) # SE = s/sqrt(n) y_vals <- dnorm(x_vals, mean(amount_of_sleep), SE) points(x_vals, y_vals, type = "l", col = "purple") # can create a non-central t-distribution as well by shifting and rescaling the central t-distribution # y_vals <- dt(x_vals, df = length(amount_of_sleep) - 1) # points( (SE * x_vals) + mean(amount_of_sleep), y_vals/SE, type = "l", col = "blue")
$\$
# getting bootstrap 95% confidence intervals using the percentile method (CI_boot_percentile <- quantile(boot_distribution, c(.025, .975))) # getting bootstrap 95% confidence intervals using a t-distribution t_critical_value <- qt(.975, df = length(amount_of_sleep) - 1) # the t* quantile value (CI_boot_t <- mean(amount_of_sleep) + c(-t_critical_value, t_critical_value) * sd(boot_distribution)) # get 95% confidence intervals using the t-distribution (CI_t <- mean(amount_of_sleep) + c(-t_critical_value, t_critical_value) * SE)
$\$
Let's rerun our simulation based hypothesis test using the mean statistic (instead of a t-statistic which was done above).
# step 2: observed statistic (obs_stat <- mean(amount_of_sleep)) # create a data set consistent with the null distribution that the mean amount of sleep is 8 hours null_sleep_data <- amount_of_sleep - mean(amount_of_sleep) + 8 # use a for loop to create shuffled treatment and control groups and shuffled statistics null_distribution <- NULL for (i in 1:10000) { # resample the data resample_data <- sample(null_sleep_data, replace = TRUE) # calculate a bootstrap mean statistic null_distribution[i] <- mean(resample_data) } hist(null_distribution, breaks = 100, main = "Null distribution", xlab = "x-bar", xlim = c(5, 11)) abline(v = obs_stat, col = "red") # still a p-value of 0 (pval_mean_stat <- (sum(null_distribution <= -obs_stat) + sum(null_distribution >= obs_stat))/length(null_distribution))
$\$
# plot the bootstrap distribution boot_hist <- hist(boot_distribution, breaks = 100) null_hist <- hist(null_distribution, breaks = 100) plot(boot_hist, col = "red", xlim = c(5, 13), main = "Null and Bootstrap distributions", xlab = "x-bar") plot(null_hist, col = "blue", add = TRUE) # put green lines at the 95% confidence interval limits abline(v = quantile(boot_distribution, c(.025, .975)), col = "green") # put a orange line at the observed statistic (x-bar) value abline(v = obs_stat, col = "orange", lwd = 3) # put an purple line at the center of the null distribution (null parameter value) abline(v = 8, col = "purple", lwd = 3)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.