FinalSoc222
Yizhan Su
2024-04-12
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "performance", "flextable", "broom", "see","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.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## Warning: package 'flextable' was built under R version 4.3.3##
## Attaching package: 'flextable'
##
## The following object is masked from 'package:purrr':
##
## compose## Warning: package 'see' was built under R version 4.3.3##
## 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## [[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] "see" "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] "scales" "see" "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] "ggeffects" "scales" "see" "broom" "flextable"
## [6] "performance" "survey" "survival" "Matrix" "grid"
## [11] "effects" "carData" "modelsummary" "fst" "infer"
## [16] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [21] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [26] "stats" "graphics" "grDevices" "utils" "datasets"
## [31] "methods" "base"
##
## [[13]]
## [1] "marginaleffects" "ggeffects" "scales" "see"
## [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"ess <- read_fst("All-ESS-Data.fst")italy_data <- ess %>%
filter(cntry == "IT")
italy_data$weight <- italy_data$dweight * italy_data$pweight
survey_design <- svydesign(ids = ~1, data = italy_data, weights = ~weight)
italy_data <- italy_data %>%
mutate(participation = case_when(
pstplonl == 1 ~ 1, # Yes to 1
pstplonl == 2 ~ 0, # No to 0
TRUE ~ NA_real_ # Everything else to NA
))
italy_data <- italy_data %>% filter(!is.na(pstplonl))df <- italy_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]
}
pstplonl_avg <- df %>%
group_by(year) %>%
summarize(pstplonl_avg = mean(pstplonl, na.rm = TRUE))
pstplonl_plot_IT <- ggplot(pstplonl_avg, aes(x = year, y = pstplonl_avg)) +
geom_point(aes(size = pstplonl_avg, color = pstplonl_avg), alpha = 0.7) +
geom_line(aes(group = 1), color = "blue", linetype = "dashed") +
labs(
x = "Survey Year",
y = "participation(pstplonl)") +
theme_minimal() +
theme(legend.position = "none")
print(pstplonl_plot_IT)
ggsave("Fig1.pdf", plot = pstplonl_plot_IT, width = 10, height = 6)ess_data <- ess
replacements <- c(2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020)
ess_data$year <- NA
for(i in 1:10){
ess_data$year[ess_data$essround == i] <- replacements[i]
}
ess_data <- ess_data %>%
mutate(participation = case_when(
pstplonl == 1 ~ 1, # Yes to 1
pstplonl == 2 ~ 0, # No to 0
TRUE ~ NA_real_
)) %>%
filter(!is.na(pstplonl))
# Now ess_data is your ess_cleaned with 'year' correctly assigned
ess_cleaned <- ess_data %>%
dplyr::select(cntry, year, pstplonl) %>%
filter(!is.na(year))
# Calculating yearly averages for each country
yearly_pstplonl_avg <- ess_cleaned %>%
group_by(cntry, year) %>%
summarize(avg_pstplonl = mean(pstplonl, na.rm = TRUE), .groups = 'drop')
# Subset for specific countries and overall ESS baseline
specific_countries_pstplonl <- yearly_pstplonl_avg %>%
filter(cntry %in% c("GB", "FR", "ES", "DE", "IT"))
ess_baseline_pstplonl <- ess_cleaned %>%
group_by(year) %>%
summarize(avg_pstplonl = mean(pstplonl, na.rm = TRUE), .groups = 'drop') %>%
mutate(cntry = "ESS Baseline")
# Combining specific countries with ESS baseline
combined_data_pstplonl <- rbind(specific_countries_pstplonl, ess_baseline_pstplonl)
# Mapping of country codes to names for better label clarity
country_labels <- c(GB = "Great Britain", FR = "France", ES = "Spain", DE = "Germany", IT = "Italy", `ESS Baseline` = "ESS Baseline")
# Plotting
pstplonl_plot <- ggplot(combined_data_pstplonl, aes(x = year, y = avg_pstplonl, group = cntry, color = cntry)) +
geom_line(size = 1) + # Thicker lines for better visibility
scale_color_manual(values = c(GB = "darkgreen", FR = "red", ES = "yellow", DE = "orange", IT = "blue", `ESS Baseline` = "black"),
labels = country_labels) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
labs(
x = "Year",
y = "participation(pstplonl)",
color = "Country") +
theme_minimal() +
theme(legend.position = "right", legend.title.align = 0.5) +
guides(color = guide_legend(override.aes = list(linewidth = 5)))## 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(pstplonl_plot)
ggsave("Fig2.pdf", plot = pstplonl_plot, width = 10, height = 6, device = "pdf")italy_data <- italy_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
italy_data$trstplt <- ifelse(italy_data$trstplt > 10, NA, italy_data$trstplt)
italy_data$trstprl <- ifelse(italy_data$trstprl > 10, NA, italy_data$trstprl)
italy_data$trstprt <- ifelse(italy_data$trstprt > 10, NA, italy_data$trstprt)
# Creating and rescaling the trust variable
italy_data$trust <- scales::rescale(italy_data$trstplt + italy_data$trstprl + italy_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
italy_data$populist <- scales::rescale(italy_data$trust, na.rm = TRUE, to=c(100,0))
table(italy_data$populist)##
## 0 3.33333333333333 6.66666666666667 10
## 11 1 7 16
## 13.3333333333333 16.6666666666667 20 23.3333333333333
## 16 14 54 38
## 26.6666666666667 30 33.3333333333333 36.6666666666667
## 77 129 136 199
## 40 43.3333333333333 46.6666666666667 50
## 319 259 349 427
## 53.3333333333333 56.6666666666667 60 63.3333333333333
## 304 396 380 354
## 66.6666666666667 70 73.3333333333333 76.6666666666667
## 351 390 331 348
## 80 83.3333333333333 86.6666666666667 90
## 380 316 292 363
## 93.3333333333333 96.6666666666667 100
## 262 235 1000italydf <- italy_data %>%
dplyr::select(participation, feelincome, adulthood, gen, populist, polID, educ.ba, weight)
# Check the first few rows to verify the selection
italydf <- na.omit(italydf)
head(italydf)## participation feelincome adulthood gen populist polID educ.ba
## 1 0 difficult Later Baby Boomers 80.00000 Moderate BA or more
## 3 0 difficult Early Millennials 90.00000 Moderate No BA
## 5 1 difficult Later Gen X 96.66667 Moderate No BA
## 6 0 difficult Early Millennials 36.66667 Moderate No BA
## 8 1 difficult Later Interwar 56.66667 Moderate BA or more
## 11 0 difficult Later Gen X 53.33333 Right No BA
## weight
## 1 2.078556
## 3 2.028642
## 5 1.966716
## 6 1.986275
## 8 2.097367
## 11 1.966716model1 <- glm(participation ~ feelincome + populist, family = binomial, data = italydf, weights = round(weight))
model2 <- glm(participation ~ polID + adulthood, family = binomial, data = italydf, weights = round(weight))
model3 <- glm(participation ~ gen + educ.ba, family = binomial, data = italydf, 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.5 (0.1)*** | ||
| populist | 0.0 (0.0)*** | ||
| polIDModerate | −0.6 (0.1)*** | ||
| polIDRight | −0.4 (0.1)*** | ||
| adulthoodLater | −0.6 (0.1)*** | ||
| genBaby Boomers | 1.5 (0.2)*** | ||
| genGen X | 1.9 (0.2)*** | ||
| genMillennials | 2.1 (0.2)*** | ||
| educ.baBA or more | 0.7 (0.1)*** | ||
| Num.Obs. | 4990 | 4990 | 4990 |
| AIC | 9224.2 | 9194.5 | 8826.2 |
| BIC | 9243.7 | 9220.6 | 8858.8 |
| Log.Lik. | −4609.101 | −4593.252 | −4408.105 |
| RMSE | 0.37 | 0.37 | 0.37 |
italydf <- italydf %>%
mutate(participation = case_when(
participation == 1 ~ "Yes", # 1 to 'Yes'
participation == 0 ~ "No", # 0 to 'No'
TRUE ~ NA_character_ # Keep NA values as is, but ensure character type
))
italyt1 <- italydf %>%
# Relabeling variables for clarity
mutate(
participation = participation,
Ideology = polID,
Adulthood = adulthood
)
table1 <- datasummary_skim(italyt1 %>% dplyr::select(participation, Ideology, Adulthood), type = "categorical",
title = "Table 1. Descriptive statistics for main variables", output = "flextable")
table1
|
| N | % |
|---|---|---|---|
participation | No | 4145 | 83.1 |
Yes | 845 | 16.9 | |
Ideology | Left | 1193 | 23.9 |
Moderate | 2321 | 46.5 | |
Right | 1476 | 29.6 | |
Adulthood | Early | 516 | 10.3 |
Later | 4474 | 89.7 |
set_flextable_defaults(fonts_ignore=TRUE)
flextable::save_as_docx(table1, path = "Table1.docx", # change name to whatever you would like
width = 7.0, height = 7.0) # play around with dimensionsdf <- df %>%
mutate(
participation = case_when(
participation == 1 ~ "Yes", # 1 to 'Yes'
participation == 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$participation)##
## No Yes
## 6749 1143table(df$Ideology)##
## Centrist Extreme Left Extreme Right Moderate Left Moderate Right
## 1399 732 450 1235 1848table(df$Adulthood)##
## Early Later
## 1306 6425test_stat <- df %>%
specify(explanatory = participation, # change variable name for explanatory variable
response = Ideology) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "Chisq") ## Warning: Removed 2406 rows containing missing values.print(test_stat$stat)## X-squared
## 59.30867null_dist <- df %>%
specify(response = participation, explanatory = Ideology) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")## Warning: Removed 2406 rows containing missing values.null_dist## Response: participation (factor)
## Explanatory: Ideology (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 2.70
## 2 2 1.66
## 3 3 3.01
## 4 4 2.55
## 5 5 4.36
## 6 6 5.02
## 7 7 1.86
## 8 8 3.30
## 9 9 0.778
## 10 10 1.26
## # ℹ 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")
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 in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
df <- df %>%
mutate(
participation = case_when(
participation == "Yes" ~ 1,
participation == "No" ~ 0,
TRUE ~ NA_integer_
)
)model1 <- glm(participation ~ Ideology, family = binomial, data = df, weights = round(weight)) # single predictor
model2 <- glm(participation ~ Ideology + Adulthood, family = binomial, data = df, weights = round(weight)) # two predictors
model3 <- glm(participation ~ 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.8 (0.1)*** | −1.3 (0.1)*** | −1.3 (0.1)*** |
| IdeologyExtreme Left | 0.7 (0.1)*** | 0.7 (0.1)*** | 0.9 (0.2)*** |
| IdeologyExtreme Right | 0.5 (0.1)*** | 0.5 (0.1)*** | 0.6 (0.2)** |
| IdeologyModerate Left | 0.4 (0.1)*** | 0.4 (0.1)*** | 0.4 (0.2)** |
| IdeologyModerate Right | 0.0 (0.1) | 0.0 (0.1) | −0.2 (0.2) |
| AdulthoodLater | −0.6 (0.1)*** | −0.6 (0.1)*** | |
| IdeologyExtreme Left × AdulthoodLater | −0.3 (0.2) | ||
| IdeologyExtreme Right × AdulthoodLater | −0.1 (0.3) | ||
| IdeologyModerate Left × AdulthoodLater | 0.0 (0.2) | ||
| IdeologyModerate Right × AdulthoodLater | 0.3 (0.2) | ||
| Num.Obs. | 5605 | 5476 | 5476 |
| AIC | 10569.6 | 10225.9 | 10225.7 |
| BIC | 10602.8 | 10265.5 | 10291.7 |
| Log.Lik. | −5279.816 | −5106.945 | −5102.830 |
| RMSE | 0.38 | 0.38 | 0.38 |
models <- list()
models[['SLR']]<- glm(participation ~ Ideology, family = binomial, data = df, weights = round(weight)) # single predictor
models[['MLR']]<- glm(participation ~ Ideology + Adulthood, family = binomial, data = df, weights = round(weight)) # two predictors
models[['Interaction']]<- glm(participation ~ 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,
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 participation')| SLR | MLR | Â Interaction | |
|---|---|---|---|
| (Intercept) | 0.2 (0.0)*** | 0.3 (0.0)*** | 0.3 (0.0)*** |
| Extreme Left | 2.0 (0.2)*** | 2.0 (0.2)*** | 2.5 (0.4)*** |
| Extreme Right | 1.7 (0.2)*** | 1.7 (0.2)*** | 1.9 (0.4)** |
| Moderate Left | 1.6 (0.1)*** | 1.5 (0.1)*** | 1.5 (0.2)** |
| Moderate Right | 1.0 (0.1) | 1.0 (0.1) | 0.8 (0.1) |
| Later Adulthood | 0.5 (0.0)*** | 0.5 (0.1)*** | |
| Extreme Left:Later Adulthood | 0.8 (0.1) | ||
| Extreme Right:Later Adulthood | 0.9 (0.2) | ||
| Moderate Left:Later Adulthood | 1.0 (0.2) | ||
| Moderate Right:Later Adulthood | 1.3 (0.2) | ||
| Num.Obs. | 5605 | 5476 | 5476 |
| AIC | 10569.6 | 10225.9 | 10225.7 |
| BIC | 10602.8 | 10265.5 | 10291.7 |
| Log.Lik. | −5279.816 | −5106.945 | −5102.830 |
| RMSE | 0.38 | 0.38 | 0.38 |
output = "Table 2"model4 <- glm(participation ~ lrscale + Adulthood + lrscale*Adulthood, family = binomial, data = df, weights = round(weight))
mydf <- predict_response(model4, terms = c("lrscale", "Adulthood"))
print(mydf, collapse_tables = TRUE)## # Predicted probabilities of participation
##
## lrscale | Adulthood | Predicted | 95% CI
## --------------------------------------------
## 0 | Later | 0.17 | 0.16, 0.18
## 2 | | 0.17 | 0.16, 0.18
## 4 | | 0.16 | 0.16, 0.17
## 6 | | 0.16 | 0.15, 0.17
## 8 | | 0.16 | 0.15, 0.16
## 88 | | 0.06 | 0.05, 0.07
## 0 | Early | 0.27 | 0.25, 0.30
## 2 | | 0.27 | 0.25, 0.29
## 4 | | 0.26 | 0.24, 0.28
## 6 | | 0.26 | 0.24, 0.28
## 8 | | 0.25 | 0.23, 0.27
## 88 | | 0.11 | 0.09, 0.14##
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
## all rows.ggplot(mydf, aes(x, predicted, colour = group)) + geom_line()
plot_comparisons(model3, variables = "Ideology", condition = "Adulthood")
Fig1 <- ggplot(mydf, aes(x, predicted, colour = group)) + geom_line()
ggsave(filename = "Fig3.pdf", plot = Fig1, device = "pdf", width = 10, height = 8)df <- df %>%
mutate(behave = ipbhprp,
ideas = ipcrtiv, # Updated variable name
safety = ipstrgv,
well_being = iphlppl, # Updated variable name and ensure valid R syntax
rules = ipfrule) %>%
mutate(across(c(behave, ideas, safety, well_being, rules),
~ na_if(.x, 7) %>% na_if(8) %>% na_if(9))) %>%
# Apply the reverse coding
mutate(across(c(behave, ideas, safety, well_being, rules), ~ 7 - .x))
# Now you can calculate 'schwartzauth' after the NA recoding
df$auth <- scales::rescale(df$behave +
df$ideas +
df$safety +
df$well_being +
df$rules, to=c(0,100), na.rm=TRUE)
df <- df %>% filter(!is.na(auth))
df <- df %>%
mutate(
# Recoding generation
gen = case_when(
yrbrn %in% 1900:1945 ~ "1",
yrbrn %in% 1946:1964 ~ "2",
yrbrn %in% 1965:1979 ~ "3",
yrbrn %in% 1980:1996 ~ "4",
TRUE ~ as.character(yrbrn)
),
gen = factor(gen, levels = c("1", "2", "3", "4"),
labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials")),
religious = case_when(
rlgdgr %in% c(77, 88, 99) ~ NA_real_,
TRUE ~ rlgdgr
)
)
model5 <- lm(auth ~ gen + religious + gen*religious, data = df, weights = weight)check_model(model5)