# 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)
# Note: Schafer does it 2 ways, we will focus now on the most straightforward
# 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 re-scaling the trust variable - taking the 3 columns, adding the values up and scaling results to be between 0 and 100, if someone answered 10, 10, 10 = 30 but scaling it so that maximum value = 100
sweden_data$trust <- scales::rescale(sweden_data$trstplt + sweden_data$trstprl + sweden_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' - don't trust establishment, etc. - inversion of trust score is populism score - create new column, take what we just did for trust column and invert the scale so minimum = 100, max = 0
# As Schafer explains there are some issues with that conceptualization, but N + I is also widely applied
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 = "purple", 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()`).
The distribution is somewhat normal and uni-modal, with a peak of respondents who scored around 50 or less for populism, indicating that Swedish respondents have either lower populist scores. This is further supported by the distribution tapering off as the populist scores increase, suggesting Swedish respondents have lower populist beliefs.
Run a linear regression with populist attitudes as the outcome and educational attainment (recoded as a binary “BA or more” and “No BA”) 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.
sweden_data_lr <- sweden_data |>
mutate(
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
) ) |>
# Cleaning NA's
mutate(
edulvla = ifelse(edulvla %in% c(77, 88), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba) # Force educ.ba to be treated as character
)
unique(sweden_data_lr$educ.ba)
## [1] "No BA" "BA or more"
sweden_pop_model <- lm(populist ~ educ.ba, data = sweden_data_lr) # populism = dependent (outcome), educ.ba = independent (explanatory)
coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)
## (Intercept) educ.baNo BA
## 43.490559 7.868121
Intercept (43.5):
educ.baNo BA Coefficient (7.9):
CONCLUSION: Comparing individuals with no BA to individuals with BAs in Sweden, we expect to see the average populist scale score to be 7.9 points higher for the ‘no BA’ group. This suggests that individuals with no BAs exhibit stronger populist attitudes than those with a BA.
tidy(sweden_pop_model) |>
knitr::kable(digits = 3)
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 43.491 | 0.299 | 145.376 | 0 |
| educ.baNo BA | 7.868 | 0.350 | 22.494 | 0 |
We see the intercept and coefficient with the ‘having a BA’ category as the reference category. Beside the coefficient is the standard error and beside that are the observed t-statistic results. A t-test was used because the populist scale (response) is a numeric variable and educational attainment (explanatory) is a categorical variable with two categories. The p-value of almost 0 indicates that the coefficient estimate is statistically significant, as it is below an alpha of 0.05. As a result, having no BA is associated with an increase in populist attitudes in Sweden.
Now using the Italian dataset and the authoritarian values scale (as we did in the tutorial), graph the average by survey year. Interpret.
table(italy_data$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
italy_data$year <- NA # creating new column for years from ESS rounds
replacements <- c(2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020)
for(i in 1:10){
italy_data$year[italy_data$essround == i] <- replacements[i]# putting in replacements one at time
}
auth_avg <- italy_data %>% # want to find average authoritarian score by year in every column
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) + # made colour change with value - lower values = darker, lighter values = higher
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)) # showing the whole scale of 0-100 bc by default it zooms into where the data is
print(plot_auth)
The average authoritarian score in Italy is highest in 2012 at ~75. Since 2012, the average authoritarian score for Italy has decreased slightly until 2016 when it became relatively stable at ~73-74. Overall, the average authoritarian score is high for Italy.
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_data_italy <- italy_data |>
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, | = "or"
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
# Recoding generational cohorts based on year of birth - creating new column called generation - re-coding to character type even though he's re-coding it as numbers
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 - converting character to factor - can give it levels/labels
gen = factor(gen,
levels = c("1", "2", "3", "4"), # ordinal variable of years being associated with generations
labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
)
table(model_data_italy$gen)
##
## Interwar Baby Boomers Gen X Millennials
## 1026 2580 2214 1811
model_italy <- lm(auth ~ gen, data = model_data_italy)
modelsummary(
list(model_italy),
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_data_italy)) # "Model 1" = label
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 |
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_rename = c("genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials"),
title = 'Table 2. Regression models 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 |
Baby Boomers (-1.1, p<0.05): Compared to the Interwar generation, Baby Boomers exhibit authoritarian scores that are, on average, 1.1 points lower on the authoritarian values scale (0-100). The single asterisk (*) denotes statistical significance below the p-value of 0.05, which is under our alpha threshold of 0.05 meaning the coefficient is statistically significant.
Gen X (-2.1, p<0.001): Generation X shows a slightly more pronounced decrease, with authoritarian scores 2.1 points lower than the Interwar generation on the authoritarian values scale. The triple asterisks (***) denote stronger statistical significance at a very low p-value under our alpha threshold of 0.05, meaning the coefficient is also statistically significant.
Millennials (-3.6, p<0.001): Millennials exhibit the most decrease in authoritarian values, with their scores being 3.6 points lower than those of the Interwar generation. This marked difference is statistically significant as it is below an alpha of 0.001, 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 explains approximately 0.6% of the variance in the authoritarian values scale for Italy, after accounting for the number of 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.
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)
# Note: Schafer does it 2 ways, we will focus now on the most straightforward
# 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 re-scaling the trust variable - taking the 3 columns, adding the values up and scaling results to be between 0 and 100, if someone answered 10, 10, 10 = 30 but scaling it so that maximum value = 100
greece_data$trust <- scales::rescale(greece_data$trstplt + greece_data$trstprl + greece_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' - don't trust establishment, etc. - inversion of trust score is populism score - create new column, take what we just did for trust column and invert the scale so minimum = 100, max = 0
# As Schafer explains there are some issues with that conceptualization, but N + I is also widely applied
greece_data$populist <- scales::rescale(greece_data$trust, na.rm = TRUE, to=c(100,0))
model_data <- greece_data %>%
mutate(
# Recoding gender
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
ID = case_when( # creating new column 'ID' = 'lrscale'
lrscale %in% 0:4 ~ "Left", # This time: Values 0 to 4 in 'lrscale' = "Left"
lrscale %in% 6:10 ~ "Right", # This time: Values 6 to 10 in 'lrscale'= "Right"
lrscale %in% c(77, 88, 99) ~ NA_character_, # Values 77, 88, 99 in 'lrscale' = NA - 5 bc can't be both left and right
TRUE ~ NA_character_ # Any other values not explicitly matched above
),
# Recoding age, setting 999 to NA
age = ifelse(agea == 999, NA, agea),
# Recoding cohort variable, setting years before 1930 and after 2000 to NA, | = "or"
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
# Recoding generational cohorts based on year of birth - creating new column called generation - re-coding to character type even though he's re-coding it as numbers
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 - converting character to factor - can give it levels/labels
gen = factor(gen,
levels = c("1", "2", "3", "4"), # ordinal variable of years being associated with generations
labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
)
table(model_data$gndr)
##
## Female Male
## 6929 5629
table(model_data$ID)
##
## Left Right
## 2688 3522
table(model_data$gen)
##
## Interwar Baby Boomers Gen X Millennials
## 2978 3484 3498 2330
unique(model_data$gndr)
## [1] "Male" "Female"
unique(model_data$ID)
## [1] "Right" NA "Left"
unique(model_data$gen)
## [1] Interwar Gen X Baby Boomers Millennials <NA>
## Levels: Interwar Baby Boomers Gen X Millennials
model_data <- model_data |> filter(!is.na(ID))
model_data <- model_data |> filter(!is.na(gen))
unique(model_data$ID)
## [1] "Right" "Left"
unique(model_data$gen)
## [1] Interwar Gen X Baby Boomers Millennials
## Levels: Interwar Baby Boomers Gen X Millennials
model1 <- lm(populist ~ ID, data = model_data)
model2 <- lm(populist ~ gndr, data = model_data)
model3 <- lm(populist ~ gen, data = model_data)
models <- list(
"Model 1" = lm(populist ~ ID, data = model_data),
"Model 2" = lm(populist ~ gndr, data = model_data),
"Model 3" = lm(populist ~ gen, data = model_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) | 72.3 (0.5)*** | 67.9 (0.5)*** | 63.0 (0.7)*** |
| IDRight | −9.0 (0.7)*** | ||
| gndrMale | −1.3 (0.7)* | ||
| genBaby Boomers | 3.9 (0.9)*** | ||
| genGen X | 6.8 (0.9)*** | ||
| genMillennials | 6.4 (1.0)*** | ||
| Num.Obs. | 4822 | 4822 | 4822 |
| R2 | 0.038 | 0.001 | 0.013 |
| R2 Adj. | 0.038 | 0.001 | 0.012 |
| AIC | 43702.6 | 43885.7 | 43832.6 |
| BIC | 43722.0 | 43905.2 | 43865.0 |
| Log.Lik. | −21848.287 | −21939.862 | −21911.303 |
| RMSE | 22.47 | 22.90 | 22.76 |
modelplot(models, coef_omit = 'Interc') # omit intercept
genGen X are the most likely to be populist with populist scores 6.8 points above Interwar generations, similar to genMillennials with populist scores 6.4 points higher than Interwar generations, with their confidence intervals denoting statistical significance by being further away from 0. genBaby Boomers have lower populist scores compared to Millennials and Gen X with 3.9 points higher than the Interwar generation, suggesting that generations are becoming more populist as they progress. Furthermore, those who identify as right-winged in Greece have populist scores 9 points lower than those who identify as left-winged, suggesting that those who identify as right-winged have more trust and is statistically significant because its confidence interval is far from 0. Finally, men have populist scores 0.3 points less than women, which has lower statistical significance when compared to ‘ID’ and ‘gen’ due to its confidence interval being closer to 0.