# 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"
ess <- read_fst("All-ESS-Data.fst")
france_data <- read.fst("france_data.fst")
# List of packages
packages <- c("tidyverse", "modelsummary", "fst", "viridis", "knitr", "rmarkdown", "questionr", "ggridges")
# 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)
## Loading required package: viridisLite
## [[1]]
## [1] "fstcore" "fst" "infer" "lubridate" "forcats" "stringr"
## [7] "dplyr" "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [13] "tidyverse" "stats" "graphics" "grDevices" "utils" "datasets"
## [19] "methods" "base"
##
## [[2]]
## [1] "modelsummary" "fstcore" "fst" "infer" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[3]]
## [1] "modelsummary" "fstcore" "fst" "infer" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[4]]
## [1] "viridis" "viridisLite" "modelsummary" "fstcore" "fst"
## [6] "infer" "lubridate" "forcats" "stringr" "dplyr"
## [11] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [16] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [21] "datasets" "methods" "base"
##
## [[5]]
## [1] "knitr" "viridis" "viridisLite" "modelsummary" "fstcore"
## [6] "fst" "infer" "lubridate" "forcats" "stringr"
## [11] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [16] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [21] "utils" "datasets" "methods" "base"
##
## [[6]]
## [1] "rmarkdown" "knitr" "viridis" "viridisLite" "modelsummary"
## [6] "fstcore" "fst" "infer" "lubridate" "forcats"
## [11] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [16] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [21] "grDevices" "utils" "datasets" "methods" "base"
##
## [[7]]
## [1] "questionr" "rmarkdown" "knitr" "viridis" "viridisLite"
## [6] "modelsummary" "fstcore" "fst" "infer" "lubridate"
## [11] "forcats" "stringr" "dplyr" "purrr" "readr"
## [16] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [21] "graphics" "grDevices" "utils" "datasets" "methods"
## [26] "base"
##
## [[8]]
## [1] "ggridges" "questionr" "rmarkdown" "knitr" "viridis"
## [6] "viridisLite" "modelsummary" "fstcore" "fst" "infer"
## [11] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [16] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"
#ANSWER: The Critical Value for independence for France is 3.841, which means the significant level at 0.05 is rejected. In other words there is close statistic association between education level and feeling close to a party. The X squared value is much higher than the critical value at 166.18. Which means that the feeling of closeness to a party does not appear the same at all education levels. As the P-value is at 0 enforces that the null hypothesis is rejected and the two variables are independant.
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 <- ess %>%
filter(cntry == "FR") %>%
mutate(
clsprty = case_when(
clsprty == 1 ~ "Yes", # Recode 1 to "Yes"
clsprty == 2 ~ "No", # Recode 2 to "No"
clsprty %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(clsprty)
),
trstplt = ifelse(trstplt %in% c(77, 88, 99), NA, trstplt),
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"
)
) %>%
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$clsprty, 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.1101248 0.3815066
## [2,] 0.1138739 0.3944946
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")
## Table of Expected Frequencies:
print(ftab)
## [,1] [,2]
## [1,] 2059.444 7134.556
## [2,] 2129.556 7377.444
# 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
# 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: 166.1849
# 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 = 166.18, df = 1, p-value < 2.2e-16
##Task 2
##Answer: Since the p-value is less than 0.05, the null hypothesis is rejected. Which means that social benefits/services cost and domicile in Europe have an significant association with one another. Meaning that both the variables are not independant from one another.
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).
df <- france_data
unique(df$geo)
## NULL
unique(df$domicil)
## [1] 3 2 1 4 5 8 7
df <- df %>% filter(!is.na(domicil))
df <- df %>% filter(!is.na(domicil))
unique(df$domicil)
## [1] 3 2 1 4 5 8 7
table(df$domicil)
##
## 1 2 3 4 5 7 8
## 3584 2361 6257 5668 1162 3 3
df <- df %>%
mutate(geo = recode(as.character(domicil),
'1' = "Urban",
'2' = "Urban", # or set to Urban | Regardless decision needs to be justified
'3' = "Rural",
'4' = "Rural",
'5' = "Rural",
'7' = NA_character_,
'8' = NA_character_,
'9' = NA_character_))
table(df$geo)
##
## Rural Urban
## 13087 5945
table(df$domicil)
##
## 1 2 3 4 5 7 8
## 3584 2361 6257 5668 1162 3 3
unique(df$Sbsc)
## NULL
unique(df$sbbsntx)
## [1] NA 1 3 2 5 4 8 7
df <- df %>% filter(!is.na(sbbsntx))
df <- df %>% filter(!is.na(sbbsntx))
unique(df$sbbsntx)
## [1] 1 3 2 5 4 8 7
table(df$sbbsntx)
##
## 1 2 3 4 5 7 8
## 757 1538 804 640 301 3 100
df <- df %>%
mutate(Sbsc = recode(as.character(sbbsntx),
'1' = "Strongly Agree",
'2' = "Agree",
'3' = "Nor Agree or Disagree",
'4' = "Disagree",
'5' = "Strongly Disagree",
'7' = NA_character_,
'8' = NA_character_,
'9' = NA_character_))
table(df$Sbsc)
##
## Agree Disagree Nor Agree or Disagree
## 1538 640 804
## Strongly Agree Strongly Disagree
## 757 301
table(df$sbbsntx)
##
## 1 2 3 4 5 7 8
## 757 1538 804 640 301 3 100
test_stat <- df %>%
specify(explanatory = geo,
response = Sbsc) %>% hypothesize(null = "independence") %>%
calculate(stat = "Chisq")
## Warning: Removed 104 rows containing missing values.
print(test_stat$stat)
## X-squared
## 19.97187
null_distribution <- df %>%
specify(explanatory = geo,
response = Sbsc) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")
## Warning: Removed 104 rows containing missing values.
p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "greater")
## 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
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.
##Answer: Looking at the output for the null distribution, the graph shows that the x-sqaured value is to extreme and that it must reject the null hypothesis. Also since the p-value is less than x, it means that it is too unusual to the null distribution in this scenerio. Meaning it is a statistically significant result.
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "greater")
## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
null_distribution
## Response: Sbsc (factor)
## Explanatory: geo (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 11.0
## 2 2 9.99
## 3 3 2.31
## 4 4 11.7
## 5 5 8.03
## 6 6 0.983
## 7 7 5.75
## 8 8 2.68
## 9 9 2.79
## 10 10 6.54
## # ℹ 990 more rows
##Answer: p-value is less than 0.05, which means that the null hypothesis is rejected. So personal responsibilty to reduce climate change and the left right scale in Europe are in significant association with each other. The two variables are not independant from one another.
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).
df <- france_data
df <- df %>%
mutate(ltrs = recode(as.character(lrscale),
'0' = "Left",
'1' = "1",
'2' = "2",
'3' = "3",
'4' = "4",
'5' = "5",
'6' = "6",
'7' = "7",
'8' = "8",
'9' = "9",
'10' = "Right",
'77' = NA_character_,
'88' = NA_character_,
'99' = NA_character_))
# check
table(df$ltrs)
##
## 1 2 3 4 5 6 7 8 9 Left Right
## 499 1269 1996 1838 5274 1594 1749 1458 409 945 686
# compare
table(df$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
df <- df %>%
mutate(ccrdprs = ifelse(ccrdprs %in% c(66,77,88,99), NA, ccrdprs ),
'0' = "0",
'1' = "1",
'2' = "2",
'3' = "3",
'4' = "4",
'5' = "5",
'6' = "6",
'7' = "7",
'8' = "8",
'9' = "9",
'10' = "10",
'66' = NA,
'77' = NA,
'88' = NA,
'99' = NA)
# compare
table(df$ccrdprs)
##
## 0 1 2 3 4 5 6 7 8 9 10
## 72 26 64 78 114 478 403 736 952 358 708
test_stat2 <- df %>%
specify(explanatory = ltrs,
response = ccrdprs) %>% hypothesize(null = "independence") %>%
calculate(stat = "F")
## Warning: Removed 15413 rows containing missing values.
print(test_stat2$stat)
## [1] 3.7871
null_distributions <- df %>%
specify(explanatory = ltrs,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "F")
## Warning: Removed 15413 rows containing missing values.
p_vals <- null_distributions %>%
get_pvalue(obs_stat = test_stat2, direction = "greater")
## 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_vals
## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0
##Answer: Null distribution on the graph indicates that the x-sqaured value is extreme and the null hypothesis is rejected. The p-value is less than x, it is too unusual to the null distribution, it is a statistically significant result.
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_distributions %>%
visualize() +
shade_p_value(obs_stat = test_stat2, direction = "greater")
## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
null_distributions
## Response: ccrdprs (numeric)
## Explanatory: ltrs (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 1.14
## 2 2 0.788
## 3 3 0.582
## 4 4 0.441
## 5 5 1.34
## 6 6 0.970
## 7 7 1.48
## 8 8 1.20
## 9 9 0.928
## 10 10 0.759
## # ℹ 990 more rows