Homework 3 (222)
Christian Dizon
March 12
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.4
## ✔ 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")Filtering to country, recoding and cleaning variables of interest
france_data <- ess %>%
# Filter for France
filter(cntry == "FR") %>%
# Recoding clsprty and cleaning trstplt (will use later)
mutate(
clsprty = case_when(
clsprty == 1 ~ "Yes", # Recode 1 to "Yes"
clsprty == 2 ~ "No", # Recode 2 to "No"
clsprty %in% c(7, 8, 9) ~ NA_character_, # Handle other specific cases where you want to set it as NA
TRUE ~ as.character(clsprty) # Keep other values as-is but ensure they are characters
),
trstplt = ifelse(trstplt %in% c(77, 88, 99), NA, trstplt),
# Recoding gender
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
# Recoding education with case_when for consistency
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
# Ensure all NAs are handled uniformly
mutate(
edulvla = ifelse(edulvla %in% c(77, 88), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba) # Force educ.ba to be treated as character
)Homework 3 (5%): due by next lecture on March 12th
Instructions: Start a new R markdown for the homework and call it “Yourlastname_Firstname_Homework_2”.
Copy everything below from Task 1 to Task 5. Keep the task prompt and questions, and provide your code and answer underneath.
Remember: you need all the steps for your code to work, including loading your data – otherwise it will not knit.
To generate a new code box, click on the +C sign above. Underneath your code, provide your answer to the task question.
When you are done, click on “Knit” above, then “Knit to Html”. Wait for everything to compile. If you get an error like “Execution halted”, it means there are issues with your code you must fix. When all issues are fixed, it will prompt a new window. Then click on “Publish” in the top right, and then Rpubs (the first option) and follow the instructions to create your Rpubs account and get your Rpubs link for your document (i.e., html link as I provide for the tutorial).
Note: Make sure to provide both your markdown file and R pubs link. If you do not submit both, you will be penalized 2 pts. out of the 5 pts. total.
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(
ID = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
lrscale %in% 4:6 ~ "Moderate",
lrscale %in% c(77, 88, 99) ~ NA_character_
),
trstplt = ifelse(trstplt %in% c(77, 88, 99), NA, trstplt),
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(77, 88), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba)
) %>%
filter(!is.na(ID) & !is.na(educ.ba))tab_dat <- france_data
mytab <- table(tab_dat$ID, 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.06135806 0.2044319
## [2,] 0.11343879 0.3779537
## [3,] 0.05605487 0.1867628ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 1087.0808 3621.919
## [2,] 2009.7951 6696.205
## [3,] 993.1241 3308.876alpha <- 0.05
c_val <- qchisq(alpha, df = 1, lower.tail = FALSE)
cat("Critical Value:", round(c_val, 3), "\n")## Critical Value: 3.841Interpreting the Critical Value: As a result of the critical value exceeding or exactly being 3.841, we reject the null hypothesis at the 0.05 significance level. This means that there is a statistically significant assocation of individuals feeling close to a party and education level.
test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")## Pearson's X^2: 148.3357Comparison to Critical Value: Since the critical value calculated to 3.841, the observed X-squared value by Pearson is significantly larger than the critical value that was calculated.
Interpretation: The difference between the observed data and the assumption of independence was statistically significant. There is a significant assocaiton between feeling close to a party and to one’s educational level. This means that the probability being close or attached to a specific party is not the same across different educational level.
p_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 = 148.34, df = 2, p-value < 2.2e-16p-value < 2.2e-16: The p-value representation in a scientific notation is extremely close to 0.
The chi-squared test tells us if two categorical variables are independent or if they are associated in some manner. In our test: The extremely small p-value (essentially 0) means that we reject the null hypothesis that the two variables are independent of each other. In other words, there’s strong evidence to suggest that there’s a significant association between the two variables in mytab.
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 <- france_data %>%
filter(cntry == "FR") %>%
mutate(
sbbsntx = case_when(
sbbsntx == 1 ~ "agree strongly",
sbbsntx == 2 ~ "agree",
sbbsntx == 3 ~ "neither agree nor disgree",
sbbsntx == 4 ~ "disagree",
sbbsntx == 5 ~ "disgree strongly",
sbbsntx %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ NA_character_
),
domicil = case_when(
domicil == 1 ~ "Big city",
domicil == 2 ~ "Suburbs",
domicil == 3 ~ "Town",
domicil == 4 ~ "Village",
domicil == 5 ~ "Country Side",
domicil %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ NA_character_
)
)%>%
filter(!is.na(sbbsntx) & !is.na(domicil))Step 1:
test_stat <- france_data %>%
specify( response = sbbsntx, # change variable name for explanatory variable
explanatory = domicil) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "chisq") # replace in between quotation marks appropriate test statistic
print(test_stat$stat)## X-squared
## 31.81446Interpretation: The t-statistic measures the difference between the observed sample mean and the population mean in units of standard error. The larger the magnitude of t (regardless of the sign), the greater the evidence against the null hypothesis. A t-statistic of 28.1 is large enough to indicate a significant difference between the levels of climate change betwen urban and rural areas.
Step 2
null_distribution <- france_data %>%
specify(explanatory = sbbsntx,
response = domicil) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>% # only line we are adding, use reps = 1000 as standard but change if too computationally demanding
calculate(stat = "chisq")Step 3: Calculate the p-value of your sample
p_val <- null_distribution %>% # replace name here if you assigned something other than null_distribution above
get_pvalue(obs_stat = test_stat, direction = "two-sided") # would only replace test_stat if assigned another name in Step 1
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.02The p-value of the taken sample becomes 0.078.
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.
null_distribution## Response: domicil (factor)
## Explanatory: sbbsntx (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 26.2
## 2 2 19.9
## 3 3 9.92
## 4 4 13.8
## 5 5 10.0
## 6 6 20.0
## 7 7 16.5
## 8 8 7.46
## 9 9 20.4
## 10 10 14.0
## # ℹ 990 more rowsThe visualization of the null distribution alongside the observed test statistic (the red line) indicates assumption under the frequentist hypothesis testing, to reject the null. If the observed statistic falls within the distribution, it will be shaded in red, for those areas are observed to have extreme stats.
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).
france_data <- ess %>%
# Filter for France
filter(cntry == "FR") %>%
# Recoding Lrsscale and cleaning trstplt (will use later)
mutate(
ID = case_when(
lrscale == 0 ~ "Left",
lrscale == 1 ~ "Left",
lrscale == 2 ~ "Left",
lrscale == 3 ~ "Left",
lrscale == 4 ~ "Moderate",
lrscale == 5 ~ "Moderate",
lrscale == 6 ~ "Moderate",
lrscale == 7 ~ "Right",
lrscale == 8 ~ "Right",
lrscale == 9 ~ "Right",
lrscale == 10 ~ "Right",
lrscale %in% c(77, 88, 99) ~ NA_character_ # Values 77, 88, 99 in 'lrscale' are set as NA
),
# Recoding domicil
ccrdprs = case_when(
ccrdprs == 1 ~ "1",
ccrdprs == 2 ~ "2",
ccrdprs == 3 ~ "3",
ccrdprs == 4 ~ "4",
ccrdprs == 5 ~ "5",
ccrdprs == 6 ~ "6",
ccrdprs == 7 ~ "7",
ccrdprs == 8 ~ "8",
ccrdprs == 9 ~ "9",
ccrdprs == 10 ~ "10",
TRUE ~ NA_character_
)
)%>%
filter(!is.na(ID) & !is.na(ccrdprs))#Step 1
test_stat <- france_data %>%
specify(explanatory = ID, # change variable name for explanatory variable
response = ccrdprs) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "chisq")
print(test_stat$stat)## X-squared
## 40.85117Interpretation: The t-statistic measures the difference between the observed sample mean and the population mean in units of standard error. The larger the magnitude of t (regardless of the sign), the greater the evidence against the null hypothesis. The range between 40.85 to a t-statistic of 41.1 is large enough to indicate a significant difference between the political levels between left and right, omitting “moderates” supporters.
#Step 2
null_distribution <- france_data %>%
specify(explanatory = ID,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>% # only line we are adding, use reps = 1000 as standard but change if too computationally demanding
calculate(stat = "chisq")#Step 3
p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "two-sided")
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.006The p-value for the political level of “left”, “right”, and “moderate” is 0.
Task 5
#Step 4
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.
## The visualization of the null distribution alongside the observed
test statistic (the red line) indicates assumption under the frequentist
hypothesis testing, to reject the null. As observed, there is no red
shaded areas to indicate extreme levels of either political “left”,
“right”, or “moderates”.
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.