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")
   france_data <- ess %>%
  filter(cntry == "FR") %>%
  
  mutate(
     lrscale = case_when(
      lrscale %in% 0:4 ~ "Left",                    
      lrscale %in% 6:10 ~ "Right",                 
      lrscale %in% c(5, 77, 88, 99) ~ NA_character_    
    ),
     ccrdprs = case_when(
      ccrdprs %in% 0:5 ~ "not at all",                    
      ccrdprs %in% 6:10 ~ "A great deal",                 
      ccrdprs %in% c(66, 77, 88, 99) ~ NA_character_    
    ),
  
    clsprty = case_when(
      clsprty == 1 ~ "Yes",  
      clsprty == 2 ~ "No",   
      clsprty %in% c(7, 8, 9) ~ NA_character_, 
      TRUE ~ as.character(clsprty)  
    ),
      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_,  
    ),
    domicil = case_when(
      domicil %in% 0:2 ~ "urban",  
      domicil %in% 3:5 ~ "Rural",                  
      domicil %in% c(7, 8, 9) ~ 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)  
  )

`

unique(france_data$lrscale)
## [1] "Left"  NA      "Right"
france_data <- france_data %>% filter(!is.na(lrscale))

unique(france_data$lrscale)
## [1] "Left"  "Right"
table(france_data$lrscale)
## 
##  Left Right 
##  6547  5896

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.

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.1352292 0.390930
## [2,] 0.1217827 0.352058
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")
## Table of Expected Frequencies:
print(ftab)
##          [,1]     [,2]
## [1,] 1682.657 4864.343
## [2,] 1515.343 4380.657
alpha <- 0.05 
c_val <- qchisq(alpha, df = 1, lower.tail = FALSE)
cat("Critical Value:", round(c_val, 3), "\n")
## Critical Value: 3.841
#  chi-squared statistic:
test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")
## Pearson's X^2: 56.7436
# p-value:
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 = 56.744, df = 1, p-value = 4.965e-14

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.

THe critical value is 3.841, the X-squared value of 56.7 is much larger than this critical value. The difference between our observed data and the expected assumption of independence is statistically significant. Thus, there is a statistically significant association between Ideological self-placement and one’s educational level. This means that the probability of Ideological self-placement is not the same across different education levels.

The Pearson’s chi-squared statistic is 56.7 and the p-value is 4.965e-14. The Chi-squared test statistic exceeds 3.841 so so the null hypothesis is rejected at the 0.05 significance level. Meaning there is a statistically significant association between Ideological self-placement and education level in the dataset. The P-value being very close to 0 shows that the null hypothesis that the two variables are independent of eachother can be rejected. So there is evidence that there is an association between the two varibales.

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

test_stat <- france_data %>%
  specify(explanatory = domicil, 
          response = sbbsntx) %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "chisq") 
## Warning: Removed 9728 rows containing missing values.
print(test_stat$stat)
## X-squared 
##  20.05435
null_distribution <- france_data %>%
  specify(explanatory = domicil,
          response = sbbsntx) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "chisq")
## Warning: Removed 9728 rows containing missing values.
p_val <- null_distribution %>% 
  get_pvalue(obs_stat = test_stat, direction = "two-sided") 

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

The value of x2 of 20.05 indicates a large difference between the frequency of data. So there is a strong association between the variables. The p value of 0.002 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.

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 2.86 
##  2         2 3.68 
##  3         3 0.827
##  4         4 4.29 
##  5         5 3.97 
##  6         6 9.17 
##  7         7 4.71 
##  8         8 0.242
##  9         9 5.54 
## 10        10 3.97 
## # ℹ 990 more rows

The visualizes the null distribution alongside the observsed test statistic indicating evidence, 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).

other_stat <- france_data %>%
  specify(explanatory = lrscale, 
          response = ccrdprs, 
          success = "A great deal") %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "chisq") 
## Warning: Removed 9913 rows containing missing values.
print(other_stat$stat)
## X-squared 
##  7.826424
null_distribution4 <- france_data %>%
  specify(explanatory = domicil,
          response = sbbsntx) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "chisq")
## Warning: Removed 9728 rows containing missing values.
p_val4 <- null_distribution4 %>% 
  get_pvalue(obs_stat = other_stat, direction = "two-sided") 

p_val4
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1   0.204

The value of x2 of 7.8 indicates that there could be a weaker association between the variables and that there is less evidence to reject the null. A P value of 0.234 shows strong evidence against the null, however still leaves possibility that the null hypothesis is possible to occur but unlikely.

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


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

null_distribution4
## Response: sbbsntx (factor)
## Explanatory: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1 6.48 
##  2         2 3.79 
##  3         3 1.95 
##  4         4 1.49 
##  5         5 6.04 
##  6         6 0.571
##  7         7 7.98 
##  8         8 4.65 
##  9         9 0.767
## 10        10 9.39 
## # ℹ 990 more rows

The visualizes the null distribution alongside the observsed test statistic indicating evidence, that the null hypothesis can occur. The statistic is very close to the null distribution as seen through the shaded red area which shows the null hypothesis is still possible to occur.