Homework 3
Melanie Fidel Rodriguez
In this tutorial, we will learn how to use the infer package for hypothesis testing as well as do tables and calculations the “long way” for the Chi-Square test statistic. Let’s get started.
Setting up environment
# List of packages
packages <- c("tidyverse", "infer", "fst") # add any you need here
# Install packages if they aren't installed already
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)
# Load the packages
lapply(packages, library, character.only = TRUE)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors## Warning: package 'infer' was built under R version 4.3.3## [[1]]
## [1] "lubridate" "forcats" "stringr" "dplyr" "purrr" "readr"
## [7] "tidyr" "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "infer" "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [7] "readr" "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [13] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "fst" "infer" "lubridate" "forcats" "stringr" "dplyr"
## [7] "purrr" "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [13] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [19] "base"ess <- read_fst("All-ESS-Data.fst")Task 1
Based on the recoding of lrscale for “left” and “right” and omitting “moderates” (see Tutorial 4), and educ.ba from this tutorial, do the long form coding for Chi-Square. Using the steps outlined in the tutorial, first generate the tables of expected proportions and frequencies. Determine and interpret the critical value for independence. Finally, determine and interpret both the Pearson’s chi-squared statistic and the p-value. Discuss main takeaways.
france_data <- ess %>%
filter(cntry == "FR") %>%
mutate(
lrscale = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
TRUE ~ NA_character_
),
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
mutate(
edulvla = ifelse(edulvla %in% c(55, 77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888, 9999), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba)
)tab_dat <- france_data
mytab <- table(tab_dat$lrscale, tab_dat$educ.ba)
rsums <- rowSums(mytab)
csums <- colSums(mytab)
N <- sum(mytab)
ptab <- tcrossprod(rsums/N, csums/N)
cat("Table of Expected Proportions:\n")## Table of Expected Proportions:print(ptab)## [,1] [,2]
## [1,] 0.1311243 0.3914592
## [2,] 0.1197912 0.3576253ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 1181.561 3527.439
## [2,] 1079.439 3222.561alpha <- 0.05
c_val <- qchisq(alpha, df = 1, lower.tail = FALSE)
cat("Critical Value:", round(c_val, 3), "\n")## Critical Value: 3.841test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")## Pearson's X^2: 100.8551p_val <- pchisq(test_stat, df = 1, lower.tail = F)
cat("p-value:", round(p_val, 4), "\n")## p-value: 0cat("Chi-squared test result:\n")## Chi-squared test result:print(chisq.test(mytab, correct = F))##
## Pearson's Chi-squared test
##
## data: mytab
## X-squared = 100.86, df = 1, p-value < 2.2e-16The p-value of 2.2e-16 is essentially zero, indicating that the two variables are independent from each other, and that the null is rejected. The chi-square reinforces the deviation from independence.
Task 2
Do Steps 1 to 3 using the infer package (as we did in the tutorial), for the variables sbbsntx (recoded to the corresponding 5 categories) as the response and domicil (recoded to urban and rural, setting 2 to urban instead of peri-urban) as the explanatory variable for the country of France.
See variable info here: https://ess.sikt.no/en/datafile/ffc43f48-e15a-4a1c-8813-47eda377c355/92?tab=1&elements=[%22874a86e2-a0f6-40b4-aef7-51c8eed98a7d/1%22])
Provide interpretations of the output (consider the variable info).
france_data <- ess %>%
filter(cntry == "FR") %>%
mutate(
geo = recode(as.character(domicil),
'1' = "Urban",
'2' = "Urban",
'3' = "Rural",
'4' = "Rural",
'5' = "Rural",
'7' = NA_character_,
'8' = NA_character_,
'9' = NA_character_
),
sbbsntx = case_when(
sbbsntx == 1 ~ "Agree Strongly",
sbbsntx == 2 ~ "Agree",
sbbsntx == 3 ~ "Neither agree nor disagree",
sbbsntx == 4 ~ "Disagree",
sbbsntx == 5 ~ "Disagree strongly",
sbbsntx %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(sbbsntx)
),
) %>%
filter(!is.na(sbbsntx), !is.na(geo))test_stat <- france_data %>%
specify(explanatory = geo, # change variable name for explanatory variable
response = sbbsntx) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "chisq")
print(test_stat$stat)## X-squared
## 19.97187null_distribution <- france_data %>%
specify(explanatory = geo,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "chisq")p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "two-sided")
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.002Since there is a p-value of 0, there is strong evidence to reject the null. This suggests a significant association between the two categorical variables being tested.
Task 3
Now do Step 4 (i.e., the null distribution visualization with confidence intervals) for the same variables and country as Task 2, and interpret the output.
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided")## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
conf_int <- null_distribution %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int)## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
null_distribution## Response: sbbsntx (factor)
## Explanatory: geo (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 7.49
## 2 2 11.9
## 3 3 10.8
## 4 4 3.55
## 5 5 2.72
## 6 6 0.227
## 7 7 2.77
## 8 8 3.04
## 9 9 2.16
## 10 10 5.69
## # ℹ 990 more rowsThere is evidence to reject the null, as is evident when visualizing the null distribution next to the observed test statistic, or rather the red line shown in the graph.
Task 4
Conduct Steps 1 to 3 using the infer package (as we did in the tutorial), for the variables ccrdprs (leaving it on the 0-10 numeric scale) as the response (or outcome) variable and lrscale (recoded as left and right, omitting “moderates” as we did in Task 1) as the explanatory variable.
Variable info here: https://ess.sikt.no/en/datafile/ffc43f48-e15a-4a1c-8813-47eda377c355/92?tab=1&elements=[%2283b08cf2-508e-49a0-9fc8-c3cc281290ec/4%22]
Provide interpretations of the output (consider the variable info).
table(france_data$lrscale)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88
## 197 110 287 417 365 1073 358 405 348 95 133 90 161table(france_data$ccrdprs)##
## 0 1 2 3 4 5 6 7 8 9 10 66 88 99
## 47 17 46 43 71 261 204 374 464 164 282 23 10 1france_data <- ess %>%
filter(cntry == "FR") %>%
mutate(
lrscale = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
TRUE ~ NA_character_,
ccrdprs %in% c(66, 77, 88, 99) ~ NA_character_,
)
) %>%
filter(!is.na(lrscale), !is.na(ccrdprs))test_stat <- france_data %>%
specify(explanatory = lrscale, # change variable name for explanatory variable
response = ccrdprs) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "t") ## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "Left" - "Right", or divided in the order "Left" / "Right" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Left", "Right")` to the calculate() function.print(test_stat$stat)## t
## 0.5130018null_distribution <- france_data %>%
specify(explanatory = lrscale,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "t")## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "Left" - "Right", or divided in the order "Left" / "Right" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Left", "Right")` to the calculate() function.p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "two-sided")p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.574A p-value of 0.6 suggests that there is not enough evidence to reject the null hypothesis, indicating that there is not a statistically significant difference in the mean ccrdprs scores between the “left” and “right” categories of lrscale.
Task 5
Now do Step 4 (i.e., the null distribution visualization with confidence intervals) for the same variables and country as Task 4, and interpret the output.
conf_int <- null_distribution %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int)
null_distribution## Response: ccrdprs (numeric)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.0266
## 2 2 0.113
## 3 3 -0.808
## 4 4 -0.484
## 5 5 1.51
## 6 6 0.0755
## 7 7 -0.947
## 8 8 -0.781
## 9 9 -1.26
## 10 10 -0.681
## # ℹ 990 more rowsThere is evidence to not reject the null, as is evident when visualizing the null distribution next to the observed test statistics, or rather the red line shown in the graph.