Brazil_Emily_Homework_3
2024-03-10
# 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")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.
france_data1 <- ess %>%
filter(cntry == "FR") %>%
mutate(
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"))
france_data1 <- france_data1 %>%
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))france_clean <- france_data1 %>%
mutate(
ID = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
lrscale %in% 4:6 ~ NA_character_,
lrscale %in% c(77, 88, 99) ~ NA_character_
)) %>%
filter(!is.na(ID))newfrance <- france_clean
tab1 <- table(newfrance$educ.ba, newfrance$lrscale)
rsums <- rowSums(tab1)
csums <- colSums(tab1)
N <- sum(tab1)ptab1 <- tcrossprod(rsums/N, csums/N)
cat("Table of Expected Proportions:\n")## Table of Expected Proportions:print(ptab1)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.02631397 0.01389489 0.0353359 0.05557956 0.04870173 0.0405987 0.01138880
## [2,] 0.07855785 0.04148187 0.1054920 0.16592748 0.14539437 0.1212035 0.03400017
## [,8]
## [1,] 0.01910199
## [2,] 0.05702718ftab1 <- N * ptab1
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab1)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 237.1152 125.2069 318.4118 500.8274 438.8513 365.8349 102.6245 172.1281
## [2,] 707.8848 373.7931 950.5882 1495.1726 1310.1487 1092.1651 306.3755 513.8719alpha <- 0.05 # significance level
c_val <- qchisq(alpha, df = 1, lower.tail = FALSE)
cat("Critical Value:", round(c_val, 3), "\n")## Critical Value: 3.841test_stat1 <- sum((tab1 - ftab1)^2 / ftab1)
cat("Pearson's X^2:", round(test_stat1, 4), "\n")## Pearson's X^2: 226.1399p_val1 <- pchisq(test_stat1, df = 1, lower.tail = F)
cat("p-value:", round(p_val1, 4), "\n")## p-value: 0cat("Chi-squared test result:\n")## Chi-squared test result:print(chisq.test(tab1, correct = F))##
## Pearson's Chi-squared test
##
## data: tab1
## X-squared = 226.14, df = 7, p-value < 2.2e-16The critical value calculated of 3.841 compared to the x-squared value of 226.1399 tells us that there is a statistically significant relationship between participants placement on left-right scale and their education level. Given that the calculated p-value was 2.2e-16 and the chi-squared test statistic gave a value of 3.841, we can reject the null hypothesis as the p-value is less than 3.841 at a significance level of 0.05.
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).
france_data2 <- ess %>%
filter(cntry == "FR") %>%
mutate(
sbbsntx = case_when (
sbbsntx == 1 ~ "Agree Strongly",
sbbsntx == 2 ~ "Agree",
sbbsntx == 3 ~ "Neither Agree Nor Disagree",
sbbsntx == 4 ~ "Disagree",
sbbsntx == 5 ~ "Disagree Strongly",
sbbsntx %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(sbbsntx)
)
) %>%
filter(!is.na(sbbsntx))
unique(france_data2$sbbsntx)## [1] "Agree Strongly" "Neither Agree Nor Disagree"
## [3] "Agree" "Disagree Strongly"
## [5] "Disagree"france_data2 <- france_data2 %>%
mutate(geo = recode(as.character(domicil),
'1' = "Urban",
'2' = "Urban",
'3' = "Rural",
'4' = "Rural",
'5' = "Rural",
'7' = NA_character_,
'8' = NA_character_,
'9' = NA_character_))
france_data2 <- france_data2 %>%
filter(!is.na(geo))
unique(france_data2$geo)## [1] "Rural" "Urban"test_stat2 <- france_data2 %>%
specify(explanatory = geo,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
calculate (stat = "Chisq")
print(test_stat2$stat)## X-squared
## 19.97187null_distribution2 <- france_data2 %>%
specify(explanatory = geo,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")p_val2 <- null_distribution2 %>%
get_pvalue(obs_stat = test_stat2, 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_val2## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0Since the p-value is 0.004 and the x-squared value is 19.97187, this indicates that the null hypothesis can be rejected based on these statistical results.
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.
conf_int2 <- null_distribution2 %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution2 %>%
visualize() +
shade_p_value(obs_stat = test_stat2, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int2)## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
null_distribution2## Response: sbbsntx (factor)
## Explanatory: geo (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 2.24
## 2 2 1.25
## 3 3 5.87
## 4 4 4.42
## 5 5 4.18
## 6 6 4.15
## 7 7 2.91
## 8 8 1.69
## 9 9 0.903
## 10 10 4.98
## # ℹ 990 more rowsSince the red line representing the x-squared value falls at the end points, this means that the null hypothesis can be rejected as a test statistic inside the shaded area would mean that the null hypothesis cannot be rejected.
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).
france_data3 <- ess %>%
filter(cntry == "FR") %>%
mutate(
ccrdprs = case_when (
ccrdprs == 0 ~ "Not at all",
ccrdprs == 1 ~ "1",
ccrdprs == 2 ~ "2",
ccrdprs == 3 ~ "3",
ccrdprs == 4 ~ "4",
ccrdprs == 5 ~ "5",
ccrdprs == 6 ~ "6",
ccrdprs == 7 ~ "7",
ccrdprs == 8 ~ "8",
ccrdprs == 9 ~ "9",
ccrdprs == 10 ~ "10",
ccrdprs %in% c(66, 77, 88, 99) ~ NA_character_,
TRUE ~ as.character(ccrdprs)
)
) %>%
filter(!is.na(ccrdprs))
unique(france_data3$ccrdprs)## [1] "10" "8" "2" "7" "6"
## [6] "5" "1" "9" "Not at all" "4"
## [11] "3"france_data3 <- france_data3 %>%
mutate(
ID = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
lrscale %in% 4:6 ~ NA_character_,
lrscale %in% c(77, 88, 99) ~ NA_character_
)) %>%
filter(!is.na(ID))
unique(france_data3$ID)## [1] "Left" "Right"test_stat3 <- france_data3 %>%
specify(explanatory = ID, # change variable name for explanatory variable
response = ccrdprs) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "Chisq") # replace in between quotation marks appropriate test statistic
print(test_stat3$stat)## X-squared
## 35.04425null_distribution3 <- france_data3 %>%
specify(explanatory = ID,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")p_val3 <- null_distribution3 %>%
get_pvalue(obs_stat = test_stat3, 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_val3## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0Since the p-value is very small and close to zero and the x-squared value is 35.04425, this indicates that the null hypothesis can be rejected based on these statistical results.
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.
conf_int3 <- null_distribution3 %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution3 %>%
visualize() +
shade_p_value(obs_stat = test_stat3, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int3)## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
null_distribution3## Response: ccrdprs (factor)
## Explanatory: ID (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 8.05
## 2 2 4.45
## 3 3 11.0
## 4 4 9.53
## 5 5 12.3
## 6 6 8.53
## 7 7 1.91
## 8 8 6.13
## 9 9 8.95
## 10 10 9.73
## # ℹ 990 more rowsSince the red line representing the x-squared value falls at the end points outside of the shaded area, this means that the null hypothesis can be rejected at the confidence level of 0.05, under assumptions of frequentist hypothesis testing.