packages <- c("tidyverse", "infer", "fst") 

new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_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.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] "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.3576253
ftab <- 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.561
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

The critical value calculated is 3.841. This means that if the Chi-squared test statistic is greater than 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 being left or right and education level in my dataset.

test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")
## Pearson's X^2: 100.8551

The observed X-squared value is 100.8551. 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 being left or right and one’s educational level. This means that the probability of being left-wing or right-wing is not the same across different education levels.

p_val <- pchisq(test_stat, df = 1, lower.tail = F)
cat("p-value:", round(p_val, 4), "\n")
## p-value: 0
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 = 100.86, df = 1, p-value < 2.2e-16

The observed p value is < 2.2e-16. 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 100.8551 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).

france_data <- ess %>%
  filter(cntry == "FR") %>%
  
 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)  
    ),
    
      domicil = case_when(
      domicil == 1 ~ "Big city",  
      domicil == 2 ~ "Suburbs",
      domicil == 3 ~ "Town",
      domicil == 4 ~ "Village",
      domicil == 5 ~ "Farm",
      domicil %in% c(7, 8, 9) ~ NA_character_, 
      TRUE ~ as.character(domicil) 
     )
    
 )

Step 1:

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 
##  36.01302

Step 2:

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.

Step 3:

p_val <- null_distribution %>% 
  get_pvalue(obs_stat = test_stat, direction = "greater") 

p_val
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1   0.004

A t-statistic of 36.01302 is quite large, indicating a significant difference in beliefs on social benefits costs between areas in which people live. It suggests that the mean belief on cost of social services is different for the multiple groups in different areas.

The p-value of 0.001, under the assumptions of frequentist hypothesis testing, 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.

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)

null_distribution
## Response: sbbsntx (factor)
## Explanatory: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1 14.8 
##  2         2 23.8 
##  3         3 16.5 
##  4         4  7.06
##  5         5 14.4 
##  6         6 17.7 
##  7         7  9.00
##  8         8 19.6 
##  9         9  9.17
## 10        10 12.2 
## # ℹ 990 more rows

The observed test statistic (red line) falls outside the null distribution and confidence interval indicating here clear evidence, under assumptions of frequentist hypothesis testing, to reject the null.

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(cntry == "FR") %>%
  
 mutate(
    ccrdprs = case_when(
      ccrdprs %in% c(66, 77, 88, 99) ~ NA_character_, 
      TRUE ~ as.character(ccrdprs)  
    ),
    
     lrscale = case_when(
      lrscale %in% 0:3 ~ "Left",      
      lrscale %in% 7:10 ~ "Right",     
      TRUE ~ NA_character_  
     )
    
 )

Step 1:

cc_stat <- france_data %>%
  specify(explanatory = lrscale, 
          response = ccrdprs) %>%
  hypothesize(null = "independence") %>%
  calculate(stat = "Chisq") 
## Warning: Removed 17167 rows containing missing values.
print(cc_stat$stat)
## X-squared 
##  35.04425

Step 2:

null_cc <- france_data %>%
  specify(explanatory = lrscale, 
          response = ccrdprs) %>%
  hypothesize(null = "independence") %>%
    generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "Chisq") 
## Warning: Removed 17167 rows containing missing values.

Step 3:

p_val <- null_cc %>% 
  get_pvalue(obs_stat = cc_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       0

A t-statistic of 35.04425 is quite large, indicating a significant difference in political allignement and between responsibility to reduce climate change. It suggests that the mean belief on reducing climate change is different for those that are left or right.

The p-value of 0, under the assumptions of frequentist hypothesis testing, 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.

conf_int <- null_distribution %>%
  get_confidence_interval(level = 0.95, type = "percentile")


null_cc %>%
  visualize() +
  shade_p_value(obs_stat = cc_stat, direction = "greater") +
  shade_confidence_interval(endpoints = conf_int)
## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf

null_cc
## Response: ccrdprs (factor)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1 12.9 
##  2         2 16.8 
##  3         3  5.69
##  4         4  4.58
##  5         5  6.45
##  6         6  9.53
##  7         7  4.73
##  8         8  6.57
##  9         9 11.5 
## 10        10 13.3 
## # ℹ 990 more rows

The observed test statistic (red line) falls outside the null distribution and confidence interval indicating here clear evidence, under assumptions of frequentist hypothesis testing, to reject the null.

End.