packages <- c("tidyverse", "infer", "fst") # add any you need here

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.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")

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% c(4,5,6,77, 88, 99) ~ 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(77, 88), NA_integer_, edulvla),
    edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
    educ.ba = as.character(educ.ba)  
  )
tab_dat <- france_data


mytab <- table(tab_dat$clsprty, tab_dat$educ.ba)


rsums <- rowSums(mytab)


csums <- colSums(mytab)


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.111530617 0.387839065
## [2,] 0.107858682 0.375070197
## [3,] 0.001184873 0.004120306
## [4,] 0.002768615 0.009627645
# Table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")
## Table of Expected Frequencies:
print(ftab)
##            [,1]       [,2]
## [1,] 2123.31989 7383.68011
## [2,] 2053.41359 7140.58641
## [3,]   22.55762   78.44238
## [4,]   52.70890  183.29110
# 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: 169.5185
# 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 = 169.52, df = 3, p-value < 2.2e-16

Discuss: Pearson’s chi-sqared statistic is 166.1849 is larger than p- value, p-value is extramely close 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).

table(ess$sbbsntx)
## 
##     1     2     3     4     5     7     8     9 
##  6577 27773 26302 25689  5868   122  8707   101
table(ess$domicil)
## 
##      1      2      3      4      5      7      8      9 
## 107787  55216 149256 148750  27534    137    493   1382
france_data2 <- 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)
  )
) %>%
    filter(!is.na(sbbsntx))

unique(france_data2$sbbsntx)
## [1] "Agree Strongly"             "Neither Agree Nor Disagree"
## [3] "Agree"                      "Disagree Strongly"         
## [5] "Disagree"
france_data2 <- france_data2 %>%
   mutate(geo = recode(as.character(domicil), 
                      '1' = "Urban", 
                      '2' = "Urban", 
                      '3' = "Rural", 
                      '4' = "Rural", 
                      '5' = "Rural",
                      '7' = NA_character_,
                      '8' = NA_character_,
                      '9' = NA_character_))

france_data2 <- france_data2 %>% 
  filter(!is.na(geo))

unique(france_data2$geo)
## [1] "Rural" "Urban"
test_stat2 <- france_data2 %>%
  specify(explanatory = geo, 
          response = sbbsntx) %>% 
  hypothesize(null = "independence") %>%
  calculate (stat = "Chisq")

print(test_stat2$stat)
## X-squared 
##  19.97187
null_distribution2 <- france_data2 %>%
  specify(explanatory = geo,
          response = sbbsntx) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "Chisq")
p_val2 <- null_distribution2 %>% # replace name here if you assigned something other than null_distribution above
  get_pvalue(obs_stat = test_stat2, direction = "two-sided") # would only replace test_stat if assigned another name in Step 1

p_val2
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1   0.002

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_int2 <- null_distribution2 %>%
  get_confidence_interval(level = 0.95, type = "percentile")


null_distribution2 %>%
  visualize() +
  shade_p_value(obs_stat = test_stat2, direction = "two-sided") +
  shade_confidence_interval(endpoints = conf_int2)
## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.

null_distribution2
## Response: sbbsntx (factor)
## Explanatory: geo (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1 9.27 
##  2         2 4.25 
##  3         3 5.19 
##  4         4 7.96 
##  5         5 3.19 
##  6         6 1.92 
##  7         7 4.09 
##  8         8 0.444
##  9         9 1.93 
## 10        10 5.91 
## # ℹ 990 more rows

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).

table(ess$ccrdprs)
## 
##     0     1     2     3     4     5     6     7     8     9    10    66    77 
##  6294  2668  4530  5706  5434 15489 12049 16523 15951  5891  9098  1939    79 
##    88    99 
##  1792   629
france_data <- ess %>%
  filter(cntry == "FR") %>%
  
  mutate(
    lrscale = case_when(
      lrscale %in% 0:3 ~ "Left",
      lrscale %in% 7:10 ~ "Right",
      TRUE ~ NA_character_,
      
      ccrdprs %in% c(66, 77, 88, 99) ~ NA_character_,
    )
  ) %>%
  filter(!is.na(lrscale), !is.na(ccrdprs))
test_stat <- france_data %>%
  specify(explanatory = lrscale, 
          response = ccrdprs) %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "t") 
## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "Left" - "Right", or divided in the order "Left" / "Right" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Left", "Right")` to the calculate() function.
print(test_stat$stat)
##         t 
## 0.5130018
null_distribution <- france_data %>%
  specify(explanatory = lrscale, 
          response = ccrdprs) %>% 
  hypothesize(null = "independence") %>%
    generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "t") 
## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "Left" - "Right", or divided in the order "Left" / "Right" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Left", "Right")` to the calculate() function.

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 = "two-sided") +
  shade_confidence_interval(endpoints = conf_int)