Johns_Ruth_Homework_3
Ruth Johns
2024-03-11
# 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## [[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")france_data <- read.fst("france_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 for France
filter(cntry == "FR") %>%
mutate(
clsprty = case_when(
clsprty == 1 ~ "Yes",
clsprty == 2 ~ "No",
clsprty %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(clsprty)
),
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)
)unique(france_data$gndr)## [1] "Female" "Male"unique(france_data$clsprty)## [1] "Yes" "No" NA# make sure it's really cleaned
france_data <- france_data %>% filter(!is.na(clsprty))
france_data <- france_data %>% filter(!is.na(trstplt))
# double-checking + quick distribution
unique(france_data$clsprty)## [1] "Yes" "No"table(france_data$clsprty)##
## No Yes
## 9098 9480tab_dat <- france_data
mytab <- table(tab_dat$clsprty, 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.1099746 0.3797444
## [2,] 0.1145921 0.3956889# Table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 2043.108 7054.892
## [2,] 2128.892 7351.108# 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.841# 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: 160.4057The x-squared value of 160.4057 greatly exceeds that of the calculated critical value of 3.841. What this means is that it is showing a statistical significance between the two. There is a statistically significant relationship between feeling close to a party and one’s educational level. The probability of feeling close to a party is not the same and can vary across different education levels.
# 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 = 160.41, df = 1, p-value < 2.2e-16The p-value is extremely small and very close to 0. This small value means that we reject the null hypothesis that the two variables are independent of each other. This shows that there’s a strong association between the variables. The larger chi-squared statistic also proves this further because if it were closer to 0, that would mean the variables would not be dependent on each other.
##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 for the country France
filter(cntry == "FR") %>%
#Recoding to the corresponding 5 categories
mutate(
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)
),
# Recoding to urban and rural
domicil = case_when(
domicil == 1 ~ "Urban",
domicil == 2 ~ "Urban",
domicil == 3 ~ "Rural",
domicil == 4 ~ "Rural",
domicil == 5 ~ "Rural",
domicil %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(domicil)
)
)test_stat <- france_data %>%
specify(explanatory = domicil, # change variable name for explanatory variable
response = sbbsntx) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "Chisq") # replace in between quotation marks appropriate test statistic## Warning: Removed 14999 rows containing missing values.print(test_stat$stat)## X-squared
## 19.97187The chi-squared statistic of 19.97187 shows a clear correlation between people residing in urban or rural areas and the range of agreement in relation to how much social benefits/services cost businesses too much in taxes/charges. This implies that the degree of agreement varies considerably amongst people depending on where they live.
null_distribution <- france_data %>%
specify(explanatory = domicil,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")## Warning: Removed 14999 rows containing missing values.p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "greater") ## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation
## based on the number of `reps` chosen in the `generate()` step.
## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information.p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0The p-value is 0 so that means that we reject the null hypothesis that the two variables are independent of each other and that there is a correlation.
##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.
conf_int <- null_distribution %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "greater") +
shade_confidence_interval(endpoints = conf_int)## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
null_distribution## Response: sbbsntx (factor)
## Explanatory: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 5.46
## 2 2 8.24
## 3 3 3.31
## 4 4 6.62
## 5 5 3.65
## 6 6 8.29
## 7 7 1.20
## 8 8 0.933
## 9 9 1.55
## 10 10 12.5
## # ℹ 990 more rowsThe distribution in the graph produced is right-skewed, meaning that its mean, mode, and median are all less than 10. The confidence intervals fall between 0 and a little bit above 10 using a 95% confidence interval meaning the true value falls between this range. And since this produces a higher certainty,the red line proves against the null hypothesis since it sits at 20 and does not fall within the null distribution.
##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 'lrscale'
mutate(
lrscale = 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
),
# Set values 66, 77, 88, 99 in 'ccrdprs' to NA
ccrdprs = ifelse(ccrdprs %in% c(66, 77, 88, 99), 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 = "F") ## Warning: Removed 15413 rows containing missing values.print(test_stat$stat)## [1] 15.41871The ANOVA value of 15.41871 indicates that the groups compared differ from one another. There is a correlation between a person’s political stance (left or right) and their sense of personal accountability for addressing climate change.
null_distribution <- france_data %>%
specify(explanatory = lrscale,
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 = "F")## Warning: Removed 15413 rows containing missing values.p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "greater") ## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation
## based on the number of `reps` chosen in the `generate()` step.
## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information.p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0The p-value is 0 so that means that we reject the null hypothesis.
##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 = "greater") +
shade_confidence_interval(endpoints = conf_int)## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
null_distribution## Response: ccrdprs (numeric)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.616
## 2 2 0.299
## 3 3 0.219
## 4 4 0.935
## 5 5 1.96
## 6 6 0.340
## 7 7 0.912
## 8 8 0.859
## 9 9 2.85
## 10 10 0.325
## # ℹ 990 more rowsThe distribution in the graph produced is right-skewed, meaning that its mean, mode, and median are all less than 5. The confidence intervals fall between 0 and almost 5 using a 95% confidence interval meaning the true value falls between this range. And since this produces a higher certainty,the red line proves against the null hypothesis since it sits above 15 and does not fall within the null distribution.