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.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")
france_data <- ess %>%
  filter(cntry == "FR") %>%
  
   mutate(
    lrscale = case_when(
      lrscale %in% c(0, 1, 2, 3) ~ "Left",
      lrscale %in% c(7, 8, 9, 10) ~ "Right",
      lrscale %in% c(77, 88, 99) ~ NA_character_
    ),
    
    gndr = case_when(
      gndr == 1 ~ "Male",
      gndr == 2 ~ "Female",
      gndr == 9 ~ NA_character_,
      TRUE ~ as.character(gndr)
    ),

    educ.ba = case_when(
      essround < 5 & edulvla == 5 ~ "BA or more",
      essround >= 5 & edulvlb > 600 ~ "BA or more",
      TRUE ~ "No BA"
    ),
  
    sbbsntx = case_when(
      sbbsntx == 1 ~ "Agree strongly",
      sbbsntx == 2 ~ "Agree",
      sbbsntx == 3 ~ "Neither agree nor disagree",
      sbbsntx == 4 ~ "Disagree",
      sbbsntx == 5 ~ "Strongly disagree",
      sbbsntx %in% c(7, 8, 9) ~ NA_character_,
      TRUE ~ as.character(sbbsntx)
    ),
   
     domicil = case_when(
      domicil %in% c(1, 2) ~ "Urban",
      domicil %in% c(3, 4, 5) ~ "Rural",
      domicil %in% c(7, 8, 9) ~ NA_character_,
      TRUE ~ as.character(domicil)
    ),
    
    ccrdprs = ifelse(ccrdprs %in% c(66, 77, 88, 99), NA, ccrdprs),

)
france_data <- france_data %>% filter(!is.na(lrscale))
france_data <- france_data %>% filter(!is.na(gndr))
france_data <- france_data %>% filter(!is.na(educ.ba))
france_data <- france_data %>% filter(!is.na(sbbsntx))
france_data <- france_data %>% filter(!is.na(domicil))
france_data <- france_data %>% filter(!is.na(ccrdprs))

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.09203083 0.3805858
## [2,] 0.10269534 0.4246880
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")
## Table of Expected Frequencies:
print(ftab)
##           [,1]     [,2]
## [1,]  90.74239 375.2576
## [2,] 101.25761 418.7424
alpha <- 0.05
c_val <- qchisq(alpha, df = 1, lower.tail = FALSE)
cat("Critical Value:", round(c_val, 3), "\n")
## Critical Value: 3.841
test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")
## Pearson's X^2: 19.2795
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 = 19.279, df = 1, p-value = 1.129e-05

The critical value is 3.841. If the chi-squared test statistic is greater than 3.841, we reject the null hypothesis at the 0.05 significance level. This means a statistically significant association between political placement and education level. If the chi-squared test statistic is less than or equal to 3.841, we do not reject the null hypothesis at the 0.05 significance level. This would mean no statistically significant association between political placement and education level. However, since the Pearson’s chi-squared statistic is 19.280, it exceeds the critical value of 3.841, suggesting evidence against the null hypothesis. Therefore, we reject the null hypothesis and conclude a statistically significant association between political placement and education level. This means political placement is not the same across different education levels. The p-value is 1.129e-05, essentially 0. This means we reject the null hypothesis, which provides evidence to suggest a statistically significant association between political placement and education level.

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

chisq_stat <- france_data %>%
  specify(explanatory = domicil,
          response = sbbsntx) %>%
  hypothesize(null = "independence") %>%
  calculate(stat = "chisq")

print(chisq_stat$stat)
## X-squared 
##  11.16404
null_distribution <- france_data %>%
  specify(explanatory = domicil,
          response = sbbsntx) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "chisq")
p_val <- null_distribution %>%
  get_pvalue(obs_stat = chisq_stat, direction = "two-sided")

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

The chi-squared test statistic calculated for the relationship between domicile and sbbsntx is 11.164. This suggests some association between area of residence and opinion on social benefits/services cost since it is greater than the critical value of 3.841. However, with a p-value of 0.062, we fail to reject the null hypothesis at the 0.05 significance level. This indicates that there is insufficient evidence to claim a statistically significant relationship between area of residence and opinion on social benefits/services costs.

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.

null_distribution %>%
  visualize() +
  shade_p_value(obs_stat = chisq_stat, direction = "two-sided")
## 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.24 
##  2         2  0.988
##  3         3  2.77 
##  4         4  4.73 
##  5         5  4.43 
##  6         6 12.1  
##  7         7  5.87 
##  8         8  3.68 
##  9         9  1.76 
## 10        10  0.884
## # ℹ 990 more rows
conf_int <- null_distribution %>%
  get_confidence_interval(level = 0.95, type = "percentile")


null_distribution %>%
  visualize() +
  shade_p_value(obs_stat = chisq_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.24 
##  2         2  0.988
##  3         3  2.77 
##  4         4  4.73 
##  5         5  4.43 
##  6         6 12.1  
##  7         7  5.87 
##  8         8  3.68 
##  9         9  1.76 
## 10        10  0.884
## # ℹ 990 more rows

The null distribution visual shows a chi-squared test for independence between the categorical variables, ‘domicil’ and ‘sbbsntx.’ It includes green shading to demonstrate the 95% confidence interval, offering a range of values within which we are reasonably confident the true effect, the association between ‘domicil’ and ‘sbbsntx,’ lies. Since the p-value of 0.062 is large, the null distribution shows shaded red regions, which depict outcomes as extreme or more extreme than the observed test statistic. The confidence interval of 3.841 and the observed test statistic of 11.164 overlap with the null distribution, suggesting insufficient evidence to reject the null hypothesis. This means the relationship between area of residence and opinion on social benefits/services cost is likely due to random chance rather than a true association between the variables.

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

t_stat <- france_data %>%
  specify(explanatory = lrscale,
          response = ccrdprs) %>%
  hypothesize(null = "independence") %>%
  calculate(stat = "t", order = c("Left", "Right"))

print(t_stat$stat)
##        t 
## 4.855399
null_distribution <- france_data %>%
  specify(explanatory = lrscale,
          response = ccrdprs) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "t", order = c("Left", "Right"))
p_val <- null_distribution %>%
  get_pvalue(obs_stat = t_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

The t-test statistic calculated for the relationship between lrscale and ccrdprs is 4.855. This suggests some association between political placement and belief in personal responsibility to reduce climate change since it is greater than the critical value of 3.841. Additionally, with a p-value of 0, there is strong evidence to reject the null hypothesis at the 0.05 significance level. This indicates that there is sufficient evidence to claim a statistically significant relationship between political placement and belief in personal responsibility to reduce climate change.

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_distribution %>%
  visualize() +
  shade_p_value(obs_stat = t_stat, direction = "two-sided")

null_distribution
## Response: ccrdprs (numeric)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate    stat
##        <int>   <dbl>
##  1         1  1.25  
##  2         2 -0.485 
##  3         3  0.588 
##  4         4 -0.843 
##  5         5 -0.0160
##  6         6  0.894 
##  7         7 -0.512 
##  8         8  0.891 
##  9         9  1.06  
## 10        10  0.423 
## # ℹ 990 more rows
conf_int <- null_distribution %>%
  get_confidence_interval(level = 0.95, type = "percentile")


null_distribution %>%
  visualize() +
  shade_p_value(obs_stat = t_stat, direction = "two-sided") +
  shade_confidence_interval(endpoints = conf_int)

null_distribution
## Response: ccrdprs (numeric)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate    stat
##        <int>   <dbl>
##  1         1  1.25  
##  2         2 -0.485 
##  3         3  0.588 
##  4         4 -0.843 
##  5         5 -0.0160
##  6         6  0.894 
##  7         7 -0.512 
##  8         8  0.891 
##  9         9  1.06  
## 10        10  0.423 
## # ℹ 990 more rows

The null distribution visual shows a t-test for independence between the numerical and categorical variables, ‘lrscale’ and ‘ccrdprs,’ respectively. It includes green shading to demonstrate the 95% confidence interval, offering a range of values within which we are reasonably confident the true effect, the association between ‘lrscale’ and ‘ccrdprs,’ lies. Since the p-value of 0 is small, the null distribution shows a singular shaded red region, which depicts outcomes as extreme or more extreme than the observed test statistic. The observed test statistic of 4.855 does not overlap with the confidence interval of 3.841 and the null distribution, suggesting sufficient evidence to reject the null hypothesis. This means the relationship between political placement and belief in personal responsibility to reduce climate change is due to a true association between the variables rather than random chance.