R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

Loading the list of packages:

# 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.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
## Warning: package 'infer' was built under R version 4.3.3
## [[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")

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.

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

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.

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.

Task1:

Loading data for France and cleaning variables of interest:

france_data<-ess %>%
  filter(cntry == "FR")%>%
   mutate(
    ID = case_when(
      lrscale %in% 0:3 ~ "Left",# Values 0-3 in are categorized "Left"
      lrscale %in% 7:10 ~ "Right",# Values 7-10 as "Right"
      # Values 4 to 6 in 'lrscale' are categorized as "Moderate" so Ommitted
      lrscale %in% c(4,5,6,77, 88, 99) ~ NA_character_# Values 77, 88, 99 in 'lrscale' are set as NA
    ),
     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
  )

#cross table created 

mytab <- table(france_data$ID, france_data$educ.ba)

Table of expected proportions and frequencies:

# 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:
print(ptab)
##           [,1]      [,2]
## [1,] 0.1311243 0.3914592
## [2,] 0.1197912 0.3576253
# Table of expected frequencies
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

critical value:

# 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

If the Chi-squared statistic is greater than the critical value (3.841 in this case), at the significance level 0.05 we will reject the null hypothesis.

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: 100.8551

Compared to the critical value, the Chi-squared score is much greater. This implies that that a strong statistical association exists between the educational levels (BA or No BA) and feeling close to a party (left and right). The probability of feeling close to a party varies with education levels. A significant difference exists in what we would expect under the assumption of independence as opposed to the observed data.

Calculating the 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 = 100.86, df = 1, p-value < 2.2e-16

A p-value of < 2.2e-16 is almost, although not exactly, zero which shows strong evidence exists to reject the null hypothesis at 5% significance level.

task 2:

Loading and cleaning variables of interest:

france_data2<-france_data%>%
  mutate(sbbsntx= ifelse(sbbsntx%in% c(7, 8, 9), NA, sbbsntx),
         domicil= ifelse(sbbsntx%in% c(7, 8, 9), NA, domicil))%>%
  mutate( # setting to urban and rural instead of peri-urban
        domicil_ur = case_when(
        domicil < 3 ~ "urban",
        TRUE ~ "rural"))

Step 1: Calculating test statistic of the sample at hand

test_stat <- france_data2 %>%
  specify(explanatory = domicil_ur, #explanatory variable
          response = sbbsntx) %>% # outcome of interest
  hypothesize(null = "independence") %>%
  calculate(stat = "t") #  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

The T-statistic signifies through standard error how much the sample mean deviates from the hypothesized population mean. A value of -4.138823 signifies a large evidence against the null hypothesis.

This implies that there is a strong statistical difference in agreeing with the statement “social benefits/services cost businesses too much in taxes/charges” based on urban/rural belonging.

Step 2: Simulating the null distribution:

null_distribution <- france_data2 %>%
  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.

step3: Calculating P value of the sample:

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

The P value is 0 but not 0- which shows that, according to frequentest testing, strong evidence exists to reject the null hypothesis at 5% significance level.

Task3:

Visualizing with confidence interval:

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 -0.208
##  2         2  1.00 
##  3         3  0.265
##  4         4  0.730
##  5         5 -0.589
##  6         6  0.382
##  7         7  0.147
##  8         8 -1.12 
##  9         9  0.653
## 10        10  1.05 
## # ℹ 990 more rows

The green highlight signifies the confidence interval. The red line highlights the T-statistic. As the T-statistic falls outside our distribution, under the under assumptions of frequentist hypothesis testing, we have enough evidence to reject the null hypothesis.

Task 4:

Loading,cleaning and conditioning both variables of interest. Lrscale turned left and right, with 4,5,6 (moderates) removed + NA values (77,88,99). NA Values for ccrdprs (66,77,88,99) also cleaned out.

france_data<-ess %>%
   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(4,5,6,77, 88, 99) ~ NA_character_# Values 77, 88, 99 in 'lrscale' are set as NA as well as removing 4,5,6 for moderates
    )) %>%
  # Ensuring NAs in ccrdprs are handled 
  mutate(
    ccrdprs = ifelse(ccrdprs %in% c(66, 77, 88,99), NA_integer_, ccrdprs),
  )

Specifying lrscale (coded for left and right) as explanatory variable and ccdprs as response variable.

mytab2<-france_data %>%
  filter(!is.na(ID))%>%
  mutate(ccrdprs=ifelse(ccrdprs%in% c(66,77,88,99), NA,ccrdprs))
test_stat <- mytab2%>%
  specify(explanatory = ID, # change variable name for explanatory variable
          response = ccrdprs) %>% # change variable name for outcome of interest
  hypothesize(null = "independence") %>%
  calculate(stat = "t") # replace in between quotation marks appropriate test statistic
## Warning: Removed 157485 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 "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 
## 25.57025

The null distribution:

null_distribution2 <- mytab2 %>%
  specify(explanatory = ID,
          response = ccrdprs) %>%
  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 157485 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 "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_distribution2 %>% 
  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

The T-statistic had a value of 5.482311 while the P-value was zero. A very low, almost zero, P-value provides adequate evidence to reject the null hypothesis. This shows that a strong difference exists within the mean scores of CCRDPRS (To what extent feel personal responsibility to reduce climate change) within individuals who identify as left leaning as opposed to those that are right leaning. An alternative hypothesis can thus be supported which states that a significant difference exists between the mean ccrdprs scores for the two groups.

task5:

conf_int <- null_distribution2 %>%
  get_confidence_interval(level = 0.95, type = "percentile")


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

The T-statistic (red line) falls outside the distribution. This, under the assumptions of statistic