Salman_Eesha_Homewrk_3
2024-03-12
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.
# 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")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-16The p-value being so small indicates that there is indeed an association between being left or right-wing and having a bachelors degree or not. The critical value is at 3.841 which means that if the Chi-squared statistic is above the critical value then there is strong evidence of association between the two variables. In this case the Chi-squared result is 510.8805 which illustrates a relation of the two variables.
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).
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.068The p value being so small indicates a larger probability that domicil and sbbsntx have a great statistical difference. The Chi-square value being above the critical value indicates that the dependent variable is related to the independent.
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.
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 62.0
## 2 2 28.5
## 3 3 33.2
## 4 4 27.4
## 5 5 22.7
## 6 6 29.4
## 7 7 40.3
## 8 8 37.8
## 9 9 36.4
## 10 10 25.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 62.0
## 2 2 28.5
## 3 3 33.2
## 4 4 27.4
## 5 5 22.7
## 6 6 29.4
## 7 7 40.3
## 8 8 37.8
## 9 9 36.4
## 10 10 25.6
## # ℹ 990 more rowsThe obserevd test statistice is within the distribution. This means that everything else shaded in red are distributions more extreme than the observed test statistics.
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).
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.002The p-value being 0 indicates a great statistical relation between being left/right wing and feeling a responsibility to help with climate change. The chi-squared value being above the critical value shows that there is enough to reject the null hypothesis; the two variables can be related.
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.
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 174.
## 2 2 152.
## 3 3 124.
## 4 4 154.
## 5 5 178.
## 6 6 163.
## 7 7 141.
## 8 8 156.
## 9 9 132.
## 10 10 225.
## # ℹ 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 174.
## 2 2 152.
## 3 3 124.
## 4 4 154.
## 5 5 178.
## 6 6 163.
## 7 7 141.
## 8 8 156.
## 9 9 132.
## 10 10 225.
## # ℹ 990 more rowsThe observes test statistice is outside the distribution which shows that we reject the null hypothesis under frequentist hypothesis testing.