In emeyers/SDS230: Tools for the class Data Exploration and Analysis

$\$

# 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")

$\$

Overview

• Review: Randomization test for two means using a t-statistic
• Parametric tests for two means
• Extra material: Hypothesis tests for a single mean
• Extra material: Connections between hypothesis tests and confidence intervals

$\$

Part 1: Permutation tests for comparing two means using a t-statistic

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.

$\$

Step 1: State the null and alternative hypotheses

$H_0: \mu_{gingko} - \mu_{placebo} = 0$

$H_A: \mu_{gingko} - \mu_{placebo} \ne 0$

$\alpha = 0.05$

$\$

Step 2a: Plot the data

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"))

$\$

Step 2b: Calculate the observed statistic using a t-statistic

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)

$\$

Step 3: Create the null distribution using a permutation/randomization test

# 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")

$\$

Step 4: Calculate a p-value

# 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)

$\$

Step 5: Make a decision

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?

Part 2: Parametric tests for comparing two means in R

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!

$\$

Step 1: State the null and alternative hypotheses

Same as before...

$H_0: \mu_{gingko} - \mu_{placebo} = 0$

$H_A: \mu_{gingko} - \mu_{placebo} > 0$

$\alpha = 0.05$

$\$

Step 2b: Calculate the observed statistic using a t-statistic

Same as before:

$$t = \frac{\bar{x}_t - \bar{x}_c}{\sqrt{\frac{s^2_t}{n_t} + \frac{s^2_c}{n_c}}}$$

$\$

Step 3: Create the null distribution

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

# get the degrees of freedom




# visualize the t-distribution density curve using the dt() function 







# how does this compare to our t-distribution created by shuffling? 

$\$

Step 4: Calculate a p-value

We can get $P(X < stat)$ for a t-distribution using the pt() function.

$\$

Step 5: Make a decision

How does our p-value and decision compare to the p-value decision we got from the permutation test?

$\$

Built in R functions for t-tests

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.

Why is the p-value slightly different than what we got when we used the pt() function?

$\$

$\$

Extra material: Hypothesis test for a single mean using a t-statistic

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!

Step 1: State the null and alternative hypotheses

$\$

Step 2a: Plot the data

# install.packages("ggplot2")
library(ggplot2)

msleep <- msleep  # make the data frame visible


# get the amount of sleep as a vector



# how much sleep humans get



# remove humans from the data



# visualize the data

$\$

Step 2: Calculate the observed statistic using a t-statistic

The formula for a t-statistic is:

$$t = \frac{\bar{x} - \mu_0}{\frac{s}{\sqrt{n}}}$$

# calculate the observed statistic 

$\$

Step 3: Create the null distribution using a t-statistic

A computational method we can use the create the null distribution consists of:

1. Transforming our sample of data to have a mean value equal to the parameter value specified in the null hypothesis.

2. 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).

3. 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



# use a for loop to create shuffled treatment and control groups and shuffled statistics 



       # resample the data


       # save the statistic of interest

$\$

Step 4: Calculate a p-value

# plot the null distribution again with a red line a the value of the observed statistic




# create a t-distribution 







# calculate the randomization method p-value




# calculate the t-statistic p-value





# using R's t-test function to get a p-value

$\$

Step 5: Make a decision

$\$

Extra material: Confidence intervals for a single mean

$\$

Creating a bootstrap distribution

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 

$\$

Plot the bootstrap distribution and a parametric approximate sampling distribution

# plot the bootstrap distribution





# create a normal distribution as a parametric approximation if using a noncentral t-distribution 





# 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")

$\$

Create confidence intervals

# getting bootstrap 95% confidence intervals using the percentile method




# getting bootstrap 95% confidence intervals using a t-distribution 




# get 95% confidence intervals using the t-distribution

$\$

Part 3: Comparing null and bootstrap distributions

Simulation based hypothesis test using the mean statistic

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



# create a data set consistent with the null distribution that the mean amount of sleep is 8 hours



# use a for loop to create shuffled treatment and control groups and shuffled statistics 



       # resample the data


       # calculate a bootstrap mean statistic




# plot a histogram of the null distribution 






# calculate the p-value

$\$

Compare null and bootstrap distribution

# plot the bootstrap distribution










# put green lines at the 95% confidence interval limits



# put a orange line at the observed statistic (x-bar) value



# put an purple line at the center of the null distribution (null parameter value)

emeyers/SDS230 documentation built on Jan. 18, 2024, 1:01 a.m.