Wilson_Marina_Homework_3
2024-03-09
# 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## 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"Homework 3:
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 <- read.fst("france_data.fst") %>%
mutate(
# Recoding 'lrscale' into a new variable 'ID'
lrscale = 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, # Values 4 to 6 in 'lrscale' are categorized as "Moderate"
lrscale %in% c(77, 88, 99) ~ NA_character_ # Values 77, 88, 99 in 'lrscale' are set as NA
),
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
)# using unique to identify
unique(france_data$educ.ba)## [1] "No BA" "BA or more"# does not return NA, but the following does
unique(france_data$lrscale)## [1] "Left" NA "Right"# 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] "Left" "Right"table(france_data$lrscale)##
## Left Right
## 4709 4302## Not necessary, but always good to rename when doing operations so you can backtrack to clean dataset
tab_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.1311243 0.3914592
## [2,] 0.1197912 0.3576253# Table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 1181.561 3527.439
## [2,] 1079.439 3222.561# 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##Interpretation: Our chi-square critical value of 3.841 is exceeded by our chi-square test statistic of 100.8551 calculated below. This indicates that we should reject the null hypothesis at the 0.05 significance level, as there is a statistically significant association between lrscale (left/right political affiliation) and educ.ba (education). # 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: 100.8551# 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.86, df = 1, p-value < 2.2e-16## Interpretation: A p-value of < 2.2e-16 means that is extremely small (essentially close to zero). Thus, we reject the null hypothesis that lrscale (left/right political affiliation) and educ.ba (education) are independent of each other. Likewise, a large chi-squared statistic of 100.86 reinforces this. If the two variables were independent, we would expect a chi-square statistic closer to zero. 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 <- read.fst("france_data.fst") %>%
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_, # Handle other specific cases where you want to set it as NA
TRUE ~ as.character(sbbsntx) # Keep other values as-is but ensure they are characters
)
)
france_data <- read.fst("france_data.fst") %>%
mutate(
domicil = case_when(
domicil == 1 ~ "Urban",
domicil == 2 ~ "Urban",
domicil == 3 ~ "Urban",
domicil == 4 ~ "Rural",
domicil == 5 ~ "Rural",
domicil %in% c(7, 8, 9) ~ NA_character_, # Handle other specific cases where you want to set it as NA
TRUE ~ as.character(domicil) # Keep other values as-is but ensure they are characters
)
) test_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 = "F") # replace in between quotation marks appropriate test statistic## Warning: Removed 14896 rows containing missing values.print(test_stat$stat)## [1] 14.6259null_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 = "F")## Warning: Removed 14896 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 = "greater") # would only replace test_stat if assigned another name in Step 1
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.001## Interpretation: In this case, a p-value of zero does not actually mean "zero"; it means that the p-value is relatively small (close to zero). The lower the p-value, the greater the statistical significance of the observed difference. Because our p-value is essentially close to zero, aka less than the threshold of 0.05, we reject the null hypothesis that there is no difference among group means between the variable domicil (urban/rural home residence location) and the variable sbbsntx (perception of whether social benefits & services cost businesses too much in taxes/charges).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 = "greater")
null_distribution## Response: sbbsntx (numeric)
## Explanatory: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.570
## 2 2 1.79
## 3 3 1.24
## 4 4 0.0370
## 5 5 1.24
## 6 6 0.214
## 7 7 0.442
## 8 8 0.106
## 9 9 5.43
## 10 10 0.0563
## # ℹ 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 = "greater") +
shade_confidence_interval(endpoints = conf_int)
null_distribution## Response: sbbsntx (numeric)
## Explanatory: domicil (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.570
## 2 2 1.79
## 3 3 1.24
## 4 4 0.0370
## 5 5 1.24
## 6 6 0.214
## 7 7 0.442
## 8 8 0.106
## 9 9 5.43
## 10 10 0.0563
## # ℹ 990 more rows## Interpretation: The null distribution is visualized to the left of the observed test statistic (red line), which indicates (under assumptions of frequentist hypothesis testing) to reject the null hypothesis. There are no red-shaded regions of the null distribution which are equally/more extreme than the observed statistic; meaning the observed statistic does not fall within the null distribution. In addition, the 95% confidence interval (from 0-5) contains the vast majority of the null distribution data, showing that the rejection of the null hypothesis is valid. 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 <- read.fst("france_data.fst") %>%
mutate(
lrscale = 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, # Values 4 to 6 in 'lrscale' are categorized as "Moderate"
lrscale %in% c(77, 88, 99) ~ NA_character_ # Values 77, 88, 99 in 'lrscale' are set as NA
)
)
france_data <- read.fst("france_data.fst") %>%
mutate(
ccrdprs = ifelse(ccrdprs %in% c(66, 77, 88, 99), NA, ccrdprs),
)testt_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 = "slope") # replace in between quotation marks appropriate test statistic## Warning: Removed 15049 rows containing missing values.## Message: The independence null hypothesis does not inform calculation of the
## observed statistic (a slope) and will be ignored.print(testt_stat$stat)## lrscale
## -0.008904258null_dis <- 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 = "slope")## Warning: Removed 15049 rows containing missing values.p_val <- null_dis %>% # replace name here if you assigned something other than null_distribution above
get_pvalue(obs_stat = testt_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## Interpretation: A p-value of 0 (or close to zero) falls below the threshold of p = 0.05, so therefore we do reject the null hypothesis that lrscale (left/right political affiliation) shares no relationship with ccrdprs (the extent to which one feels personally responsible to reduce climate change). 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_dis %>%
visualize() +
shade_p_value(obs_stat = testt_stat, direction = "two-sided")
null_dis## Response: ccrdprs (numeric)
## Explanatory: lrscale (numeric)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.000829
## 2 2 -0.00339
## 3 3 0.00131
## 4 4 0.000410
## 5 5 0.000491
## 6 6 -0.000486
## 7 7 -0.000101
## 8 8 -0.000217
## 9 9 0.0000581
## 10 10 -0.00122
## # ℹ 990 more rowsconf_int <- null_dis %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_dis %>%
visualize() +
shade_p_value(obs_stat = testt_stat, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int)
null_dis## Response: ccrdprs (numeric)
## Explanatory: lrscale (numeric)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.000829
## 2 2 -0.00339
## 3 3 0.00131
## 4 4 0.000410
## 5 5 0.000491
## 6 6 -0.000486
## 7 7 -0.000101
## 8 8 -0.000217
## 9 9 0.0000581
## 10 10 -0.00122
## # ℹ 990 more rows## Interpretation: The null distribution is visualized to the right of the observed test statistic (red line),which indicates (under assumptions of frequentist hypothesis testing) that we should not reject the null hypothesis. Our confidence interval shows that the majority of the data falls beyond the observed statistic (lining up with the peak of the null distribution), further confirming that the placement of the null distribution is correct and the null hypothesis should not be rejected. This contradicts the inference from our p-value. (I originally got 0.32 as the p-value for this question, but when I knit the file it switched to 0? Not sure what happened there, so this may be incorrect).