# 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"
france_data <- read.fst("france_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 <- france_data %>%
  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$ID, 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.1311243 0.3914592
## [2,] 0.1197912 0.3576253
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
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
test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")
## Pearson's X^2: 100.8551
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

If the chi-squared statistic were to be greater than the critical value, there is a 5 percent signifigant level to reject the null hypothesis with strong evidence. From the curated values, the chi-statistic of 100.86 is greater than the critical value of 3.841 which means that there is a 5% or less probability that the chi-statistic would be observed if the null hypothesis were true. Thus, the null hypothesis can be rejected. Therefore from these outputed values it can be determined that there is an association between the variables of Domcil and Lrscale and are indepedent. Since the P-value ouputed is 2.23-16 this allows us with confidence to reject the null due to the alpha value being 0.05. Therefore from this observation, it can be determined that the variables of political party affiliation (left or right wing) and education level (Ba or no BA) are unlikey to have occured by chance but rather a meaningful association.

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(france_data$domicil)
## 
##    1    2    3    4    5    7    8 
## 3584 2361 6257 5668 1162    3    3
table(france_data$sbbsntx)
## 
##    1    2    3    4    5    7    8 
##  757 1538  804  640  301    3  100
france_data <- france_data%>%
  mutate(
    domicil = ifelse(domicil == 2, 0, ifelse(domicil %in% c(7, 8, 9), NA, clsprty)), # If 'clsprty' is 2, set it to 0. If it's 7, 8, or 9, set it to NA.
  sbbsntx = ifelse(sbbsntx %in% c(77, 88, 99), NA, sbbsntx),  
  )
france_data <- france_data %>%
  filter(!is.na(sbbsntx)) %>%
  filter(!is.na(domicil))
table(france_data$domicil)
## 
##    0    1    2    7    8 
##  487 1884 1705   35   31
table(france_data$sbbsntx)
## 
##    1    2    3    4    5    7    8 
##  756 1538  804  640  301    3  100
france_data <- france_data %>%
  mutate(
    domicil = recode(as.character(domicil), 
                 '1' = "Urban", 
                 '2' = "Urban",
                 '3' = "Rural", 
                 '4' = "Rural", 
                 '5' = "Rural",
                 '7' = NA_character_,
                 '8' = NA_character_,
                 '9' = NA_character_)
  ) %>%
  filter(!is.na(domicil))
france_data <- france_data %>%
  mutate(
    sbbsntx = recode(as.character(sbbsntx), 
                 '1' = "Agree Strongly", 
                 '2' = "Agree",
                 '3' = "Neither Agree or Disagree", 
                 '4' = "Disagree", 
                 '5' = "Disagree Strongly",
                 '7' = NA_character_,
                 '8' = NA_character_,
                 '9' = NA_character_)
  ) %>%
  filter(!is.na(sbbsntx))
test_stat <- france_data %>%
  specify(explanatory = domicil, 
          response = sbbsntx) %>%
  hypothesize(null = "independence") %>%
  calculate(stat = "chisq") 

print(test_stat$stat)
## X-squared 
##  4.055423

This value measures how much our observed data deviates from what we would expect if the variables of domicil and sbbsntx were independent. With the chi-squared statistic of 4.89 indicates that there is evidence against the null hypothesis. This means the variables are indepednet of one another

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 = test_stat, direction = "two-sided")

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

The p-value represents the probability of observing our data if domicil and sbbsntx are independent. Since the p-value is greater than 0.05 we fail to reject the null hypothesis meaning that there is likely no signifigant associaion between sbbsntx and domicilin the France ESS data. ## 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 = test_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  5.47 
##  2         2  1.75 
##  3         3  1.95 
##  4         4  2.72 
##  5         5 11.1  
##  6         6  0.558
##  7         7  4.43 
##  8         8  3.51 
##  9         9  6.02 
## 10        10  5.62 
## # ℹ 990 more rows
null_distribution
## Response: sbbsntx (factor)
## Explanatory: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate   stat
##        <int>  <dbl>
##  1         1  5.47 
##  2         2  1.75 
##  3         3  1.95 
##  4         4  2.72 
##  5         5 11.1  
##  6         6  0.558
##  7         7  4.43 
##  8         8  3.51 
##  9         9  6.02 
## 10        10  5.62 
## # ℹ 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: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate   stat
##        <int>  <dbl>
##  1         1  5.47 
##  2         2  1.75 
##  3         3  1.95 
##  4         4  2.72 
##  5         5 11.1  
##  6         6  0.558
##  7         7  4.43 
##  8         8  3.51 
##  9         9  6.02 
## 10        10  5.62 
## # ℹ 990 more rows

Since the red line that represents the chi-statistic, and falls within the distribution this graph demonstrates that we do not reject the null hypothesis. However if that red line were to fall outside of the distribution then we would reject the null hypothesis

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 <- read.fst("france_data.fst")
table(france_data$ccrdprs)
## 
##   0   1   2   3   4   5   6   7   8   9  10  66  77  88  99 
##  72  26  64  78 114 478 403 736 952 358 708  27   1  29   1
table(france_data$lrscale)
## 
##    0    1    2    3    4    5    6    7    8    9   10   77   88 
##  945  499 1269 1996 1838 5274 1594 1749 1458  409  686  381  940
france_data <- france_data %>%
  mutate(
    ccrdprs = ifelse(ccrdprs %in% c(66, 88, 99), NA, ccrdprs)
  )
france_data <- france_data %>%
    mutate(lrscale = ifelse(lrscale %in% 77:99, NA, lrscale)) %>%
    filter(!lrscale %in% 4:6) %>%
  mutate(
    lrscale = case_when(
      lrscale %in% 0:3 ~ "Left",                   
      lrscale %in% 7:10 ~ "Right",                  
      lrscale %in% c(77, 88, 99) ~ NA_character_ 
    ))
france_data <- france_data %>% filter(!is.na(ccrdprs))
france_data <- france_data %>% filter(!is.na(lrscale))
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 
## 5.482311
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.
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

Through the outcomes executed, it can be determined that due to the P value being zero, this suggests to reject the null hypothesis. The t statistic is intended to represent the difference between the two variables of interest. Since the t statistic in this case is 5.48 which is realitvely large, this suggest that there is a big bifference is means between the two groups lrscale and ccrdprs.

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 = test_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 -0.637
##  2         2  0.102
##  3         3 -0.639
##  4         4 -0.617
##  5         5 -0.596
##  6         6 -0.406
##  7         7 -0.554
##  8         8  1.97 
##  9         9  1.06 
## 10        10 -0.469
## # ℹ 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)

null_distribution
## Response: ccrdprs (numeric)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate   stat
##        <int>  <dbl>
##  1         1 -0.637
##  2         2  0.102
##  3         3 -0.639
##  4         4 -0.617
##  5         5 -0.596
##  6         6 -0.406
##  7         7 -0.554
##  8         8  1.97 
##  9         9  1.06 
## 10        10 -0.469
## # ℹ 990 more rows

The red line that represents the test statistic. Since it is outside of the confidence interval we reject the null hypothesis.