Graves-Rice_Olivia_Homework_3
Olivia Graves-Rice
2024-03-08
Homework 3 (5%): due by next lecture on March 12th
# 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"library(fst)
france_data <- read_fst("france_data.fst")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_data <- ess %>%
# Filter for France
filter(cntry == "FR") %>%
mutate(
lrscale = case_when(
lrscale == 0:3 ~ "Left",
lrscale == 7:10 ~ "Right",
lrscale %in% c(7, 8, 9) ~ NA_character_, # Handle other specific cases where you want to set it as NA
TRUE ~ as.character(lrscale)
),
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
# Ensure all NAs are handled uniformly
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) # Force educ.ba to be treated as character
)## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `lrscale = case_when(...)`.
## Caused by warning in `lrscale == 0:3`:
## ! longer object length is not a multiple of shorter object length
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.tab_dat <- france_data
mytab <- table(tab_dat$lrscale, 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.009915507 0.03424447
## [2,] 0.005074649 0.01752595
## [3,] 0.007027495 0.02427035
## [4,] 0.013202340 0.04559597
## [5,] 0.020559894 0.07100622
## [6,] 0.025276980 0.08729728
## [7,] 0.072530355 0.25049286
## [8,] 0.021921385 0.07570831
## [9,] 0.005239679 0.01809590
## [10,] 0.012927291 0.04464605
## [11,] 0.016007837 0.05528511
## [12,] 0.014852632 0.05129547# Table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 161.89049 559.1095
## [2,] 82.85380 286.1462
## [3,] 114.73792 396.2621
## [4,] 215.55460 744.4454
## [5,] 335.68139 1159.3186
## [6,] 412.69725 1425.3028
## [7,] 1184.20310 4089.7969
## [8,] 357.91046 1236.0895
## [9,] 85.54823 295.4518
## [10,] 211.06388 728.9361
## [11,] 261.35996 902.6400
## [12,] 242.49893 837.5011# 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: 439.5129# 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 = 439.51, df = 11, p-value < 2.2e-16#Upon reviewing the data calculated, we fail to reject the null because the p-value is 2.2 meaning greater than 0.05. This would mean there is no significant association between the two variables. The Pearson’s x2 is also greater than the critical value since the x-squared is 439.51 and the critical value is 3.841.
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_data <- france_data %>%
mutate(
geo = recode(as.character(sbbsntx),
'1' = "Agree strongly",
'2' = "Agree",
'3' = "Neither agree nor disagree",
'4' = "Disagree",
'5' = "Disagree strongly",
'7' = NA_character_,
'8' = NA_character_,
'9' = NA_character_),
sbbsntx = ifelse(sbbsntx %in% c(7, 8, 9), NA, sbbsntx),
sbbsntx = ifelse(sbbsntx == 2, 0, sbbsntx),
) %>%
filter(!is.na(sbbsntx), !is.na(geo))
cond <- france_data %>%
drop_na(sbbsntx) %>%
count(sbbsntx, geo) %>%
group_by(geo) %>%
mutate(prob = n / sum(n))
cond## # A tibble: 5 × 4
## # Groups: geo [5]
## sbbsntx geo n prob
## <dbl> <chr> <int> <dbl>
## 1 0 Agree 1538 1
## 2 1 Agree strongly 757 1
## 3 3 Neither agree nor disagree 804 1
## 4 4 Disagree 640 1
## 5 5 Disagree strongly 301 1france_data <- france_data %>%
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_),
domicil = ifelse(domicil %in% c(7, 8, 9), NA, domicil),
domicil = ifelse(domicil == 2, 0, domicil)
)
cond <- france_data %>%
drop_na(domicil) %>%
count(domicil, geo) %>%
group_by(geo) %>%
mutate(prob = n / sum(n))
print(cond)## # A tibble: 5 × 4
## # Groups: geo [2]
## domicil geo n prob
## <dbl> <chr> <int> <dbl>
## 1 0 Urban 479 0.412
## 2 1 Urban 684 0.588
## 3 3 Rural 1416 0.492
## 4 4 Rural 1221 0.425
## 5 5 Rural 239 0.0831test_stat <- france_data %>%
specify(explanatory = domicil, # change variable name for explanatory variable
response = sbbsntx) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "slope") # replace in between quotation marks appropriate test statistic## Warning: Removed 1 rows containing missing values.## Message: The independence null hypothesis does not inform calculation of the
## observed statistic (a slope) and will be ignored.print(test_stat$stat)## domicil
## -0.07032579null_distribution <- france_data %>%
specify(explanatory = domicil,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>% # only line we are adding, use reps = 1000 as standard but change if too computationally demanding
calculate(stat = "slope")## Warning: Removed 1 rows containing missing values.p_val <- null_distribution %>% # replace name here if you assigned something other than null_distribution above
get_pvalue(obs_stat = test_stat, direction = "two-sided") # would only replace test_stat if assigned another name in Step 1## 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#The p-value for this varibale is 0 meaning they are more to the left when believing social benefits/services cost businesses too much in taxes/charges meaning they most likely strongly agree with the variable.
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_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)
null_distribution## Response: sbbsntx (numeric)
## Explanatory: domicil (numeric)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.00632
## 2 2 -0.0277
## 3 3 0.00747
## 4 4 0.0151
## 5 5 0.0128
## 6 6 -0.00464
## 7 7 -0.0289
## 8 8 -0.0164
## 9 9 -0.0226
## 10 10 -0.00554
## # ℹ 990 more rows#The null distribution present ranges from -0.06 to 0.06 and is also to the right of the median. This graph also shows that the null distribution does not include values -0.035 to 0.035 since they do not fall in the net that is provided. With this the graph is centered around 0 whih makes sense due to teh pvalue being 0. the data is also consitantly around 125 - 150 mark.
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_data <- ess %>%
filter(cntry == "FR") %>%
mutate(
lrscale = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
lrscale %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(lrscale)
)
)test_stat <- france_data %>%
specify(explanatory = lrscale, # change variable name for explanatory variable
response = ccrdprs) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "F") # replace in between quotation marks appropriate test statistic## Warning: Removed 14991 rows containing missing values.print(test_stat$stat)## [1] 1.54944null_distribution <- france_data %>%
specify(explanatory = lrscale,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>% # only line we are adding, use reps = 1000 as standard but change if too computationally demanding
calculate(stat = "F")## Warning: Removed 14991 rows containing missing values.p_val <- null_distribution %>% # replace name here if you assigned something other than null_distribution above
get_pvalue(obs_stat = test_stat, direction = "two-sided") # would only replace test_stat if assigned another name in Step 1
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.272#The data present shows that France is leaning to the left when it comes to believeing it is a persons responsibility to reduce climate change due to the data showing the p-value being 0.25. Since the value is under 1 it would seem that they believe that poeple should nothave to be responsible at all.
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: F usually corresponds to right-tailed tests. Proceed with caution.
null_distribution## Response: ccrdprs (numeric)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.770
## 2 2 0.754
## 3 3 0.415
## 4 4 0.628
## 5 5 1.57
## 6 6 1.10
## 7 7 2.93
## 8 8 0.373
## 9 9 0.774
## 10 10 1.21
## # ℹ 990 more rows#Based on the visualization it seems that the null distribution is to the left of the median since the p-value is 0.25, with the highest value being 300 and the lowest being around 20. This also depicts a two tailed test.