# 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")##task 1
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# does not return NA, but the following does
unique(france_data$educ.ba)## [1] "No BA" "BA or more"# 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-16since the p-value is lower it shows that there is no association between left or right win and having bachelors degree or not. The critical value is 3.841 which indicates that the chi-squared statistic is above the critical value
##Task2
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.06since the p value is small this shows that the probbility is lardger that domicil and sbbsntx have great statistical difference. the depenedent variable is related to the independent variable and this is because the chi square value is above the critical value
##task3
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 20.3
## 2 2 22.5
## 3 3 24.3
## 4 4 26.2
## 5 5 27.6
## 6 6 15.1
## 7 7 28.8
## 8 8 25.1
## 9 9 38.7
## 10 10 23.6
## # ℹ 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 20.3
## 2 2 22.5
## 3 3 24.3
## 4 4 26.2
## 5 5 27.6
## 6 6 15.1
## 7 7 28.8
## 8 8 25.1
## 9 9 38.7
## 10 10 23.6
## # ℹ 990 more rowsanything shaded in red means distribution is more extreme then the observed test statistics
##Task4
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")
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.002p-value = 0 shows a great statistical relation between left and right wing and feeling a responsility to help with climate change. the varaibles can be realted, the chi squared value is above the critical value which allows enough the reject the null hypothesis
##task5
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 171.
## 2 2 197.
## 3 3 133.
## 4 4 137.
## 5 5 170.
## 6 6 152.
## 7 7 166.
## 8 8 183.
## 9 9 167.
## 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 171.
## 2 2 197.
## 3 3 133.
## 4 4 137.
## 5 5 170.
## 6 6 152.
## 7 7 166.
## 8 8 183.
## 9 9 167.
## 10 10 192.
## # ℹ 990 more rowsnull hypothesis is rejected under frequentest hypothesis testing because the test statistic is outside the distribution