# 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 for France
  filter(cntry == "FR") %>%
  
  # Recoding Lrsscale and cleaning trstplt (will use later)
  mutate(
    ID = 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% 4:6 ~ "Moderate",                # Values 4 to 6 in 'lrscale' are categorized as "Moderate"
      lrscale %in% c(77, 88, 99) ~ NA_character_    # Values 77, 88, 99 in 'lrscale' are set as NA
    ),
    trstplt = ifelse(trstplt %in% c(77, 88, 99), NA, trstplt),
    # 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"
    )
  ) %>%
  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), 
  ) %>%

  filter(!is.na(ID) & !is.na(educ.ba))

Long form coding for Chi-Square and calculating/interpreting critical value for independence

tab_dat <- france_data

# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$ID, 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)
# 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.06135806 0.2044319
## [2,] 0.11343879 0.3779537
## [3,] 0.05605487 0.1867628
# Table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")
## Table of Expected Frequencies:
print(ftab)
##           [,1]     [,2]
## [1,] 1087.0808 3621.919
## [2,] 2009.7951 6696.205
## [3,]  993.1241 3308.876
# 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

The critical value refers tot he cutoof value for the Chi-Square, given the significance level of 0.05. Therefore, a cirtical value of 3.81 means that the Chi-Square calculated above this value would lead to rejection of the null hypothsis meaning that there is a relationship between the variables.

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

#step 1: calculate the test statistics of sample

france_data <- france_data %>%
  # Filter for France
  filter(cntry == "FR") %>%
  
  mutate(
     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 ~ NA_character_   
      
    ),# Recoding domicil
    domicil = case_when(
      domicil == 1 ~ "Big city",  
      domicil == 2 ~ "Suburbs",
      domicil == 3 ~ "Small City",
      domicil == 4 ~ "Village",
      domicil == 5 ~ "Country Side",
      domicil %in% c(7, 8, 9) ~ NA_character_,
      TRUE ~ NA_character_  ))%>%

  filter(!is.na(sbbsntx) & !is.na(domicil))
test_stat <- france_data %>%
  specify( response = sbbsntx, 
          explanatory = domicil) %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "chisq") 

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

#step 2: simulate the null distribution

null_distribution <- france_data %>%
  specify(explanatory = sbbsntx,
          response = domicil) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "Chisq")

#step 3: 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.02

Therefore, the p vlaue is low and the null evidence is rejected, therefore there is a low likelyhood of the differnces in groups happened due to random variability.

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: domicil (factor)
## Explanatory: sbbsntx (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1 15.9 
##  2         2 12.1 
##  3         3 12.7 
##  4         4 14.9 
##  5         5 13.0 
##  6         6 44.0 
##  7         7 18.4 
##  8         8 15.0 
##  9         9  4.07
## 10        10 16.3 
## # ℹ 990 more rows

Therefore, since the obsereved statisitc is on the far right, this indicates that the relationships between the varibales is very significant. As well, the line being on the far right indicates null distribution. Statistically, the line indicates a low p-value as well meaning that the observered factors are not random variation.

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 <- france_data %>%
  filter(cntry == "FR") %>%
  # Recoding Lrsscale
  mutate(
    ID = case_when(
      lrscale == 0 ~ "Left",
      lrscale == 1 ~ "Left",
      lrscale == 2 ~ "Left", 
      lrscale == 3 ~ "Left",  
      lrscale == 4 ~ "Moderate",  
      lrscale == 5 ~ "Moderate",
      lrscale == 6 ~ "Moderate",
      lrscale == 7 ~ "Right", 
      lrscale == 8 ~ "Right",
      lrscale == 9 ~ "Right",
      lrscale == 10 ~ "Right", 
      lrscale %in% c(77, 88, 99) ~ NA_character_    # Values 77, 88, 99 in 'lrscale' are set as NA
    ),
    # Recoding domicil
    ccrdprs = case_when(
      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",
      TRUE ~ NA_character_ 
    )
  )%>%
  
filter(!is.na(ID) & !is.na(ccrdprs))
null_lrs <- france_data %>%
  specify(explanatory = ID, 
          response = ccrdprs) %>%
  hypothesize(null = "independence") %>%
  calculate(stat = "chisq") 
print(null_lrs$stat)
## X-squared 
##  31.35275

Such a large t-statisitc indicates a signifcant difference between variables, therefore the greater the t-statisitc, the greater for the chances of a null hypothesis

#step 2: simulate the null

null_dist <- france_data %>%
  specify(explanatory = ID,
          response = ccrdprs) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "Chisq")

#step 3: p-value

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

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

Therefore,the p-value of 0.014 inidcates that the relationship between the groups is unlikely to be random, as well since the p value is low, it also suggests that the null is rejected.

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

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

null_dist
## Response: ccrdprs (factor)
## Explanatory: ID (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1  19.8
##  2         2  15.9
##  3         3  21.0
##  4         4  11.2
##  5         5  29.2
##  6         6  20.6
##  7         7  23.3
##  8         8  13.0
##  9         9  17.4
## 10        10  18.2
## # ℹ 990 more rows

Therefore, the Chi-square is once again far to the right of the null distribution. this means that the observed statisitcs are greater than what is expected if there were no connection between the variables, ‘ccrdprs’ and ‘irscale’. This suggests that there is a connection between the variables which is supported by the 95% confidence interval as the observed variable moves aways from the assumption of independence between the variables. The p-value is very low, once again suggesting the null hypothesis. To conclude, the graph emphasizes the statistical significance of the observed differences.