In emeyers/SDS100: Tools for the class Introductory Statistics

knitr::opts_chunk$set(echo = TRUE)

library(SDS100)

download_data("popularkids.txt")
download_data("analgesics.txt")

library(SDS100)

$\$

Differences in categorical distributions:chi-square test for association

The data is from:

Chase, M. A., and Dummer, G. M. (1992), "The Role of Sports as a Social Determinant for Children," Research Quarterly for Exercise and Sport, 63, 418-424

The subjects, students in grades 4-6 in selected schools in Michigan, were asked the following question: What would you most like to do at school?

A. Make good grades. B. Be good at sports. C. Be popular.

Demographic information was also collected for each student.

# chi-square test for association
# grades, popularity, sports preference across grades 


# data code book
# https://math.tntech.edu/machida/ISR/3070/DASL/a/popularkids.html?dataset=3070/DASL/a/popularkids.txt


library(SDS100)

#download_data("popularkids.txt")


kids <- read.table("popularkids.txt", header = TRUE)

grade <- kids$Grade
goals <- kids$Goals


# Step 1: 
# H0: There is no difference between preference (of grades, popularity and sports) between grades
# HA: There is difference in preference between grades




# Step 2:

# observed table
(observed_counts <- table(goals, grade))


# visualize the data
par(mfrow = c(3, 1))
barplot(observed_counts[, 1], main = "4th grade")
barplot(observed_counts[, 2], main = "5th grade")
barplot(observed_counts[, 3], main = "6th grade")


# expected counts
row_sums <- rowSums(observed_counts)
col_sums <- colSums(observed_counts)
n_tot <- sum(observed_counts)

expected_counts <- outer(row_sums, col_sums)/n_tot

# can visualize the expected counts
par(mfrow = c(3, 1))
barplot(expected_counts[, 1], main = "4th grade")
barplot(expected_counts[, 2], main = "5th grade")
barplot(expected_counts[, 3], main = "6th grade")


chi_stat <- sum(((observed_counts - expected_counts)^2)/expected_counts)





# step 3: null dist
degree_free <- (nrow(observed_counts) - 1) * (ncol(observed_counts) - 1)

par(mfrow = c(1, 1))
x_vals <- seq(0, 20, length.out = 1000)
y_vals <- dchisq(x_vals, degree_free)
plot(x_vals, y_vals, type = "l")
abline(v = chi_stat, col = "red")





# step 4: p-value
pchisq(chi_stat, degree_free, lower.tail = FALSE)




# step 5: decision


# built in R function to do it!
chisq.test(observed_counts)





# Extra: can also look at rural, suburan and urban
# Is there a difference between these? 


living_location <- kids$Urban.Rural


# observed table
(observed_counts <- table(goals, living_location))


# visualize the data
par(mfrow = c(3, 1))
barplot(observed_counts[, 1], main = "Rural")
barplot(observed_counts[, 2], main = "Suburban")
barplot(observed_counts[, 3], main = "Urban")


chisq.test(observed_counts)

$\$

One-way ANOVA

A pharmaceutical company tested three formulations of a pain relief medicine for migraine headache sufferers. For the experiment, 27 volunteers were selected and 9 were randomly assigned to one of three drug formulations. The subjects were instructed to take the drug during their next migraine headache episode and to report their pain on a scale of 1 = no pain to 10 = extreme pain 30 minutes after taking the drug.

data from: https://dasl.datadescription.com/datafile/analgesics/?_sfm_methods=Analysis+of+Variance&_sfm_cases=4+59943

$\$

Step 1

$H_0$: all $\mu$ the same $H_A: at least one $\mu$ is different

# download_data("analgesics.txt")

drugs <- read.table("analgesics.txt", header = TRUE)




# step 2

drug_type <- drugs$Drug
pain_rating <- drugs$Pain

by(pain_rating, drug_type, mean)


# visualize the data
boxplot(pain_rating ~ drug_type)



# see if there is the same variability between groups
by(pain_rating, drug_type, sd)


# get the observed F-statistic using teh get_F_stat() function 
F_stat <- get_F_stat(pain_rating, drug_type)






# step 3: visualize null distribution

df1 <- 2    # K - 1
df2 <- 27 - 3   # N - K

x_vals <- seq(0, 15, length.out = 1000)
y_vals <- df(x_vals, df1, df2)
plot(x_vals, y_vals, type = "l")
abline(v = F_stat, col = "red")





# step 4

pf(F_stat, df1, df2, lower.tail = FALSE)




# step 5:




# built in R functions
anova_model <- aov(pain_rating ~ drug_type)

summary(anova_model)

$\$

Inference for regression - Bechdel data

Step 1

$H_0: \beta = 0$ $H_A: \beta > 0$

# load the library with the data
library(fivethirtyeight)

# remove missing values
bechdel <- na.omit(bechdel)

# extract variables of interest
budget <- bechdel$budget/10^6
revenue <- bechdel$domgross/10^6


lm_fit <- lm(revenue ~ budget)

obs_stat <- coef(lm_fit)[2]



# step 3: using randomization methods

null_dist <- do_it(10000) * {

  lm_shuff <- lm(revenue ~ shuffle(budget))
  coef(lm_shuff)[2]

}


# visualize the null distribution: where is the observed stat
hist(null_dist, breaks = 100)



# step 4 
pnull(obs_stat, null_dist, lower.tail = FALSE)




# step 5





# using parametric methods with R's built-in functions

summary(lm_fit)





# could do the bootstrap to get a CI for beta0 using SDS100::resample_pairs()

boot_dist <- do_it(1000) * {

  resampled_data <- resample_pairs(budget, revenue)
  lm_boot <- lm(vector2 ~ vector1, data = resampled_data)
  coef(lm_boot)[2]

}

quantile(boot_dist, c(.025, .975))





# we can also use R's built in functions to test if rho = 0...
cor.test(budget, revenue)