Final_Report
2024-04-06
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.
Setting up the environment:
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "performance", "flextable", "broom", "scales", "ggeffects", "marginaleffects", "interactions", "see") # 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## 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 'performance' was built under R version 4.3.3##
## 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## Warning: package 'ggeffects' was built under R version 4.3.3## Warning: package 'marginaleffects' was built under R version 4.3.3## Warning: package 'interactions' was built under R version 4.3.3##
## Attaching package: 'interactions'
##
## The following object is masked from 'package:ggeffects':
##
## johnson_neyman## Warning: package 'see' 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] "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"
##
## [[13]]
## [1] "interactions" "marginaleffects" "ggeffects" "scales"
## [5] "broom" "flextable" "performance" "survey"
## [9] "survival" "Matrix" "grid" "effects"
## [13] "carData" "modelsummary" "fst" "infer"
## [17] "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble"
## [25] "ggplot2" "tidyverse" "stats" "graphics"
## [29] "grDevices" "utils" "datasets" "methods"
## [33] "base"
##
## [[14]]
## [1] "see" "interactions" "marginaleffects" "ggeffects"
## [5] "scales" "broom" "flextable" "performance"
## [9] "survey" "survival" "Matrix" "grid"
## [13] "effects" "carData" "modelsummary" "fst"
## [17] "infer" "lubridate" "forcats" "stringr"
## [21] "dplyr" "purrr" "readr" "tidyr"
## [25] "tibble" "ggplot2" "tidyverse" "stats"
## [29] "graphics" "grDevices" "utils" "datasets"
## [33] "methods" "base"Loading data for Italy:
Italy_data <- read.fst("italy_data.fst")code our main variable of interest:
Utilizing the variable pbldmn ()
Italy_data$weight <- Italy_data$dweight * Italy_data$pweight
survey_design <- svydesign(ids = ~1, data = Italy_data, weights = ~weight)
Italy_data <- Italy_data %>%
mutate(protest = case_when(
pbldmn == 1 ~ 1, # Yes to 1
pbldmn == 2 ~ 0, # No to 0
TRUE ~ NA_real_ # Everything else to NA
))
Italy_data <- Italy_data %>% filter(!is.na(protest))Coding predictors:
Italy_data <- Italy_data %>%
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"))
)Categorical descriptive table
Italy_datadf <- Italy_data %>%
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
))
Italy_datadf <- Italy_datadf %>%
# Relabeling variables for clarity
mutate(
Protest = protest,
Ideology = polID,
Generation = gen,
)
table1 <- datasummary_skim(Italy_datadf %>% dplyr::select(Protest, Ideology, Generation), type = "categorical",
title = "Table 1. Descriptive statistics for main variables", output = "flextable")
table1
|
| N | % |
|---|---|---|---|
Protest | No | 6715 | 90.1 |
Yes | 741 | 9.9 | |
Ideology | Left | 1391 | 18.7 |
Moderate | 2527 | 33.9 | |
Right | 1483 | 19.9 | |
NA | 2055 | 27.6 | |
Generation | Interwar | 1200 | 16.1 |
Baby Boomers | 2343 | 31.4 | |
Gen X | 1927 | 25.8 | |
Millennials | 1475 | 19.8 |
print(table1, preview ="pdf")## a flextable object.
## col_keys: ` `, ` `, `N`, `%`
## header has 1 row(s)
## body has 10 row(s)
## original dataset sample:
## N %
## 1 Protest No 6715 90.1
## 2 Yes 741 9.9
## 3 Ideology Left 1391 18.7
## 4 Moderate 2527 33.9
## 5 Right 1483 19.9Numeric descriptive table
Italy_data <- Italy_data %>%
mutate(
Protest = protest,
Ideology = case_when(
lrscale > 10 ~ NA_real_, # Set values above 10 to NA
TRUE ~ lrscale # Keep other values as is
),
Generation = case_when(
yrbrn > 2000 ~ NA_real_, # Set values above 100 to NA
yrbrn < 1920 ~ NA_real_, # Set values below 15 to NA
TRUE ~ yrbrn # Keep other values as is
)
)
df <- Italy_data
df <- df %>% filter(!is.na(Ideology))
df <- df %>% filter(!is.na(Generation))
table1b <- datasummary_skim(df %>% dplyr::select(Protest, Ideology, Generation), title = "Table 2. 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.1 | 0.3 | 0.0 | 0.0 | 1.0 |
Ideology | 11 | 0 | 5.1 | 2.5 | 0.0 | 5.0 | 10.0 |
Generation | 81 | 0 | 1965.1 | 18.3 | 1920.0 | 1965.0 | 2000.0 |
Hypothesis testing with infer:
step 1:
test_stat <- Italy_datadf %>%
specify(explanatory = Generation, # change variable name for explanatory variable
response = Protest) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "Chisq") ## Warning: Removed 511 rows containing missing values.print(test_stat$stat)## X-squared
## 49.58601null_dist <- Italy_datadf %>%
specify(response = Protest, explanatory = Ideology) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")## Warning: Removed 2055 rows containing missing values.null_dist## Response: Protest (factor)
## Explanatory: Ideology (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 3.88
## 2 2 7.08
## 3 3 5.32
## 4 4 0.136
## 5 5 0.468
## 6 6 0.814
## 7 7 3.83
## 8 8 1.54
## 9 9 0.431
## 10 10 0.617
## # ℹ 990 more rowsp_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 0conf_int <- null_dist%>%
get_confidence_interval(level = 0.95, type = "percentile")
figure2 <- 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)
print(figure2)## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
ggsave("output_figure.pdf", plot = figure2, device = "pdf")## Saving 7 x 5 in image## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Infprint(figure2, preview = "pdf")## 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 <- 49.58601
conf_int <- null_dist %>%
get_confidence_interval(level = 0.95, type = "percentile")
# Visualize the null distribution with manual x-axis limits
figure1 <- 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 = "Figure 1. 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.print(figure1)## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
print(figure1, preview= "PDF")## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
Regression modelling:
Italy_datadf <- Italy_datadf %>%
mutate(
Protest = case_when(
Protest == "Yes" ~ 1,
Protest == "No" ~ 0,
TRUE ~ NA_integer_
)
)model1 <- glm(Protest ~ Generation, family = binomial, data = Italy_datadf, weights = round(weight)) # single predictor
model2 <- glm(Protest ~ Generation + Ideology, family = binomial, data = Italy_datadf, weights = round(weight)) # two predictors
model3 <- glm(Protest ~ Generation + Ideology + Generation*Ideology, family = binomial, data = Italy_datadf, weights = round(weight)) # both predictors + interactionModel Summary:
modelsummary(
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) | -2.9 (0.1)*** | -2.0 (0.1)*** | -2.0 (0.1)*** |
GenerationBaby Boomers | 0.8 (0.1)*** | 0.7 (0.1)*** | 0.7 (0.1)*** |
GenerationGen X | 0.8 (0.1)*** | 0.7 (0.1)*** | 0.7 (0.1)*** |
GenerationMillennials | 1.1 (0.1)*** | 1.1 (0.1)*** | 1.1 (0.1)*** |
IdeologyModerate | -0.9 (0.1)*** | -1.0 (0.2)*** | |
IdeologyRight | -0.9 (0.1)*** | -0.9 (0.2)*** | |
GenerationBaby Boomers × IdeologyModerate | 0.1 (0.2) | ||
GenerationGen X × IdeologyModerate | 0.1 (0.2) | ||
GenerationMillennials × IdeologyModerate | 0.0 (0.2) | ||
GenerationBaby Boomers × IdeologyRight | 0.0 (0.2) | ||
GenerationGen X × IdeologyRight | 0.0 (0.2) | ||
GenerationMillennials × IdeologyRight | -0.1 (0.2) | ||
Num.Obs. | 6945 | 5104 | 5104 |
AIC | 12779.9 | 10545.5 | 10555.3 |
BIC | 12807.3 | 10584.7 | 10633.8 |
Log.Lik. | -6385.936 | -5266.741 | -5265.666 |
RMSE | 0.29 | 0.32 | 0.32 |
Improving display:
models <- list()
models[['SLR']]<- glm(Protest ~ Generation, family = binomial, data = Italy_datadf, weights = round(weight)) # single predictor
models[['MLR']]<- glm(Protest ~ Generation + Ideology, family = binomial, data = Italy_datadf, weights = round(weight)) # two predictors
models[['Interaction']]<- glm(Protest ~ Generation + Ideology + Generation*Ideology, family = binomial, data = Italy_datadf, weights = round(weight)) # both predictors + interaction
pdf("output_table.pdf", width = 10, height = 6)
modelsummary(models,
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
exponentiate = FALSE, ## NOTE: this is to exponentiate directly, otherwise leave out or set to FALSE
coef_rename = c("GenerationBaby Boomers" = "Baby Boomers", "GenerationGen X" = "Gen X", "GenerationMillennials" = "Millennials", "IdeologyModerate" = "Moderate", "IdeologyRight" = "Right"),
title = 'Table 3. Regression models predicting protest')
| SLR | MLR | Interaction |
|---|---|---|---|
(Intercept) | -2.9 (0.1)*** | -2.0 (0.1)*** | -2.0 (0.1)*** |
Baby Boomers | 0.8 (0.1)*** | 0.7 (0.1)*** | 0.7 (0.1)*** |
Gen X | 0.8 (0.1)*** | 0.7 (0.1)*** | 0.7 (0.1)*** |
Millennials | 1.1 (0.1)*** | 1.1 (0.1)*** | 1.1 (0.1)*** |
Moderate | -0.9 (0.1)*** | -1.0 (0.2)*** | |
Right | -0.9 (0.1)*** | -0.9 (0.2)*** | |
Baby Boomers:Moderate | 0.1 (0.2) | ||
Gen X:Moderate | 0.1 (0.2) | ||
Millennials:Moderate | 0.0 (0.2) | ||
Baby Boomers:Right | 0.0 (0.2) | ||
Gen X:Right | 0.0 (0.2) | ||
Millennials:Right | -0.1 (0.2) | ||
Num.Obs. | 6945 | 5104 | 5104 |
AIC | 12779.9 | 10545.5 | 10555.3 |
BIC | 12807.3 | 10584.7 | 10633.8 |
Log.Lik. | -6385.936 | -5266.741 | -5265.666 |
RMSE | 0.29 | 0.32 | 0.32 |
pdf("output_table.pdf")
print('Table 3. Regression models predicting protest')## [1] "Table 3. Regression models predicting protest"dev.off()## png
## 2Interaction plot:
model3<- glm(Protest ~ Generation + Ideology + Generation*Ideology, family = binomial, data = Italy_datadf, weights = round(weight))
cat_plot(model3, pred = Generation, modx = Ideology , geom= "line", jnplot = TRUE) 
title = ('Figure 2. Catplot')
ggsave("catplot.pdf", plot = last_plot(), device = "pdf")## Saving 7 x 5 in imagecheck_model(model1)