Research Report
2024-04-06
Initial set-up
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "performance", "flextable", "broom", "scales", "ggeffects", "marginaleffects") # 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
## Loading required package: carData
##
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
##
## 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
##
##
##
## Attaching package: 'flextable'
##
##
## The following object is masked from 'package:purrr':
##
## compose
##
##
##
## Attaching package: 'scales'
##
##
## The following object is masked from 'package:purrr':
##
## discard
##
##
## The following object is masked from 'package:readr':
##
## col_factor## [[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] "performance" "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"
##
## [[8]]
## [1] "flextable" "performance" "survey" "survival" "Matrix"
## [6] "grid" "effects" "carData" "modelsummary" "fst"
## [11] "infer" "lubridate" "forcats" "stringr" "dplyr"
## [16] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [21] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [26] "datasets" "methods" "base"
##
## [[9]]
## [1] "broom" "flextable" "performance" "survey" "survival"
## [6] "Matrix" "grid" "effects" "carData" "modelsummary"
## [11] "fst" "infer" "lubridate" "forcats" "stringr"
## [16] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [21] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [26] "utils" "datasets" "methods" "base"
##
## [[10]]
## [1] "scales" "broom" "flextable" "performance" "survey"
## [6] "survival" "Matrix" "grid" "effects" "carData"
## [11] "modelsummary" "fst" "infer" "lubridate" "forcats"
## [16] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [21] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [26] "grDevices" "utils" "datasets" "methods" "base"
##
## [[11]]
## [1] "ggeffects" "scales" "broom" "flextable" "performance"
## [6] "survey" "survival" "Matrix" "grid" "effects"
## [11] "carData" "modelsummary" "fst" "infer" "lubridate"
## [16] "forcats" "stringr" "dplyr" "purrr" "readr"
## [21] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [26] "graphics" "grDevices" "utils" "datasets" "methods"
## [31] "base"
##
## [[12]]
## [1] "marginaleffects" "ggeffects" "scales" "broom"
## [5] "flextable" "performance" "survey" "survival"
## [9] "Matrix" "grid" "effects" "carData"
## [13] "modelsummary" "fst" "infer" "lubridate"
## [17] "forcats" "stringr" "dplyr" "purrr"
## [21] "readr" "tidyr" "tibble" "ggplot2"
## [25] "tidyverse" "stats" "graphics" "grDevices"
## [29] "utils" "datasets" "methods" "base"Loading dataset
ess <- read_fst("All-ESS-Data.fst")spain_data <- ess %>%
filter(cntry == "ES")
spain_data$weight <- spain_data$dweight * spain_data$pweight
survey_design <- svydesign(ids = ~1, data = spain_data, weights = ~weight)
spain_data <- spain_data %>%
mutate(protest = case_when(
pbldmn == 1 ~ 1, # Yes to 1
pbldmn == 2 ~ 0, # No to 0
TRUE ~ NA_real_ # Everything else to NA
))
spain_data <- spain_data %>% filter(!is.na(protest))Spain by itself
df <- spain_data
df$year <- NA
replacements <- c(2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020)
for(i in 1:10){
df$year[df$essround == i] <- replacements[i]
}
protest_avg <- df %>%
group_by(year) %>%
summarize(protest_avg = mean(protest, na.rm = TRUE))
protest_plot_SE <- ggplot(protest_avg, aes(x = year, y = protest_avg)) +
geom_point(aes(size = protest_avg, color = protest_avg), alpha = 0.7) +
geom_line(aes(group = 1), color = "blue", linetype = "dashed") +
labs(
x = "Survey Year",
y = "Pr(Protest)") +
theme_minimal() +
theme(legend.position = "none")
print(protest_plot_SE)
spain_data <- spain_data %>%
# Recoding for feeling about household's income
mutate(feelincome = case_when(
hincfel == 1 ~ "comfortable", # Living comfortably to "comfortable"
hincfel %in% 2:4 ~ "difficult", # Other categories to "difficult"
TRUE ~ NA_character_ # Everything else to NA
)) %>%
# Recoding age, adding age group, and cohort with generational cohorts
mutate(
age = ifelse(agea == 999, NA, agea),
adulthood = ifelse(age >= 18 & age <= 30, "Early", ifelse(age >= 31, "Later", NA)),
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
gen = case_when(
yrbrn %in% 1900:1945 ~ "Interwar",
yrbrn %in% 1946:1964 ~ "Baby Boomers",
yrbrn %in% 1965:1979 ~ "Gen X",
yrbrn %in% 1980:1996 ~ "Millennials",
TRUE ~ as.character(cohort)
),
gen = factor(gen, levels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
) %>%
# Recoding political ID and education
mutate(
polID = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
lrscale %in% 4:6 ~ "Moderate",
lrscale %in% c(77, 88, 99) ~ NA_character_
),
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
),
edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = factor(educ.ba, levels = c("No BA", "BA or more"))
)
# Setting values greater than 10 to NA
spain_data$trstplt <- ifelse(spain_data$trstplt > 10, NA, spain_data$trstplt)
spain_data$trstprl <- ifelse(spain_data$trstprl > 10, NA, spain_data$trstprl)
spain_data$trstprt <- ifelse(spain_data$trstprt > 10, NA, spain_data$trstprt)
# Creating and rescaling the trust variable
spain_data$trust <- scales::rescale(spain_data$trstplt + spain_data$trstprl + spain_data$trstprt, na.rm = TRUE, to = c(0, 100))
# Here's how you would "flip it" because populist is conceived by N + I as 'anti-trust'
# As Schafer explains there are some issues with that conceptualization, but N + I is also widely applied
spain_data$populist <- scales::rescale(spain_data$trust, na.rm = TRUE, to=c(100,0))
table(df$educ.ba)## < table of extent 0 > table(spain_data$populist)##
## 0 3.33333333333333 6.66666666666667 10
## 20 2 9 18
## 13.3333333333333 16.6666666666667 20 23.3333333333333
## 35 41 123 96
## 26.6666666666667 30 33.3333333333333 36.6666666666667
## 134 203 266 376
## 40 43.3333333333333 46.6666666666667 50
## 481 554 615 1017
## 53.3333333333333 56.6666666666667 60 63.3333333333333
## 634 757 660 765
## 66.6666666666667 70 73.3333333333333 76.6666666666667
## 730 742 595 698
## 80 83.3333333333333 86.6666666666667 90
## 592 827 601 612
## 93.3333333333333 96.6666666666667 100
## 474 298 1448spaindf <- spain_data %>%
dplyr::select(protest, feelincome, adulthood, gen, populist, polID, educ.ba, weight)
# Check the first few rows to verify the selection
spaindf <- na.omit(spaindf)
head(spaindf)## protest feelincome adulthood gen populist polID educ.ba
## 1708 0 difficult Later Baby Boomers 83.33333 Moderate No BA
## 1709 0 difficult Later Baby Boomers 90.00000 Left No BA
## 1710 1 comfortable Later Gen X 70.00000 Moderate No BA
## 1711 0 difficult Later Interwar 23.33333 Right No BA
## 1712 0 comfortable Later Baby Boomers 33.33333 Right No BA
## 1713 1 difficult Later Gen X 36.66667 Left No BA
## weight
## 1708 1.415770
## 1709 1.500677
## 1710 1.500677
## 1711 2.377397
## 1712 2.377397
## 1713 2.377397model1 <- glm(protest ~ feelincome + populist, family = binomial, data = spaindf, weights = round(weight))
model2 <- glm(protest ~ polID + adulthood, family = binomial, data = spaindf, weights = round(weight))
model3 <- glm(protest ~ gen + educ.ba, family = binomial, data = spaindf, weights = round(weight))
modelsummary(
list(model1, model2, model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_omit = "Intercept")| Â (1) | Â Â (2) | Â Â (3) | |
|---|---|---|---|
| feelincomedifficult | −0.3 (0.0)*** | ||
| populist | 0.0 (0.0) | ||
| polIDModerate | −0.9 (0.0)*** | ||
| polIDRight | −1.3 (0.1)*** | ||
| adulthoodLater | −0.3 (0.0)*** | ||
| genBaby Boomers | 1.0 (0.1)*** | ||
| genGen X | 0.8 (0.1)*** | ||
| genMillennials | 1.0 (0.1)*** | ||
| educ.baBA or more | 0.7 (0.0)*** | ||
| Num.Obs. | 12048 | 12048 | 12048 |
| AIC | 26139.4 | 25101.0 | 25371.6 |
| BIC | 26161.6 | 25130.6 | 25408.6 |
| Log.Lik. | −13066.687 | −12546.498 | −12680.789 |
| RMSE | 0.42 | 0.41 | 0.42 |
Descriptive table
spain_data <- spain_data %>%
mutate(
protest = protest,
ideology = case_when(
lrscale > 10 ~ NA_real_, # Set values above 10 to NA
TRUE ~ lrscale # Keep other values as is
),
age = case_when(
agea > 100 ~ NA_real_, # Set values above 100 to NA
agea < 15 ~ NA_real_, # Set values below 15 to NA
TRUE ~ agea # Keep other values as is
)
)
df <- spain_data
df <- df %>% filter(!is.na(ideology))
df <- df %>% filter(!is.na(age))
table1b <- datasummary_skim(df %>% dplyr::select(protest, ideology, age), title = "Table 1. Descriptive statistics for main variables", 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.table1b
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
protest | 2 | 0 | 0.2 | 0.4 | 0.0 | 0.0 | 1.0 |
ideology | 11 | 0 | 4.5 | 2.1 | 0.0 | 5.0 | 10.0 |
age | 84 | 0 | 47.2 | 18.3 | 15.0 | 46.0 | 99.0 |
df <- df %>%
mutate(
protest = case_when(
protest == 1 ~ "Yes", # 1 to 'Yes'
protest == 0 ~ "No", # 0 to 'No'
TRUE ~ NA_character_ ), # Keep NA values as is, but ensure character type
ideology = case_when(
lrscale >= 0 & lrscale <= 2 ~ "Extreme Left",
lrscale > 2 & lrscale <= 4 ~ "Moderate Left",
lrscale == 5 ~ "Centrist",
lrscale >= 6 & lrscale <= 8 ~ "Moderate Right",
lrscale >= 9 & lrscale <= 10 ~ "Extreme Right",
TRUE ~ NA_character_ # For any cases that do not match above conditions
),
age = ifelse(agea == 999, NA, agea), # Replace 999 with NA in age
adulthood = case_when(
age >= 18 & age <= 30 ~ "Early",
age >= 31 ~ "Later",
TRUE ~ NA_character_ # For ages outside the specified ranges or NA
)
)
table(df$protest)##
## No Yes
## 11331 3354table(df$ideology)##
## Centrist Extreme Left Extreme Right Moderate Left Moderate Right
## 4560 2370 499 4162 3094table(df$adulthood)##
## Early Later
## 2760 11570Hypothesis testing with infer
Step 1
test_stat <- df %>%
specify(explanatory = protest, # change variable name for explanatory variable
response = ideology) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "Chisq")
print(test_stat$stat)## X-squared
## 786.4693The Chi-squared test of independence is suitable for my data because it allows me to determine if there is a relationship between if people protest and their ideological affiliations, which are both categorical variables.
Step 2
null_dist <- df %>%
specify(response = protest, explanatory = ideology) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")
null_dist## Response: protest (factor)
## Explanatory: ideology (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 9.74
## 2 2 3.80
## 3 3 2.86
## 4 4 6.12
## 5 5 1.34
## 6 6 5.40
## 7 7 0.814
## 8 8 5.60
## 9 9 13.8
## 10 10 2.95
## # ℹ 990 more rowsStep 3
p_val <- null_dist %>% #
get_pvalue(obs_stat = test_stat, direction = "greater")## 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 0Step 4
conf_int <- null_dist%>%
get_confidence_interval(level = 0.95, type = "percentile")
null_dist %>%
visualize(data, bins = 10, 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
# Assuming 'test_stat' is your observed test statistic value
test_stat <- 960
conf_int <- null_dist %>%
get_confidence_interval(level = 0.95, type = "percentile")
# Visualize the null distribution with manual x-axis limits
null_dist %>%
visualize(bins = 15, method = "simulation", dens_color = "black") +
geom_vline(aes(xintercept = test_stat), color = "red", linetype = "dashed", size = 1) +
shade_p_value(obs_stat = test_stat, direction = "greater") +
shade_confidence_interval(endpoints = conf_int) +
coord_cartesian(xlim = c(0, max(null_dist$stat) * 1.1)) + # Set manual x-axis limits
labs(title = "Simulation-Based Null Distribution",
x = "Chi-Square",
y = "Count")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
Modelling
# If running your GLM, need to put your outcome as 0 and 1
df <- df %>%
mutate(
protest = case_when(
protest == "Yes" ~ 1,
protest == "No" ~ 0,
TRUE ~ NA_integer_
)
)Regression models table
model1 <- glm(protest ~ ideology, family = binomial, data = df, weights = round(weight)) # single predictor
model2 <- glm(protest ~ ideology + adulthood, family = binomial, data = df, weights = round(weight)) # two predictors
model3 <- glm(protest ~ ideology + adulthood + ideology*adulthood, family = binomial, data = df, weights = round(weight)) # both predictors + interactionmodelsummary(
list(model1, model2, model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL)| Â (1) | Â Â (2) | Â Â (3) | |
|---|---|---|---|
| (Intercept) | −1.6 (0.0)*** | −1.4 (0.0)*** | −1.4 (0.1)*** |
| ideologyExtreme Left | 1.3 (0.0)*** | 1.3 (0.0)*** | 1.3 (0.1)*** |
| ideologyExtreme Right | −0.2 (0.1)* | −0.2 (0.1)+ | −0.1 (0.2) |
| ideologyModerate Left | 0.7 (0.0)*** | 0.7 (0.0)*** | 0.8 (0.1)*** |
| ideologyModerate Right | −0.2 (0.0)*** | −0.2 (0.0)*** | −0.4 (0.1)*** |
| adulthoodLater | −0.4 (0.0)*** | −0.3 (0.1)*** | |
| ideologyExtreme Left × adulthoodLater | 0.0 (0.1) | ||
| ideologyExtreme Right × adulthoodLater | −0.1 (0.3) | ||
| ideologyModerate Left × adulthoodLater | −0.1 (0.1) | ||
| ideologyModerate Right × adulthoodLater | 0.2 (0.1)+ | ||
| Num.Obs. | 14490 | 14141 | 14141 |
| AIC | 29616.6 | 28644.3 | 28643.5 |
| BIC | 29654.5 | 28689.6 | 28719.1 |
| Log.Lik. | −14803.315 | −14316.137 | −14311.760 |
| RMSE | 0.41 | 0.40 | 0.40 |
Improving display
models <- list()
models[['SLR']]<- glm(protest ~ ideology, family = binomial, data = df, weights = round(weight)) # single predictor
models[['MLR']]<- glm(protest ~ ideology + adulthood, family = binomial, data = df, weights = round(weight)) # two predictors
models[['Interaction']]<- glm(protest ~ ideology + adulthood + ideology*adulthood, family = binomial, data = df, weights = round(weight)) # both predictors + interaction
modelsummary(models,
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
exponentiate = TRUE, ## NOTE: this is to exponentiate directly, otherwise leave out or set to FALSE
coef_rename = c("ideologyExtreme Left" = "Extreme Left", "ideologyExtreme Right" = "Extreme Right", "ideologyModerate Left" = "Moderate Left", "ideologyModerate Right" = "Moderate Right", "adulthoodLater" = "Later Adulthood"),
title = 'Table 2. Regression models predicting protest')| SLR | MLR | Â Interaction | |
|---|---|---|---|
| (Intercept) | 0.2 (0.0)*** | 0.3 (0.0)*** | 0.3 (0.0)*** |
| Extreme Left | 3.6 (0.1)*** | 3.5 (0.1)*** | 3.6 (0.3)*** |
| Extreme Right | 0.8 (0.1)* | 0.8 (0.1)+ | 0.9 (0.2) |
| Moderate Left | 2.0 (0.1)*** | 2.0 (0.1)*** | 2.2 (0.2)*** |
| Moderate Right | 0.8 (0.0)*** | 0.8 (0.0)*** | 0.7 (0.1)*** |
| Later Adulthood | 0.7 (0.0)*** | 0.7 (0.0)*** | |
| Extreme Left:Later Adulthood | 1.0 (0.1) | ||
| Extreme Right:Later Adulthood | 0.9 (0.2) | ||
| Moderate Left:Later Adulthood | 0.9 (0.1) | ||
| Moderate Right:Later Adulthood | 1.3 (0.2)+ | ||
| Num.Obs. | 14490 | 14141 | 14141 |
| AIC | 29616.6 | 28644.3 | 28643.5 |
| BIC | 29654.5 | 28689.6 | 28719.1 |
| Log.Lik. | −14803.315 | −14316.137 | −14311.760 |
| RMSE | 0.41 | 0.40 | 0.40 |
Visual
plot_comparisons(model3, variables = "ideology", condition = "adulthood")