Homework 3
2024-03-01
# 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"france_data <- read.fst("france_data.fst")france_data <- france_data %>%
# 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 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
)#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.
Long form coding for Chi-Square
# Renaming dataset for clarity
tab_dat <- france_data
# Create a crosstable and save it into the object mytab
mytab <- table(tab_dat$clsprty, 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.1101248 0.3815066
## [2,] 0.1138739 0.3944946# Table of expected frequencies
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 2059.444 7134.556
## [2,] 2129.556 7377.444# 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# Calculating 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: 166.1849# Calculating the p-value
p_val <- pchisq(test_stat, df = 1, lower.tail = FALSE)
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 = FALSE))##
## Pearson's Chi-squared test
##
## data: mytab
## X-squared = 166.18, df = 1, p-value < 2.2e-16We can better understand what we would anticipate in the event that there was no association between the variables thanks to the table of expected proportions and frequencies. 3.841 is the critical value for independence at the 0.05 significance level. We reject the null hypothesis if the computed chi-squared statistic is greater than this crucial value. With a Pearson’s chi-squared value of 160.4057, the observed and expected frequencies differ significantly from one another. As the p-value is nearly zero, there is substantial evidence to refute the null hypothesis. As a result, we reject the null hypothesis that there is no correlation between education level and feeling close to a party.
#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).
# Step 1: Create a contingency table
contingency_table <- table(france_data$sbbsntx, france_data$domicil)
#Perform chi-squared test
chi_squared_test <- chisq.test(contingency_table)## Warning in chisq.test(contingency_table): Chi-squared approximation may be
## incorrectprint(chi_squared_test)##
## Pearson's Chi-squared test
##
## data: contingency_table
## X-squared = NaN, df = 36, p-value = NA# Convert 'sbbsntx' and 'domicil' to factors
france_data$sbbsntx <- factor(france_data$sbbsntx)
france_data$domicil <- factor(france_data$domicil)
# Step 2: Simulate the null distribution
null_distribution <- france_data %>%
specify(explanatory = domicil,
response = sbbsntx) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")## Warning: Removed 14895 rows containing missing values.# Step 3: Calculate the p-value of your sample
p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, 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_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0#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.
# Step 4: Visualize with Confidence Intervals
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)## 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 13.9
## 2 2 27.3
## 3 3 13.7
## 4 4 16.0
## 5 5 26.4
## 6 6 29.8
## 7 7 18.8
## 8 8 22.6
## 9 9 28.7
## 10 10 36.5
## # ℹ 990 more rowsIn Task 2, the output will include the p-value, which indicates the significance level of the relationship between the explanatory variable (domicil) and the response variable (sbbsntx), as well as the test statistic (z-statistic).
We can better comprehend the distribution of the null hypothesis in relation to the observed statistic for Task 3 thanks to the visualization featuring confidence intervals. The observed statistic provides evidence to reject the null hypothesis if it is outside the shaded regions. On the other hand, it suggests that the observed statistic is consistent with the null hypothesis if it falls inside the shaded regions.
#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).
# Step 1: Calculate the test statistic of your sample (Correlation Test)
cor_test_result <- france_data %>%
specify(explanatory = lrscale, response = ccrdprs) %>%
hypothesize(null = "independence") %>%
calculate(stat = "correlation")## Warning: Removed 14991 rows containing missing values.## Message: The independence null hypothesis does not inform calculation of the
## observed statistic (a correlation) and will be ignored.print(cor_test_result$stat)## [1] 0.03615757# Step 2: Simulate the null distribution
null_distribution <- france_data %>%
specify(explanatory = lrscale, response = ccrdprs) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "correlation")## Warning: Removed 14991 rows containing missing values.# Step 3: Calculate the p-value of your sample
p_val <- null_distribution %>%
get_pvalue(obs_stat = cor_test_result, direction = "two-sided")
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.032We can evaluate the significance of the correlation between the variables ccrdprs and lrscale using the p-value that was obtained from the test. We reject the null hypothesis of independence if the p-value is small enough (usually below a threshold like 0.05), indicating that the two variables have a significant correlation. On the other hand, if the p-value is high, the null hypothesis cannot be rejected, indicating that there is not enough data to draw the conclusion that the variables are correlated.
#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.
# Step 4: Visualize the null distribution and the observed statistic
null_distribution %>%
visualize() +
shade_p_value(obs_stat = cor_test_result, direction = "both")
The visual aid facilitates the interpretation of test results by giving
a clear picture of how the observed correlation compares to the null
distribution. There is stronger evidence against the null hypothesis and
a strong correlation between the variables if the observed correlation
is located in the extreme tails of the null distribution. On the other
hand, the evidence opposing the null hypothesis is weakened if the
observed correlation is found to be within the majority of the null
distribution, implying that the correlation may have happened by
accident.