# 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.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
## [[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",                  
      lrscale %in% c(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"
  )
)
# using unique to identify
unique(france_data$gndr)
## [1] 2 1
# does not return NA, but the following does
unique(france_data$clsprty)
## [1] 1 2 8 7
# to 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] 1 2 8 7
table(france_data$clsprty)
## 
##    1    2    7    8 
## 9507 9194  101  236
tab_dat <- france_data

# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$lrscale, 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
# 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
# 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

The Chi-squared test statistic (100.8551) exceeds 3.841, so we should reject the null hypothesis at the 0.05 significance level. This indicates that there is a statistically significant association between placement on the left-right scale (lrscale) and level of education (educ.ba) in the dataset. This means that an individuals’ political orientation appears to be linked to their educational attainment within this dataset.

The observed X-squared value (100.86) is much larger than the critical value. The difference between the observed data and what we would expect under the assumption of independence is statistically significant. Thus, there is a statistically significant association between placement on the left-right scale and one’s educational level. This means that placement on the left-right scale is not the same across different education levels.

The extremely small p-value (2.2e-16) means that we should reject the null hypothesis that the two variables are independent of each other. There is strong evidence to suggest that there is a significant association between placement on the left-right scale and level of education in the dataset.

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(
    geo = recode(as.character(domicil), 
                 '1' = "Urban", 
                 '2' = "Urban", 
                 '3' = "Rural", 
                 '4' = "Rural", 
                 '5' = "Rural",
                 '7' = NA_character_,
                 '8' = NA_character_,
                 '9' = NA_character_
    ),
    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(clsprty), !is.na(geo)) 

Step 1: Calculate the test statistic of sample

test_stat <- france_data %>%
  specify(explanatory = geo, 
          response = sbbsntx) %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "Chisq") 
## Warning: Removed 14993 rows containing missing values.
print(test_stat$stat)
## X-squared 
##  19.97187

The Chi-squared test statistic of 19.97187 suggests a significant association between individuals’ beliefs about if social benefits/services costs businesses too much in taxes/charges and where they live. This suggests that individuals’ opinions on the cost of social benefits/services to businesses vary significantly depending on their place of residence.

Step 2: Simulate the null distribution

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

Step 3: Calculate the p-value of sample

p_val <- null_distribution %>% 
  get_pvalue(obs_stat = test_stat, direction = "two-sided") 
## 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 p-value of 0 indicates that 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.

Step 4: Visualize

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: geo (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate   stat
##        <int>  <dbl>
##  1         1 7.10  
##  2         2 4.21  
##  3         3 4.72  
##  4         4 0.600 
##  5         5 1.49  
##  6         6 0.0996
##  7         7 2.15  
##  8         8 2.47  
##  9         9 2.55  
## 10        10 1.86  
## # ℹ 990 more rows
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: geo (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate   stat
##        <int>  <dbl>
##  1         1 7.10  
##  2         2 4.21  
##  3         3 4.72  
##  4         4 0.600 
##  5         5 1.49  
##  6         6 0.0996
##  7         7 2.15  
##  8         8 2.47  
##  9         9 2.55  
## 10        10 1.86  
## # ℹ 990 more rows

The visualization of the null distribution alongside the observsed test statistic, indicates 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(
    lrscale = case_when(
      lrscale %in% 0:3 ~ "Left",                    
      lrscale %in% 7:10 ~ "Right",                  
      lrscale %in% c(77, 88, 99) ~ NA_character_ 
    ),
      
    ccrdprs = case_when(
     ccrdprs == 0 ~ "Not at all",
     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 ~ "A great deal",
     ccrdprs %in% c(66, 77, 88, 99) ~ NA_character_,  
    TRUE ~ as.character(ccrdprs)
  ))

Step 1: Calculate the test statistic of sample

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

The Chi-squared test statistic of 35.04425 suggests a significant association between the extent to which individuals feel personal responsibility to reduce climate change and their placement on the left-right scale. This indicates that political orientation influences attitudes towards climate change responsibility.

Step 2: Simulate the null distribution

null_distribution2 <- 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: Calculate the p-value of sample

p_val2 <- null_distribution2 %>% 
  get_pvalue(obs_stat = test_stat2, direction = "two-sided") 
## 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

A p-value of 0 indicates that 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.

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

null_distribution2
## Response: ccrdprs (factor)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1 15.3 
##  2         2  8.24
##  3         3  4.82
##  4         4  8.69
##  5         5  6.93
##  6         6 11.8 
##  7         7  6.03
##  8         8 13.8 
##  9         9  6.15
## 10        10  5.63
## # ℹ 990 more rows
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: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1 15.3 
##  2         2  8.24
##  3         3  4.82
##  4         4  8.69
##  5         5  6.93
##  6         6 11.8 
##  7         7  6.03
##  8         8 13.8 
##  9         9  6.15
## 10        10  5.63
## # ℹ 990 more rows

The visualization of the null distribution alongside the observsed test statistic, indicates clear evidence, under assumptions of frequentist hypothesis testing, to reject the null.