# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "interactions","broom") # 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
##
##
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools
## [[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"
##
## [[8]]
## [1] "broom" "interactions" "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"
Graph the histogram distribution using the populist scale (as we did in the tutorial) for the country of Sweden. What do you note?
sweden_data <- read.fst("sweden_data.fst")
# Populist scale (Norris/Inglehart use low levels of trust in parliaments, politicians and parties as an indicator of populism)
# Setting values greater than 10 to NA
sweden_data$trstplt <- ifelse(sweden_data$trstplt > 10, NA, sweden_data$trstplt)
sweden_data$trstprl <- ifelse(sweden_data$trstprl > 10, NA, sweden_data$trstprl)
sweden_data$trstprt <- ifelse(sweden_data$trstprt > 10, NA, sweden_data$trstprt)
# Creating and rescaling the trust variable
sweden_data$trust <- scales::rescale(sweden_data$trstplt + sweden_data$trstprl + sweden_data$trstprt, na.rm = TRUE, to = c(0, 100))
# "flip it" - populist is conceived by N + I as 'anti-trust'
sweden_data$populist <- scales::rescale(sweden_data$trust, na.rm = TRUE, to=c(100,0))
ggplot(sweden_data, aes(x = populist)) +
geom_histogram(bins = 30, fill = "blue", color = "black") +
theme_minimal() +
labs(title = "Distribution of Populist Scale for Sweden",
x = "Populist Scale",
y = "Count")
## Warning: Removed 2503 rows containing non-finite values (`stat_bin()`).
There is a peak around 50 on the populist scale which signifies a moderate level of populism in Sweden.
Run a linear regression with populist attitudes as the outcome, educational attainment (recoded as a binary “BA or more” and “No BA”), and using the Swedish data. Print the intercept and coefficients only and interpret. Then do a tidy model table displaying the p-value (with digits = 3) and interpret.
model.s_data <- sweden_data %>%
mutate(
# Recoding education
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
),
# Handling NAs for education levels
edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
# Explicitly making 'No BA' the reference category
educ.ba = factor(educ.ba, levels = c("No BA", "BA or more")),
)
sweden_model <- lm(populist ~ educ.ba, data = model.s_data)
tidy(sweden_model) |>
knitr::kable(digits = 3)
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 51.359 | 0.181 | 283.363 | 0 |
| educ.baBA or more | -7.868 | 0.350 | -22.494 | 0 |
Intercept (51.359) Interpretation: This value represents the estimated average score on the populist attitude scale for individuals who do not have a BA degree (reference group). In other words, individuals who do not have a BA degree have an estimated average populist attitude score of approximately 51.359.
Coefficient for “BA or more” (-7.868) Interpretation: The negative coefficient indicates that, on average, individuals with a BA degree or higher tend to have lower populist attitude scores compared to those without a BA degree. Specifically, individuals with a BA degree or higher have an estimated average populist attitude score that is 7.868 points lower than individuals without a BA degree.
P-value (0) Interpretation: The p-value indicates it is statistically significant. Therefore, we should reject the null hypothesis and conclude that there is a statistically significant relationship between having a BA degree or higher and lower populist attitude scores.
Now using the Italian dataset and the authoritarian values scale (as we did in the tutorial), graph the average by survey year. Interpret.
df <- read.fst("italy_data.fst")
df <- df %>%
mutate(behave = ipbhprp,
secure = impsafe,
safety = ipstrgv,
tradition = imptrad,
rules = ipfrule) %>%
mutate(across(c("behave", "secure", "safety", "tradition", "rules"),
~ na_if(.x, 7) %>% na_if(8) %>% na_if(9))) %>%
# Apply the reverse coding
mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))
df$auth <- scales::rescale(df$behave +
df$secure +
df$safety +
df$tradition +
df$rules, to=c(0,100), na.rm=TRUE)
df <- df %>% filter(!is.na(auth))
table(df$secure)
##
## 1 2 3 4 5 6
## 54 139 626 1856 2876 2881
table(df$auth)
##
## 0 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84
## 7 4 4 9 13 20 27 57 113 155 163 323 367 596 647 789 823 897 992 733
## 88 92 96 100
## 632 465 299 297
ggplot(df, aes(x = auth)) +
geom_histogram(bins = 30, fill = "blue", color = "black") +
theme_minimal() +
labs(title = "Distribution of Authoritarian Scale for Italy",
x = "Authoritarian Scale",
y = "Count")
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]
}
auth_avg <- df %>%
group_by(year) %>%
summarize(auth_avg = mean(auth, na.rm = TRUE))
plot_auth <- ggplot(auth_avg, aes(x = year, y = auth_avg)) +
geom_point(aes(color = auth_avg), alpha = 1.0) +
geom_line(aes(group = 1), color = "blue", linetype = "dashed") +
labs(title = "Authoritarian Scale Average by Year (Italy)",
x = "Survey Year",
y = "Authoritarian Scale Average") +
theme_minimal() +
theme(legend.position = "none") +
scale_y_continuous(limits = c(0, 100))
print(plot_auth)
The average authoritarian values scale in Italy has remained relatively constant at around 75. This suggests that there has not been a significant shift in authoritarian attitudes during this period.
Now run a linear regression model, using the modelsummary package, with authoritarian values as the outcome, and generations as the predictor/potential explanatory variable, and using the Italian data again. Rename the coefficients using the coef_rename function. Interpret the coefficients, whether they are statistically significant, and the adjusted R-squared.
model.i_data <- df %>%
mutate(
# Recoding age, setting 999 to NA
age = ifelse(agea == 999, NA, agea),
# Recoding cohort variable, setting years before 1930 and after 2000 to NA
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
# Recoding generational cohorts based on year of birth
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(cohort) # Keeping other values as character if they do not fit the ranges
),
# Converting cohort to a factor with labels
gen = factor(gen,
levels = c("1", "2", "3", "4"),
labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
)
table(model.i_data$gen)
##
## Interwar Baby Boomers Gen X Millennials
## 1026 2580 2214 1811
italy_model <- lm(auth ~ gen, data = model.i_data)
modelsummary(
list(italy_model),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_omit = "Intercept")
| (1) | |
|---|---|
| genBaby Boomers | −1.1 (0.6)* |
| genGen X | −2.1 (0.6)*** |
| genMillennials | −3.6 (0.6)*** |
| Num.Obs. | 7631 |
| R2 | 0.006 |
| R2 Adj. | 0.006 |
| AIC | 62934.0 |
| BIC | 62968.7 |
| Log.Lik. | −31462.009 |
| RMSE | 14.94 |
models <- list(
"Model 1" = lm(auth ~ gen, data = model.i_data))
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}"),
statistic = NULL)
| Model 1 | |
|---|---|
| (Intercept) | 75.3 (0.5)*** |
| genBaby Boomers | −1.1 (0.6)* |
| genGen X | −2.1 (0.6)*** |
| genMillennials | −3.6 (0.6)*** |
| Num.Obs. | 7631 |
| R2 | 0.006 |
| R2 Adj. | 0.006 |
| AIC | 62934.0 |
| BIC | 62968.7 |
| Log.Lik. | −31462.009 |
| RMSE | 14.94 |
Coefficients Interpretation:
genBaby Boomers (-1.1, p<0.05): Compared to the Interwar generation (reference category that is omitted), Baby Boomers exhibit scores that are, on average, 1.1 points lower on the authoritarian values scale (0-100). The single asterisk (*) denotes statistical significance at the common alpha threshold of 0.05.
genGen X (-2.1, p<0.001): Generation X shows an even more pronounced decrease, with scores 2.1 points lower than the Interwar generation on the authoritarian values scale. The triple asterisks (***) denote statistical significance at a very low p-value (under our alpha threshold of 0.05).
genMillennials (-3.6, p<0.001): Millennials exhibit the most considerable decrease in authoritarian values, with their scores being 3.6 points lower than those of the Interwar generation. This marked difference is statistically significant, indicating that Millennials diverge most strongly from the Interwar generation in comparing their authoritarian values.
Adjusted R-Squared (Adj. R² = 0.006):
The adjusted R-squared value of 0.006 means that the ‘gen’ variable (with the four coded generations) explain approximately 0.6% of the variance in the authoritarian values scale for Italy, after accounting for the predictors. While statistically significant, this suggests that the model explains only a small portion of the variability in authoritarian values, pointing to other factors outside generational differences that may play critical roles in shaping these attitudes.
modelplot(models, coef_omit = 'Interc')
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_rename = c("genMillennials" = "Millennials", "genGen X" = "Gen X", "genBaby Boomers" = "Baby Boomers"),
title = 'Table 2. Regression model for authoritarian attitudes')
| Model 1 | |
|---|---|
| (Intercept) | 75.3 (0.5)*** |
| Baby Boomers | −1.1 (0.6)* |
| Gen X | −2.1 (0.6)*** |
| Millennials | −3.6 (0.6)*** |
| Num.Obs. | 7631 |
| R2 | 0.006 |
| R2 Adj. | 0.006 |
| AIC | 62934.0 |
| BIC | 62968.7 |
| Log.Lik. | −31462.009 |
| RMSE | 14.94 |
Now visualize, using the modelplot() function from the modelsummary package, coefficient estimates and 95% confidence intervals for the following three models:
Model 1 = populist attitudes scale as the outcome; left/right as the single predictor (omitting moderates as we did in the tutorial), and using data for Greece.
Model 2 = populist attitudes scale as the outcome; gender as the single predictor, and using data for Greece.
Model 3 = populist attitudes scale as the outcome; generations as the single predictor (using the four generational category recode we did in the tutorial), and using data for Greece. Interpret.
greece_data <- read.fst("greece_data.fst")
# Populist scale (Norris/Inglehart use low levels of trust in parliaments, politicians and parties as an indicator of populism)
# Setting values greater than 10 to NA
greece_data$trstplt <- ifelse(greece_data$trstplt > 10, NA, greece_data$trstplt)
greece_data$trstprl <- ifelse(greece_data$trstprl > 10, NA, greece_data$trstprl)
greece_data$trstprt <- ifelse(greece_data$trstprt > 10, NA, greece_data$trstprt)
# Creating and rescaling the trust variable
greece_data$trust <- scales::rescale(greece_data$trstplt + greece_data$trstprl + greece_data$trstprt, na.rm = TRUE, to = c(0, 100))
# "flip it" - populist is conceived by N + I as 'anti-trust'
greece_data$populist <- scales::rescale(greece_data$trust, na.rm = TRUE, to=c(100,0))
model.g_data <- greece_data %>%
mutate(
ID = case_when(
lrscale %in% 0:4 ~ "Left",
lrscale %in% 6:10 ~ "Right",
lrscale %in% c(77, 88, 99) ~ NA_character_,
TRUE ~ NA_character_
),
# Recoding gender
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
# Recoding age, setting 999 to NA
age = ifelse(agea == 999, NA, agea),
# Recoding cohort variable, setting years before 1930 and after 2000 to NA
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
# Recoding generational cohorts based on year of birth
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(cohort) # Keeping other values as character if they do not fit the ranges
),
# Converting cohort to a factor with labels
gen = factor(gen,
levels = c("1", "2", "3", "4"),
labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
)
table(model.g_data$gen)
##
## Interwar Baby Boomers Gen X Millennials
## 2978 3484 3498 2330
model1 <- lm(populist ~ lrscale, data = model.g_data)
model2 <- lm(populist ~ gndr, data = model.g_data)
model3 <- lm(populist ~ gen, data = model.g_data)
modelsummary(
list(model1),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_omit = "Intercept")
| (1) | |
|---|---|
| lrscale | 0.1 (0.0)*** |
| Num.Obs. | 9823 |
| R2 | 0.021 |
| R2 Adj. | 0.021 |
| AIC | 88667.6 |
| BIC | 88689.2 |
| Log.Lik. | −44330.802 |
| RMSE | 22.07 |
models <- list(
"Model 1" = lm(populist ~ lrscale, data = model.g_data),
"Model 2" = lm(populist ~ gndr, data = model.g_data),
"Model 3" = lm(populist ~ gen, data = model.g_data))
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL)
| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 67.5 (0.3)*** | 70.0 (0.3)*** | 65.8 (0.5)*** |
| lrscale | 0.1 (0.0)*** | ||
| gndrMale | −1.0 (0.5)* | ||
| genBaby Boomers | 3.5 (0.7)*** | ||
| genGen X | 5.6 (0.7)*** | ||
| genMillennials | 5.8 (0.7)*** | ||
| Num.Obs. | 9823 | 9823 | 9562 |
| R2 | 0.021 | 0.000 | 0.009 |
| R2 Adj. | 0.021 | 0.000 | 0.009 |
| AIC | 88667.6 | 88873.2 | 86473.9 |
| BIC | 88689.2 | 88894.8 | 86509.7 |
| Log.Lik. | −44330.802 | −44433.590 | −43231.944 |
| RMSE | 22.07 | 22.30 | 22.25 |
modelplot(models, coef_omit = 'Interc')
Model 1 Interpretation:
Intercept (67.5): This value represents the average score on the populist scale for the reference group in Greece, which is the “Left” group. Therefore, the expected average populist scale score would be 67.5 for the left, since it is the reference category. This score indicates a certain level of populism, where higher scores suggest stronger populist attitudes.
IDRight Coefficient (0.1, p < 0.001): The positive coefficient associated with the “Right” group indicates that, on average, individuals who identify as ideologically right in Greece have populist scale scores that are 0.1 points higher than those who identify as left.
-Comparing right-leaning individuals to left-leaning ones in Greece, we expect to see the average populist scale score to be 0.1 points higher for the right-leaning group. This effect size is small (0.1 points) and may not necessarily indicate a substantial difference in populist attitudes between the two groups.
Model 2 Interpretation:
Intercept (70.0): The intercept of 70.0 represents the average value of populist attitudes for individuals who identify as female in Greece.
Male Gender Coefficient (-1.0, p<0.05): The coefficient of -1.0 implies that, compared to females, males tend to have an average decrease in populist attitudes by 1.0 points.
-Comparing males to females in Greece, we expect to see the average populist scale score to be 0.1 points lower for men. This effect size is relatively small (-1.0 units) and may not necessarily indicate a substantial difference in populist attitudes between the two genders.
Model 3 Interpretation:
Intercept (65.8): The intercept value of 65.8 represents the average score on the populist attitudes scale for the reference group (Interwar generation).
Baby Boomers, Gen X, and Millennials: The coefficient estimates for each of these generational categories (3.5, 5.6, and 5.8; all p < 0.001) indicate that, on average, each generation has higher scores on the populist attitudes scale compared to the reference group (Interwar generation). This suggests that Baby Boomers, Gen X, and Millennials, exhibit stronger populist attitudes compared to the Interwar generation.
-Comparing Baby Boomers, Gen X, and Millennials to the Interwar generation in Greece, we expect to see differences of 3.5, 5.6, and 5.8 units, in the average populist scale score. These effect sizes are relatively small and may not necessarily indicate substantial differences in populist attitudes between the different generations.