Kashvi Shah | 1009413706
# List of packages
packages <- c("tidyverse", "fst", "modelsummary") # 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.1
## ✔ 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] "fst" "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [7] "readr" "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [13] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "modelsummary" "fst" "lubridate" "forcats" "stringr"
## [6] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [11] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [16] "utils" "datasets" "methods" "base"rm(list=ls()); gc() ## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 1208570 64.6 2196262 117.3 NA 2196262 117.3
## Vcells 2061502 15.8 8388608 64.0 16384 2885288 22.1ess <- read_fst("All-ESS-Data.fst")unique(ess$cntry)## [1] "AT" "BE" "CH" "CZ" "DE" "DK" "ES" "FI" "FR" "GB" "GR" "HU" "IE" "IL" "IT"
## [16] "LU" "NL" "NO" "PL" "PT" "SE" "SI" "EE" "IS" "SK" "TR" "UA" "BG" "CY" "RU"
## [31] "HR" "LV" "RO" "LT" "AL" "XK" "ME" "RS" "MK"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") %>%
# 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
)## 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$clsprty, 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.1101248 0.3815066
## [2,] 0.1138739 0.3944946# Table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 2059.444 7134.556
## [2,] 2129.556 7377.444# 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.841At the 0.05 significance level, we would reject the null hypothesis if the calculated Chi-squared test statistic exceeded 3.841, as determined by the chisq.test function. Therefore, this would suggest that, within your dataset, there is a statistically significant correlation between feeling near to an event and education level.
# 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: 166.1849In comparison, we compute a critical value of 3.841 and discover that our observed Chi-squared value (166.1849) significantly above this barrier.
Statistical significance is demonstrated by the significant difference between our observed data and the expected results based on the independence assumption. As such, there is a statistically significant correlation between educational attainment and the sense of being near to a party. This suggests that there are differences in the chances of feeling connected to a party between educational 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 = 166.18, df = 1, p-value < 2.2e-16The p-value that corresponds to the Chi-squared test is < 2.2e-16. Strong evidence is present against the null hypothesis, as indicated by the extremely small p-value, which is almost zero. Stated differently, the variables within the dataset have a statistically significant correlation with one another.
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.
france_data <- ess %>%
# Filter for 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_, # Handle other specific cases where you want to set it as NA
TRUE ~ as.character(sbbsntx) # Keep other values as-is but ensure they are characters
),
# 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_, # Handle other specific cases where you want to set it as NA
TRUE ~ as.character(domicil) # Keep other values as-is but ensure they are characters
)
)library(infer)
test_stat <- france_data %>%
specify(explanatory = domicil,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
calculate(stat = "chisq")## Warning: Removed 14999 rows containing missing values.print(test_stat$stat)## X-squared
## 19.97187The 19.97187 chi-squared statistic shows a significant relationship between living in an urban or rural location and the degree of agreement regarding the perceived cost of social benefits/services on enterprises in the form of taxes or charges. This implies that the degree of agreement varies considerably amongst people depending on where they live. Moreover, this discovery offers strong proof to reject the null hypothesis.
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 %>% # replace name here if you assigned something other than null_distribution above
get_pvalue(obs_stat = test_stat, direction = "greater") # would only replace test_stat if assigned another name in Step 1## 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 of 0 indicates substantial evidence to reject the null hypothesis, according to frequentist hypothesis testing principles. There is insufficient evidence to refute the null hypothesis, which holds that there is no correlation between place of residence (rural or urban) and agreement with the claim that “social benefits/services cost businesses too much in taxes/charges.” This high p-value indicates that, in the event that there was no true underlying link, the observed variations in attitudes between urban and rural individuals might have easily happened by coincidence. Based on this study, we do not have sufficient evidence to reject the null hypothesis, which suggests that the difference in agreement with the statement between urban and rural areas may not be as statistically significant as first believed.
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: sbbsntx (factor)
## Explanatory: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 1.04
## 2 2 2.35
## 3 3 10.1
## 4 4 2.99
## 5 5 5.29
## 6 6 5.26
## 7 7 1.40
## 8 8 2.58
## 9 9 1.12
## 10 10 5.32
## # ℹ 990 more rowsconf_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: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 1.04
## 2 2 2.35
## 3 3 10.1
## 4 4 2.99
## 5 5 5.29
## 6 6 5.26
## 7 7 1.40
## 8 8 2.58
## 9 9 1.12
## 10 10 5.32
## # ℹ 990 more rowsThe image’s distribution is right-skewed, which means that its median, mean, and mode are all smaller than 10. Furthermore, the confidence intervals fall between 0 and about 10, suggesting that this range is most likely where the actual value is. Additionally, because the red line is at 20 and does not fall into the null distribution, it provides evidence against the null hypothesis.
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.
table(france_data$lrscale)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88
## 945 499 1269 1996 1838 5274 1594 1749 1458 409 686 381 940table(france_data$ccrdprs)##
## 0 1 2 3 4 5 6 7 8 9 10 66 77 88 99
## 72 26 64 78 114 478 403 736 952 358 708 27 1 29 1france_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
),
ccrdprs = ifelse(ccrdprs %in% c(66, 77, 88, 99), NA, ccrdprs), # Set values 66, 77, 88, 99 in 'ccrdprs' to NA
)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 under comparison differ from one another. There is an obvious correlation between a person’s political stance (left or right) and their sense of personal accountability for addressing climate change. In other words, a person’s political stance essentially determines how much they believe has to be done to combat 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 %>% # replace name here if you assigned something other than null_distribution above
get_pvalue(obs_stat = test_stat, direction = "greater") # would only replace test_stat if assigned another name in Step 1## 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 0According to frequentist hypothesis testing assumptions, there is strong evidence to reject the null hypothesis with a p-value of 0.
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.644
## 2 2 0.138
## 3 3 0.609
## 4 4 0.367
## 5 5 1.15
## 6 6 1.11
## 7 7 0.110
## 8 8 0.787
## 9 9 1.29
## 10 10 0.314
## # ℹ 990 more rowsThe image’s distribution is right-skewed, which means that its median, mean, and mode are all smaller than five. Furthermore, it is possible that the true number falls within this range because the confidence intervals run from 0 to almost 5. Additionally, because the red line is slightly over 15 and does not fit into the null distribution, it offers evidence against the null hypothesis.