Lin_Zetong__Homework_3
2024-03-10
Setting up environment
# 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## [[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")Filtering to country, recoding and cleaning variables of interest
france_data <- ess %>%
# Filter for France
filter(cntry == "FR") %>%
# Recoding clsprty and cleaning trstplt (will use later)
mutate(
clsprty = case_when(
clsprty == 1 ~ "Yes", # Recode 1 to "Yes"
clsprty == 2 ~ "No", # Recode 2 to "No"
clsprty %in% c(7, 8, 9) ~ NA_character_, # Handle other specific cases where you want to set it as NA
TRUE ~ as.character(clsprty) # Keep other values as-is but ensure they are characters
),
trstplt = ifelse(trstplt %in% c(77, 88, 99), NA, trstplt),
# Recoding 'lrscale' into a new variable 'ID'
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" omit 'moderates for now
lrscale %in% c(77, 88, 99) ~ NA_character_ # Values 77, 88, 99 in 'lrscale' are set as NA
),
# Recoding gender
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
# Recoding education with case_when for consistency
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
)Remove NAs
# using unique to identify
unique(france_data$gndr)## [1] "Female" "Male"# does not return NA, but the following does
unique(france_data$clsprty)## [1] "Yes" "No" NA# here's how to make sure it's really cleaned
france_data <- france_data %>% filter(!is.na(clsprty))
france_data <- france_data %>% filter(!is.na(trstplt))
# double-checking + quick distribution
unique(france_data$clsprty)## [1] "Yes" "No"table(france_data$clsprty)##
## No Yes
## 9098 9480Task 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.
tab_dat <- france_data
# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$ID, 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.1320008 0.3911177
## [2,] 0.1203336 0.3565479# Table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 1173.355 3476.645
## [2,] 1069.645 3169.355# 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.841The critical value is 3.841. If the Chi-squared test statistic is greater than 3.841, then we reject the null hypothesis at significance level of 0.05. Moreover, we could conclude that Ideological self-placement is associated with education level in our dataset.
test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")## Pearson's X^2: 100.1051Since the test statistic is way greater than 3.841. We reject the null hypothesis at the 0.05 significance level. This would mean that there is a statistically significant association between Ideological self-placement and education level in the dataset.
# 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 = 100.11, df = 1, p-value < 2.2e-16The p-vale is smaller than 2.2e-16 in this quewstion and it is extremely close to 0. It represents that we have strong evidence that Ideological self-placement is associated with education level significantly.
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).、
#recoding
table(ess$sbbsntx)##
## 1 2 3 4 5 7 8 9
## 6577 27773 26302 25689 5868 122 8707 101table(ess$domicil)##
## 1 2 3 4 5 7 8 9
## 107787 55216 149256 148750 27534 137 493 1382task2 <- ess %>%
filter(cntry == "FR") %>%
mutate(
sbbsntx = case_when (
sbbsntx %in% c(1) ~ "Strong Agree ",
sbbsntx %in% c(2) ~ "Agree",
sbbsntx %in% c(3) ~ "neutral",
sbbsntx %in% c(4) ~ "Disagree",
sbbsntx %in% c(5) ~ "Strong Disagree ",
sbbsntx %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(sbbsntx)
)
) %>%
filter(!is.na(sbbsntx))
unique(task2$sbbsntx)## [1] "Strong Agree " "neutral" "Agree" "Strong Disagree "
## [5] "Disagree"task2 <- task2 %>%
mutate(exp = recode(as.character(domicil),
'1' = "Urban",
'2' = "Urban",
'3' = "Rural",
'4' = "Rural",
'5' = "Rural",
'7' = NA_character_,
'8' = NA_character_,
'9' = NA_character_))
task2 <- task2 %>%
filter(!is.na(exp))
unique(task2$exp)## [1] "Rural" "Urban"stat2 <- task2 %>%
specify(explanatory = exp,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
calculate (stat = "Chisq")
print(stat2$stat)## X-squared
## 19.97187null2 <- task2 %>%
specify(explanatory = exp,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")p2 <- null2 %>% # replace name here if you assigned something other than null_distribution above
get_pvalue(obs_stat = stat2, 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.p2## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0Since the p-value is less than 0.05, we reject the null hypothesis. Moreover, chi-sqaue value is greater than test statistic, we could conclude that whether the education level of people in urban are different from rural area.They are associated statistically significant.
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 <-null2 %>%
get_confidence_interval(level = 0.95, type = "percentile")
null2 %>%
visualize() +
shade_p_value(obs_stat = stat2, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int)## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
According to the plot,the red line(test statistic) is away from the
distribution. there is clear evidence shown that we reject the null
hypothesis.
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(ess$ccrdprs)##
## 0 1 2 3 4 5 6 7 8 9 10 66 77
## 6294 2668 4530 5706 5434 15489 12049 16523 15951 5891 9098 1939 79
## 88 99
## 1792 629france_data <- ess %>%
filter(cntry == "FR") %>%
mutate(
lrscale = case_when(
lrscale %in% 0:4 ~ "Left",
lrscale %in% 7:10 ~ "Right",
TRUE ~ NA_character_,
ccrdprs %in% c(66, 77, 88, 99) ~ NA_character_,
)
) %>%
filter(!is.na(lrscale), !is.na(ccrdprs))test_stat <- france_data %>%
specify(explanatory = lrscale,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
calculate(stat = "t") ## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "Left" - "Right", or divided in the order "Left" / "Right" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Left", "Right")` to the calculate() function.print(test_stat$stat)## t
## 0.7354004null_distribution <- france_data %>%
specify(explanatory = lrscale,
response = ccrdprs) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "t")## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "Left" - "Right", or divided in the order "Left" / "Right" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Left", "Right")` to the calculate() function.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.436Since the t-statistics is small, we don’t have sufficient evidence to reject the null hypothesis. Therefore, To what extent feel personal responsibility to reduce climate change are not significantly associated with education level. Moreover, p-val is greater than 0.05, it also indicates that we don’t have sufficient evidence to reject null hypohesis.
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_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)
Acoording to the plot, the test statistic fall within null distribution
. Therefore. our observed statistic is not extreme. we don’t reject the
hypotheiss and there is no evidence shown that to what extent feel
personal responsibility to reduce climate change are significantly
associated with education level.