SOC222 research report
Jio Sohn
2024-04-11
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "interactions") # add any you need here
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_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## Warning: package 'infer' was built under R version 4.3.3## Warning: package 'effects' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Warning: package 'survey' was built under R version 4.3.3## Loading required package: grid
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Loading required package: survival
##
## Attaching package: 'survey'
##
## The following object is masked from 'package:graphics':
##
## dotchart## Warning: package 'interactions' 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"
##
## [[4]]
## [1] "modelsummary" "fst" "infer" "lubridate" "forcats"
## [6] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [11] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "effects" "carData" "modelsummary" "fst" "infer"
## [6] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [11] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [16] "stats" "graphics" "grDevices" "utils" "datasets"
## [21] "methods" "base"
##
## [[6]]
## [1] "survey" "survival" "Matrix" "grid" "effects"
## [6] "carData" "modelsummary" "fst" "infer" "lubridate"
## [11] "forcats" "stringr" "dplyr" "purrr" "readr"
## [16] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [21] "graphics" "grDevices" "utils" "datasets" "methods"
## [26] "base"
##
## [[7]]
## [1] "interactions" "survey" "survival" "Matrix" "grid"
## [6] "effects" "carData" "modelsummary" "fst" "infer"
## [11] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [16] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"sweden_data <- read.fst("sweden_data.fst")sweden_clean <- sweden_data %>%
mutate(
cptppola = ifelse(cptppola %in% c(7, 8, 9), NA, cptppola),
hinctnta = ifelse(hinctnta %in% c(77, 88, 99), NA, hinctnta),
gndr = ifelse(gndr %in% c(9), NA, gndr),
)#remove NA
sweden_clean <- sweden_clean %>% filter(!is.na(cptppola), !is.na(hinctnta), !is.na(gndr))table(sweden_clean$hinctnta)##
## 1 2 3 4 5 6 7 8 9 10
## 272 333 394 406 479 614 730 653 577 615#descriptive statistics table
table_summary <- datasummary_skim(sweden_clean %>% select(cptppola, hinctnta, gndr))
table_summary| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max | ||
|---|---|---|---|---|---|---|---|---|
| cptppola | 5 | 0 | 2.5 | 1.1 | 1.0 | 2.0 | 5.0 | |
| hinctnta | 10 | 0 | 6.2 | 2.7 | 1.0 | 7.0 | 10.0 | |
| gndr | 2 | 0 | 1.5 | 0.5 | 1.0 | 2.0 | 2.0 |
#descriptive statistics table, better aesthetics
table1 <- datasummary_skim(sweden_clean %>% dplyr::select(cptppola, hinctnta, gndr), title = "Table 1. Descriptive statistics", output = "flextable")## Warning: The histogram argument is only supported for (a) output types "default",
## "html", "kableExtra", or "gt"; (b) writing to file paths with extensions
## ".html", ".jpg", or ".png"; and (c) Rmarkdown, knitr or Quarto documents
## compiled to PDF (via kableExtra) or HTML (via kableExtra or gt). Use
## `histogram=FALSE` to silence this warning.table1
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
cptppola | 5 | 0 | 2.5 | 1.1 | 1.0 | 2.0 | 5.0 |
hinctnta | 10 | 0 | 6.2 | 2.7 | 1.0 | 7.0 | 10.0 |
gndr | 2 | 0 | 1.5 | 0.5 | 1.0 | 2.0 | 2.0 |
sweden_data <- sweden_data %>%
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
TRUE ~ NA_character_ # Set anything that is not 1 or 2 to NA
),
hinctnta = case_when(
hinctnta %in% 1:3 ~ "LowIncome", #1st to 3rd decile categorized as low income
hinctnta %in% 4:7 ~ "MidIncome", #4th to 7th decile categorized as mid income
hinctnta %in% 8:10 ~ "HighIncome", #8th to 10th decile categorized as high income
hinctnta %in% c(77, 88, 99) ~ NA_character_,
TRUE ~ NA_character_
),
cptppola = case_when(
cptppola == 1 ~ "not at all confident",
cptppola == 2 ~ "a little confident",
cptppola == 3 ~ "quite confident",
cptppola == 4 ~ "very confident",
cptppola == 5 ~ "completely confident",
TRUE ~ NA_character_
)
)#income and confidence data for graph
income_conf <- sweden_data %>%
filter(!is.na(hinctnta), !is.na(cptppola)) %>%
group_by(hinctnta, cptppola) %>%
summarise(count = n(), .groups = 'drop') %>%
mutate(percentage = count / sum(count) * 100)
income_conf## # A tibble: 15 × 4
## hinctnta cptppola count percentage
## <chr> <chr> <int> <dbl>
## 1 HighIncome a little confident 514 10.1
## 2 HighIncome completely confident 92 1.81
## 3 HighIncome not at all confident 236 4.65
## 4 HighIncome quite confident 694 13.7
## 5 HighIncome very confident 311 6.12
## 6 LowIncome a little confident 279 5.49
## 7 LowIncome completely confident 33 0.650
## 8 LowIncome not at all confident 342 6.73
## 9 LowIncome quite confident 267 5.26
## 10 LowIncome very confident 79 1.56
## 11 MidIncome a little confident 690 13.6
## 12 MidIncome completely confident 69 1.36
## 13 MidIncome not at all confident 552 10.9
## 14 MidIncome quite confident 687 13.5
## 15 MidIncome very confident 234 4.61#ggplot for personal better understanding of data
trialplot <- ggplot(income_conf, aes(x = cptppola, y = percentage, fill = cptppola)) +
geom_bar(stat = "identity", position = position_dodge()) + # Dodged bar chart
facet_wrap(~ hinctnta, scales = "fixed") +
scale_fill_brewer(palette = "Set3") +
labs(
title = "Sweden",
x = "Political Orientation",
y = "Percentage of Respondents",
fill = "Orientation"
) +
theme_minimal() +
theme(legend.position = "bottom") +
theme(axis.text.x = element_text(angle = 60, size=10, hjust = 1))
trialplot
sweden_data <- read.fst("sweden_data.fst")sweden_data$cptppola <- ifelse(sweden_data$cptppola > 5, NA, sweden_data$cptppola)
sweden_data$cptppola <- ifelse(sweden_data$cptppola > 5, NA, sweden_data$cptppola)
sweden_data$cptppola <- ifelse(sweden_data$cptppola > 5, NA, sweden_data$cptppola)
table(sweden_data$cptppola)##
## 1 2 3 4 5
## 1192 1564 1716 649 201model_data <- sweden_data %>%
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
TRUE ~ NA_character_ # Set anything that is not 1 or 2 to NA
),
inc = case_when(
hinctnta %in% 1:3 ~ "LowIncome",
hinctnta %in% 4:7 ~ "MidIncome",
hinctnta %in% 8:10 ~ "HighIncome",
hinctnta %in% c(77, 88, 99) ~ NA_character_,
TRUE ~ NA_character_
)
)#checking
table(model_data$inc)##
## HighIncome LowIncome MidIncome
## 4713 2319 4537#checking
table(model_data$cptppola)##
## 1 2 3 4 5
## 1192 1564 1716 649 201model1 = lm(cptppola ~ inc, data = model_data)
model2 = lm(cptppola ~ gndr + inc, data = model_data)
model3 = lm(cptppola ~ gndr + inc + gndr*inc, data = model_data)
modelsummary(
list(model1, model2, model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_rename = c("incLowIncome" = "Low Income", "incMidIncome" = "Mid Income", "gndrMale" = "Male", "gndrMale:incLowIncome" = "Male*Low income", "gndrMale:incMidIncome" = "Male*Mid Income"),
title = 'Table 2. Regression model')| (1) | (2) | (3) | |
|---|---|---|---|
| (Intercept) | 2.7 (0.0)*** | 2.6 (0.0)*** | 2.6 (0.0)*** |
| Low Income | −0.6 (0.0)*** | −0.5 (0.0)*** | −0.6 (0.1)*** |
| Mid Income | −0.4 (0.0)*** | −0.4 (0.0)*** | −0.4 (0.0)*** |
| Male | 0.2 (0.0)*** | 0.2 (0.0)*** | |
| Male*Low income | 0.1 (0.1) | ||
| Male*Mid Income | 0.0 (0.1) | ||
| Num.Obs. | 5079 | 5073 | 5073 |
| R2 | 0.040 | 0.053 | 0.053 |
| R2 Adj. | 0.040 | 0.052 | 0.052 |
| AIC | 15005.1 | 14917.3 | 14919.5 |
| BIC | 15031.3 | 14950.0 | 14965.2 |
| Log.Lik. | −7498.565 | −7453.668 | −7452.741 |
| RMSE | 1.06 | 1.05 | 1.05 |
interaction_plot <- effect("gndr*inc", model3, na.rm=TRUE)
plot(interaction_plot,
main="Figure 2. Interaction Model",
xlab="Income level",
ylab="confidence in political participation")
T-statistic
##Step 1 (calculate critical value, chi-squared value)
sweden_task <- sweden_data %>%
mutate(
inc = case_when(
hinctnta %in% 1:3 ~ "LowIncome",
hinctnta %in% 4:7 ~ "MidIncome",
hinctnta %in% 8:10 ~ "HighIncome",
hinctnta %in% c(77, 88, 99) ~ NA_character_,
TRUE ~ NA_character_
),
conf = case_when(
cptppola == 1 ~ "not at all confident",
cptppola == 2 ~ "a little confident",
cptppola == 3 ~ "quite confident",
cptppola == 4 ~ "very confident",
cptppola == 5 ~ "completely confident",
TRUE ~ NA_character_
)
) %>%
filter(!is.na(inc), !is.na(conf)) mytab <- table(sweden_task$inc, sweden_task$conf)
mytab##
## a little confident completely confident not at all confident
## HighIncome 514 92 236
## LowIncome 279 33 342
## MidIncome 690 69 552
##
## quite confident very confident
## HighIncome 694 311
## LowIncome 267 79
## MidIncome 687 234rsums <- rowSums(mytab)
csums <- colSums(mytab)
N <- sum(mytab)
expected <- outer(rsums, csums) / N
expected## a little confident completely confident not at all confident
## HighIncome 539.2993 70.54893 410.9293
## LowIncome 291.9866 38.19650 222.4847
## MidIncome 651.7141 85.25458 496.5859
## quite confident very confident
## HighIncome 599.3022 226.9203
## LowIncome 324.4733 122.8588
## MidIncome 724.2245 274.2209alpha <- 0.05
df5 <- (nrow(mytab) - 1) * (ncol(mytab) - 1)
critical_value <- qchisq(1 - alpha, df5)
observed <- mytab
chi_squared <- sum((observed - expected)^2 / expected)
p_value <- pchisq(chi_squared, df5, lower.tail = FALSE)
cat("Critical Value:", critical_value, "\n")## Critical Value: 15.50731cat("Chi-Squared Statistic:", chi_squared, "\n")## Chi-Squared Statistic: 238.9603test_stat <- sweden_task %>%
specify(explanatory = inc,
response = conf) %>%
hypothesize(null = "independence") %>%
calculate(stat = "chisq")
#stat = "t" or "z" presented errors
#chi-squared is also useful because variables are categorical (gender, income levels, confidence levels)print(test_stat$stat)## X-squared
## 238.9603##Step 2 (null)
null_distribution <- sweden_task %>%
specify(explanatory = inc,
response = conf) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "chisq")
null_distribution## Response: conf (factor)
## Explanatory: inc (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 5.85
## 2 2 6.54
## 3 3 10.1
## 4 4 22.2
## 5 5 12.8
## 6 6 7.61
## 7 7 10.5
## 8 8 5.57
## 9 9 5.23
## 10 10 8.81
## # ℹ 990 more rows##Step 3 (p_value)
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#cat("P-Value:", p_value, "\n")## P-Value: 3.759289e-47##Step 4 (visualize)
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided")## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
conf_int <- null_distribution%>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution %>%
visualize(data, bins = 15, method = "simulation", dens_color = "black") +
shade_p_value(obs_stat = test_stat, direction = "greater") +
shade_confidence_interval(endpoints = conf_int)## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf