Setting up environment

# To ensure there are no problems with the "Vector Memory" error, this code ensures R clears its environment which may help free up some memory
rm(list=ls()); gc() 
##          used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 527616 28.2    1172772 62.7         NA   669402 35.8
## Vcells 972131  7.5    8388608 64.0      16384  1851829 14.2
# 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)

# Loading 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"
#Loading the data
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.

table(ess$lrscale)
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  16299  10633  24208  41241  41150 139328  40165  44145  34501  11695  18027 
##     77     88     99 
##   9428  57354   2381
## Recoding lrscale for "left" and "right" and omitting "moderates" from Tutorial 4, and educ.ba in Tutorial 6

tab_dat <- ess %>% #Here we rename ess to tab_dat for when we do operations I can backtrack to clean dataset
  filter(cntry == "FR")%>%
  mutate(
    # Recoding 'lrscale' into a new variable 'ID'
    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",                # Omitted Values 4 to 6 in 'lrscale' are categorized as "Moderate"
      lrscale %in% c(4, 5, 6, 77, 88, 99) ~ NA_character_    # Values 77, 88, 99 in 'lrscale' are set as NA
    ),
      # 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"
    )) %>%
      # Ensure all NAs are handled uniformly
  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)  # Force educ.ba to be treated as character
  )

# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$clsprty, tab_dat$educ.ba)
# This is done so that the cross table is saved alongside all recoding from lrscale and educ.ba, and to ensure clarity when reading 
# Now, we do the long form coding for Chi square
# As we will be solving for expected proportions, expected frequencies, and the critical value for independence, it is best to organize the output for these numeric values for ease of reading

# 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:
# Generate table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")
## Table of Expected Frequencies:
# Generate 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
# Here we print the Table of Expected Proportions
print(ptab)
##             [,1]        [,2]
## [1,] 0.111530617 0.387839065
## [2,] 0.107858682 0.375070197
## [3,] 0.001184873 0.004120306
## [4,] 0.002768615 0.009627645
# Now we print the Expected Frequences
print(ftab)
##            [,1]       [,2]
## [1,] 2123.31989 7383.68011
## [2,] 2053.41359 7140.58641
## [3,]   22.55762   78.44238
## [4,]   52.70890  183.29110
# Now we print  the Critical Value for Independce
cat("Critical Value:", round(c_val, 3), "\n")
## Critical Value: 3.841
# Critical Vlaue is 3.841

My interpretation is that if the critical value for independance surpasses 3.841, then we will have to reject the null hypothesis at 0.05 significance. Therefore there is a statistically significant association between feeling close to a party and one’s educational level. This means that the probability of feeling close to a party is not the same across different educaiton levels.

# As we have just solved for the Critical Value for Independence, now it is time to determine and interpret the Pearson's chi-squared statistic and the p-value. We start with Pearson's where we then follow up with the p-value. 
# 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: 169.5185
# Pearson's X^2: 169.5185
# Now, we do the p-value, as well as using chisq.test for validation
# 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 = 169.52, df = 3, p-value < 2.2e-16
# p-value < 2.2e-16

The main takeaways is that p-value < 2.2e-16, in scientific notation, is extremely close to 0. The chi-squared test tells us if two categorical variables are independent or if they area associated in some manner. This very small p-value means we reject the null hypothesis that the two varaibles are independent of each other. The large chi-squared statistic of 169.52 reinforces this point.

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 statistic of your sample
table(ess$sbbsntx) #using the table(ess$) but inputting the variable sbbsntx to gain insight into its table
## 
##     1     2     3     4     5     7     8     9 
##  6577 27773 26302 25689  5868   122  8707   101
# It displays 9 categories, which we will need to make 7 , 8, 9 as NA
# In task one 'mytab' was the nanme used to save the crosstable in, for simplicity reasons in Task 2, we shall use the name 'mytab2' to represent task 2 and for organizational purposes
mytab2 <- tab_dat %>%
  # Filter for France
  filter(cntry == "FR") %>%
# Now, we sub in the variables as we did earlier for edulva 
  mutate(sbbsntx = ifelse(sbbsntx %in% c(7, 8, 9), NA, sbbsntx),
         domicil = ifelse(sbbsntx %in% c(7, 8, 9), NA, domicil))%>%
  # The %>% is key for the code to work, without it 'domicil' is unanble to be found
 # Here, for variable domicil, we indicate what integers shall represent urban and whihc shall represent rural
# Also recodes education with case_when for consistency
 mutate(
         domicil_ur = case_when(
         domicil < 3 ~ "urban",
         TRUE ~ "rural"))
# Now it is time to calculate the test statistic of the sample with  test_stat
test_stat <- mytab2 %>%
  specify(explanatory = domicil_ur, # change variable name for explanatory variable
          response = sbbsntx) %>% # change variable name for outcome of interest
  hypothesize(null = "independence") %>%
  calculate(stat = "t") # replace in between quotation marks appropriate test statistic
## Warning: Removed 14998 rows containing missing values.
## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "rural" - "urban", or divided in the order "rural" / "urban" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("rural", "urban")` to the calculate() function.
print(test_stat$stat)
##         t 
## -4.138823
# -4.138823

The t- statistic measures the difference between the observed sample mean and the population mean (or another sample mean, in the case of two-sample tests) in units of standard error. The larger the magnitude of t (regardless of the sign), the greater the evidence against the null hypothesis.

My interpretation is that due to the t-statistic being -4.138823 showcases a large difference.

# Step 2: Simulate the null doistribution
null_distribution <- mytab2 %>%
  specify(explanatory = domicil_ur,
          response = sbbsntx) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>% # only line we are adding, use reps = 1000 as standard but change if too computationally demanding
  calculate(stat = "t")
## Warning: Removed 14998 rows containing missing values.
## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "rural" - "urban", or divided in the order "rural" / "urban" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("rural", "urban")` to the calculate() function.
# Step 3: Calculate the p-value of your sample
p_val <- null_distribution %>% # replace name here if you assigned something other than null_distribution above
  get_pvalue(obs_stat = test_stat, direction = "two-sided") # would only replace test_stat if assigned another name in Step 1
## 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
# p-value is 0

My interpretation is that the p-value is 0, so we can reject the null hypothesis.

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

null_distribution
## Response: sbbsntx (numeric)
## Explanatory: domicil_ur (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate   stat
##        <int>  <dbl>
##  1         1  1.89 
##  2         2 -0.357
##  3         3  1.06 
##  4         4 -0.853
##  5         5  1.54 
##  6         6  1.60 
##  7         7  1.00 
##  8         8  0.558
##  9         9  1.40 
## 10        10  0.770
## # ℹ 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: sbbsntx (numeric)
## Explanatory: domicil_ur (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate   stat
##        <int>  <dbl>
##  1         1  1.89 
##  2         2 -0.357
##  3         3  1.06 
##  4         4 -0.853
##  5         5  1.54 
##  6         6  1.60 
##  7         7  1.00 
##  8         8  0.558
##  9         9  1.40 
## 10        10  0.770
## # ℹ 990 more rows

My interpretation of this graph is that due to the red line being sufficiently far away from the center of the null distribution, where it falls on the left tail, it therefore proves that we can 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).

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.