Set up:

# List of packages
packages <- c("tidyverse", "infer", "fst") 

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

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.

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 and Right.

france_data <- ess %>%
  # Filter for France
  filter(cntry == "FR") %>%
  
  # Recoding lfscale and educ.ba
  mutate(
    lfscale = 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% c(4, 5, 6, 77, 88, 99) ~ NA_character_    # Values 4, 5, 6, 77, 88, 99 in 'lrscale' are set as NA
    
),
 # 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" # anything else is 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
  )
unique(france_data$educ.ba)
## [1] "No BA"      "BA or more"
unique(france_data$lfscale)
## [1] "Left"  NA      "Right"
france_data <- france_data %>% filter(!is.na(educ.ba))
france_data <- france_data %>% filter(!is.na(lfscale))

# double-checking + quick distribution
unique(france_data$educ.ba)
## [1] "No BA"      "BA or more"
table(france_data$lfscale)
## 
##  Left Right 
##  4709  4302
tab_dat <- france_data

# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$lfscale, 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.1311243 0.3914592
## [2,] 0.1197912 0.3576253
# Table of expected frequencies
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
# 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

Interpretation: The code shows that the critical value is 3.841. So, if the Chi-square test statistic is greater than 3.841, it would indicate the null hypothesis (no group differences) can be rejected as one’s position on the political spectrum and education level may be correlated at the 0.05 significance level.

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: 100.8551

Interpretation: Since the Chi-squared value is 100.855, greater than the critical value of 3.841, there seems to be a statistically significant relationship between the probability of leaning politically left or right and various education levels (we can reject the null hypothesis).

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 = 100.86, df = 1, p-value < 2.2e-16

Interpretation: The p-value is <2,2e-16, which is very close to 0, suggesting there is sufficient evidence to reject the null hypothesis (which states the two variables are independent of 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 France
  filter(cntry == "FR") %>%
  
  # Recoding sbbsntx and domicil
  mutate(
    sbbsntx = case_when(
      sbbsntx == 1 ~ "Agree strongly",  
      sbbsntx == 2 ~ "Agree",   
      sbbsntx == 3 ~ "Neither agree or disagree",   
      sbbsntx == 4 ~ "Disagree",
      sbbsntx == 5 ~ "Disagree strongly",  
      sbbsntx %in% c(7, 8, 9) ~ NA_character_,  # excluding 7, 8 and 9
    ),

   
    domicil = case_when(
       domicil == '1' ~ "Urban", 
       domicil == '2' ~ "Urban",
       domicil == '3' ~ "Rural", 
       domicil == '4' ~ "Rural", 
       domicil == '5' ~ "Rural",
       domicil == '7' ~ NA_character_,
       domicil == '8' ~ NA_character_,
       domicil == '9' ~ NA_character_)
  )

Double-checking if variables are cleaned

table(france_data$sbbsntx) # it is clean because there are 5 categories 
## 
##                     Agree            Agree strongly                  Disagree 
##                      1538                       757                       640 
##         Disagree strongly Neither agree or disagree 
##                       301                       804
table(france_data$domicil) # it is cleaned because there are 2 categories
## 
## Rural Urban 
## 13087  5945

Infer - Step 1 (calculate test statistic that measures the difference between the observed sample mean and the population mean):

test_stat <- france_data %>%
  specify(explanatory = domicil,
          response = sbbsntx) %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "Chisq") # Chose Chi-square test of independence because domicil is a categorical variable with three or more categories and sbbsntx is a categorical variable with two categories.
## Warning: Removed 14999 rows containing missing values.
print(test_stat$stat)
## X-squared 
##  19.97187

Therefore, the test statistic is 19.97, meaning there may be a statistically significant relationship between domicil and sbbsntx as the larger the t-statistic (often higher than 2), the larger the difference between the groups- indicating greater evidence to reject the null hypothesis.

Step 2 (simulate null distribution):

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 (calculate p-value):

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

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

Therefore, the p-value is less than 0.05 (the significance level), suggesting there is a relationship between one’s geographic location (rural or urban/domicil) and one’s view on whether social services are too costly for businesses with regard to taxes (sbbsntx).

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 = "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.33
##  2         2  6.68
##  3         3  5.05
##  4         4  7.83
##  5         5  2.46
##  6         6  7.34
##  7         7  9.62
##  8         8  6.42
##  9         9  3.87
## 10        10  3.86
## # ℹ 990 more rows

Interpretation: The graph shows that there is sufficient evidence to reject the null hypothesis (no group differences) as the test statistic of 19.97 (red line) falls outside the blue area (signifies the null distribution with confidence intervals- where the test statistic would lie if the null hypothesis was supported).

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.

Provide interpretations of the output (consider the variable info).

france_data <- ess %>%
  # Filter for France
  filter(cntry == "FR") %>%
  
  # Recoding ccrdprs and lfscale 
  mutate(
    ccrdprs = case_when(
      ccrdprs == 0 ~ "0",  
      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",
      ccrdprs %in% c(66, 77, 88, 99) ~ NA_character_, 
 ),

  lfscale = 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% c(4, 5, 6, 77, 88, 99) ~ NA_character_)    # Values 4, 5, 6, 77, 88, 99 in 'lrscale' are set as NA
  )

Double-check if variables are cleaned:

table(france_data$ccrdprs)
## 
##   0   1  10   2   3   4   5   6   7   8   9 
##  72  26 708  64  78 114 478 403 736 952 358
table(france_data$lfscale)
## 
##  Left Right 
##  4709  4302

Infer - Step 1:

test_stat2 <- france_data %>%
  specify(explanatory = lfscale,
          response = ccrdprs) %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "Chisq") # I used the Chi-square test because the response variable refers to various categories (for instance 0 refers to "Not at all," the numbers would hold possible inherent value without this) and the explanatory variable is one categorical variable with two categories.
## Warning: Removed 17167 rows containing missing values.
print(test_stat2$stat)
## X-squared 
##  35.04425

Step 2:

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

Step 3:

p_val2 <- null_distribution2 %>% 
  get_pvalue(obs_stat = test_stat2, direction = "two-sided") # 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_val2
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0

Therefore, the p-value is 0, suggesting there is strong evidence that there is a relationship between one’s position on the political spectrum (left or right/lfscale) and one’s sense of personal responsibility towards climate change reduction (ccrdprs), as it is less than the significance level of 0.05.

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_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: ccrdprs (factor)
## Explanatory: lfscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1 17.3 
##  2         2  4.54
##  3         3  9.45
##  4         4 14.4 
##  5         5  5.75
##  6         6  5.52
##  7         7  3.93
##  8         8 10.1 
##  9         9 10.4 
## 10        10  4.34
## # ℹ 990 more rows

Interpretation: The graph shows that there is sufficient evidence to reject the null hypothesis (no group differences) as the test statistic of 35.04 (red line) falls outside the blue area signifying where the test statistic would lie if the null distribution was supported.