Homework3
Quincy Hou
2024-03-12
# 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")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.
Filtering to country, recoding and cleaning variables of interest
france_data <- ess %>%
# Filter for France
filter(cntry == "FR") %>%
# Recoding Lrsscale and cleaning trstplt (will use later)
mutate(
ID = case_when(
lrscale %in% 0:3 ~ "Left", # Values 0 to 3 in 'lrscale' are categorized as "Left"
lrscale %in% 7:10 ~ "Right", # Values 7 to 10 in 'lrscale' are categorized as "Right"
lrscale %in% 4:6 ~ "Moderate", # Values 4 to 6 in 'lrscale' are categorized as "Moderate"
lrscale %in% c(77, 88, 99) ~ NA_character_ # Values 77, 88, 99 in 'lrscale' are set as NA
),
trstplt = ifelse(trstplt %in% c(77, 88, 99), NA, trstplt),
# 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
) %>%
filter(!is.na(ID) & !is.na(educ.ba))Beginning of Long form coding for Chi-Square
Calculating and interpreting critical value for independence
In our case for education, recoded as having attained at least a BA vs. not, and the recoded lrscale as the corresponding categorical values of Left, Right and Moderate.
## Not necessary, but always good to rename when doing operations so you can backtrack to clean dataset
tab_dat <- france_data
# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$ID, tab_dat$educ.ba)
# Calculate row-wise sums of counts in table
rsums <- rowSums(mytab)
# Calculate column-wise sums of counts in table
csums <- colSums(mytab)
# Get the total count of the table (N)
N <- sum(mytab)
# Generate table of expected proportions
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.1867628# Table of expected frequencies
ftab <- 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.876# Critical Value for Independence:
alpha <- 0.05 # significance level
c_val <- qchisq(alpha, df = 1, lower.tail = FALSE)
cat("Critical Value:", round(c_val, 3), "\n")## Critical Value: 3.841Interpreting the Critical Value: If the Chi-squared test statistic that you compute (from the chisq.test function) exceeds 3.841, then we reject the null hypothesis at the 0.05 significance level. This would mean that there is a statistically significant association between feeling close to a party and education level in your dataset.
Calculating and interpreting Pearson’s chi-squared statistic
# Pearson's chi-squared statistic:
test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")## Pearson's X^2: 148.3357comparison to the critical value: We previously calculated a critical value of 3.841. Our observed X-squared value (148.3357) is much larger than this critical value.
Interpretation: The difference between our observed data and what we’d expect under the assumption of independence is statistically significant. Thus, there is a statistically significant association between one’s political preference and one’s educational level. This means that the probability of being on one political side is not the same across different education levels.
Calculating the p-value
# p-value:
p_val <- pchisq(test_stat, df = 1, lower.tail = F)
cat("p-value:", round(p_val, 4), "\n")## p-value: 0# Using chisq.test for validation:
cat("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: This is another representation of the p-value in scientific notation because it’s extremely close to 0.
The chi-squared test tells us if two categorical variables are independent or if they are associated in some manner.
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.
The large chi-squared statistic of 148.3357 further reinforces this point. If the two variables were independent, we’d expect a chi-squared statistic closer to 0.
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).
Step 1: Calculate the test statistic of your sample
Important Note: to assign the appropriate test statistic, please refer to the summary table presented during the lecture.
Filtering to country, recoding and cleaning variables of interest
france_data <- ess %>%
# Filter for France
filter(cntry == "FR") %>%
# Recoding Lrsscale and cleaning trstplt (will use later)
mutate(
sbbsntx = case_when(
sbbsntx == 1 ~ "Agree Strongly", # Recode 1 to "Yes"
sbbsntx == 2 ~ "Agree", # Recode 2 to "No"
sbbsntx == 3 ~ "Neither Agree nor disagree",
sbbsntx == 4 ~ "Disagree",
sbbsntx == 5 ~ "Disagree Strongly",
sbbsntx %in% c(7, 8, 9) ~ NA_character_, # Handle other specific cases where you want to set it as NA
TRUE ~ NA_character_ # Keep other values as-is but ensure they are characters
),
# Recoding domicil
domicil = case_when(
domicil == 1 ~ "Big city",
domicil == 2 ~ "Suburbs and Outskirts",
domicil == 3 ~ "Town or Small City",
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))Using combination of unique() and is.na to make sure there are no NAs before proceeding. Even after you thought you cleaned, and table() does not return NAs, there can still be NAs that introduce issues. One way to double check is the the unique() function. With unique(), you can identify all the unique values for a given variable, including NAs. Then, we will add a line using is.na to make sure it’s clean.
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
## 36.01302Interpretation: The t-statistic measures the difference between the observed sample mean and the population mean (or another sample mean, in the case of two-sample tests) 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 36.01 is quite large, indicating a significant difference in sbbsntx and domicil.
Step 2: Simulate the null distribution
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.002Since the p_value has a value of under 0.01, under the assumptions of frequentist hypothesis testing, this means there is strong evidence to reject the null.
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.
Step 4: Visualize
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 8.12
## 2 2 13.2
## 3 3 17.7
## 4 4 22.1
## 5 5 12.0
## 6 6 11.6
## 7 7 15.7
## 8 8 10.5
## 9 9 20.5
## 10 10 9.85
## # ℹ 990 more rowsInterpretation
The observed redline is the Chi-Square statistic which is far to the right. This would suggests that the observed statistic is greater than what would be expected if the sbbsntx and domicil were independent. This therefore suggest that there is a siginificant association between the two variables.
The P-value is also significantly small, being less then 0.01, which means there is strong evidence against null hypothesis of independence. Therefore further proving a further association between the two variables.
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).
Step 1: Calculate the test statistic of your sample
Filtering to country, recoding and cleaning variables of interest
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))null_lrs <- 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") # replace in between quotation marks appropriate test statistic
print(null_lrs$stat)## X-squared
## 40.85117Interpretation: The t-statistic measures the difference between the observed sample mean and the population mean (or another sample mean, in the case of two-sample tests) 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 40.85 is quite large, indicating a significant difference in sbbsntx and domicil.
Step 2: Simulate the null distribution
null_dist <- 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: Calculate the p-value of your sample
p_val <- null_dist %>% # 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.016Since the p_value has a value of under 0.01, under the assumptions of frequentist hypothesis testing, this means there is strong evidence to reject the null.
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.
Step 4: Visualize
conf_int <- null_dist %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_dist %>%
visualize() +
shade_p_value(obs_stat = null_lrs, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int)## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
null_dist## Response: ccrdprs (factor)
## Explanatory: ID (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 11.4
## 2 2 17.4
## 3 3 16.5
## 4 4 11.5
## 5 5 13.3
## 6 6 22.9
## 7 7 14.4
## 8 8 8.60
## 9 9 11.9
## 10 10 26.9
## # ℹ 990 more rowsInterpretation
The observed redline is the Chi-Square statistic which is far to the right. This would suggests that the observed statistic is greater than what would be expected if the ccrdprs and lrscale were independent. This therefore suggest that there is a siginificant association between the two variables.
The confidence interval reoresents the 95% confidence interval of the null distribution.It indicates that the observed statistic is significatly different from what would have been expencted under the null hypothesis.
The P-value is also significantly small, being under 0.01, which means there is strong evidence against null hypothesis of independence. Therefore further proving a further association between the two variables.
Therefore there is a statistical association between the political orientation and the levels of personal responsibility people feel to reduce climate change. As observed through chi-square statistics which falls well outside of the range of what we would expect from no association.