# 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("france_data.fst")task1
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 ~ NA_character_, # Values 4 to 6 in 'lrscale' are categorized as "Moderate"
lrscale %in% c(77, 88, 99) ~ NA_character_
)
)
table(france_data$ID)##
## Left Right
## 4709 4302france_data <- ess %>%
mutate(
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
),
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
)
educ_table <- table(france_data$educ.ba)
educ_table##
## BA or more No BA
## 4252 14786# using unique to identify
unique(france_data$lrscale)## [1] 3 2 1 4 5 6 10 8 88 7 9 0 77# here's how to make sure it's really cleaned
france_data <- france_data %>% filter(!is.na(lrscale))
france_data <- france_data %>% filter(!is.na(educ.ba))
# double-checking + quick distribution
unique(france_data$lrscale)## [1] 3 2 1 4 5 6 10 8 88 7 9 0 77table(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 940tab_dat <- france_data
# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$lrscale, 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)
# 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.011086193 0.03855137
## [2,] 0.005853979 0.02035676
## [3,] 0.014887173 0.05176899
## [4,] 0.023415916 0.08142703
## [5,] 0.021562351 0.07498140
## [6,] 0.061871513 0.21515338
## [7,] 0.018699885 0.06502740
## [8,] 0.020518255 0.07135064
## [9,] 0.017104411 0.05947926
## [10,] 0.004798151 0.01668520
## [11,] 0.008047755 0.02798544
## [12,] 0.004469671 0.01554294
## [13,] 0.011027535 0.03834740# Table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 211.0589 733.9411
## [2,] 111.4481 387.5519
## [3,] 283.4220 985.5780
## [4,] 445.7922 1550.2078
## [5,] 410.5040 1427.4960
## [6,] 1177.9099 4096.0901
## [7,] 356.0084 1237.9916
## [8,] 390.6265 1358.3735
## [9,] 325.6338 1132.3662
## [10,] 91.3472 317.6528
## [11,] 153.2132 532.7868
## [12,] 85.0936 295.9064
## [13,] 209.9422 730.0578# 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: 510.8805# 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 = 510.88, df = 12, p-value < 2.2e-16##Based on the data, it can be extrapolated that a left or right wing is correlated with obtaining a bachelor’s degree (p-value). In this case, the Chi-squared value is 510.8805, which shows a strong connection between the two variables.
Task 2:
table_domicil <- table(france_data$domicil)
print(table_domicil)##
## 1 2 3 4 5 7 8
## 3584 2361 6257 5668 1162 3 3table_sbbsntx <- table(france_data$sbbsntx)
print(table_sbbsntx)##
## 1 2 3 4 5 7 8
## 757 1538 804 640 301 3 100france_clean <- france_data %>%
mutate(
ID = case_when(
sbbsntx %in% 1 ~ "Agree Strongly",
sbbsntx %in% 2 ~ "Agree",
sbbsntx %in% 3 ~ "Neither agree nor disagree",
sbbsntx %in% 4 ~ "Disagree",
sbbsntx %in% 5 ~ "Disagree Strongly",
sbbsntx %in% c(77, 88, 99) ~ NA_character_
)
)
table(france_clean$ID)##
## Agree Agree Strongly
## 1538 757
## Disagree Disagree Strongly
## 640 301
## Neither agree nor disagree
## 804france_clean <- france_data %>%
mutate(
ID = case_when(
domicil %in% 1:2 ~ "Urban",
domicil %in% 3:5 ~ "Rural",
domicil %in% c(77, 88, 99) ~ NA_character_
)
)
table(france_data$ID)## < table of extent 0 >library(infer)
library(dplyr)
france_data <- france_data %>%
filter(!is.na(domicil), !is.na(sbbsntx))
france_data <- france_data %>%
mutate(
domicil = factor(domicil),
sbbsntx = factor(sbbsntx)
)
test_stat <- france_data %>%
specify(explanatory = domicil,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
calculate(stat = "chisq")
print(test_stat$stat)## X-squared
## 49.95025null_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.086##Checking the data, a low p-value suggests a link between domicile and sbbs tx. If Chi-squared values go over the critical point, it means there’s a connection between the dependent and independent variables.
Task 3:
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 57.0
## 2 2 21.0
## 3 3 21.9
## 4 4 30.0
## 5 5 26.1
## 6 6 24.4
## 7 7 17.2
## 8 8 32.4
## 9 9 17.6
## 10 10 24.7
## # ℹ 990 more rowsconf_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 57.0
## 2 2 21.0
## 3 3 21.9
## 4 4 30.0
## 5 5 26.1
## 6 6 24.4
## 7 7 17.2
## 8 8 32.4
## 9 9 17.6
## 10 10 24.7
## # ℹ 990 more rows##The statistic that was observed is contained in the distribution. Every red-shaded region represents a distribution that is more severe.
Task 4:
table_ccrdprs <- table(france_data$ccrdprs)
print(table_ccrdprs)##
## 0 1 2 3 4 5 6 7 8 9 10 66 77 88 99
## 48 18 49 44 73 273 210 383 476 167 287 24 1 16 1table_lrscale <- table(france_data$lrscale)
print(table_lrscale)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88
## 200 111 294 425 374 1101 364 411 356 96 133 98 180france_clean <- france_clean %>%
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 ~ NA_character_, # Values 4 to 6 in 'lrscale' are categorized as "Moderate"
lrscale %in% c(77, 88, 99) ~ NA_character_
)
)
table(france_clean$ID)##
## Left Right
## 4709 4302france_clean <- france_clean %>%
mutate(
ccrdprs = ifelse(ccrdprs %in% c(66, 77, 88, 99), NA, ccrdprs),
)
france_clean <- france_clean %>% filter(!is.na(ccrdprs))
table(france_clean$ccrdprs)##
## 0 1 2 3 4 5 6 7 8 9 10
## 72 26 64 78 114 478 403 736 952 358 708library(infer)
library(dplyr)
france_data <- france_data %>%
filter(!is.na(ccrdprs), !is.na(lrscale))
france_data <- france_data %>%
mutate(
lrscale = factor(lrscale),
ccrdprs = factor(ccrdprs)
)
test_stat <- france_data %>%
specify(explanatory = lrscale,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
calculate(stat = "chisq")
print(test_stat$stat)## X-squared
## 266.9943null_distribution <- france_data %>%
specify(explanatory = lrscale,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "chisq")
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##climate change is statistically correlated with left/right wing with a p-value of 0. If the chi-squared value is greater than the critical threshold variables can be related, which implies there is sufficient evidence to reject the null hypothesis.
Task 5:
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: ccrdprs (factor)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 162.
## 2 2 162.
## 3 3 158.
## 4 4 166.
## 5 5 153.
## 6 6 180.
## 7 7 182.
## 8 8 177.
## 9 9 148.
## 10 10 192.
## # ℹ 990 more rowsconf_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: ccrdprs (factor)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 162.
## 2 2 162.
## 3 3 158.
## 4 4 166.
## 5 5 153.
## 6 6 180.
## 7 7 182.
## 8 8 177.
## 9 9 148.
## 10 10 192.
## # ℹ 990 more rows##When testing hypotheses the null hypothesis is rejected when the observed test is outside of the distribution.