Setting up environment

# 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"

Task 1

Graph the histogram distribution using the populist scale (as we did in the tutorial) for the country of Sweden. What do you note?

Reading in Sweden Data

sweden_data <- read.fst("sweden_data.fst")

Creating the Populist Scale for Sweden Data

# 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))

Graphing Histogram Distribution of Sweden Data

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()`).

Interpretation

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.

Task 2

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.

Re-coding ‘educ.ba’

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
  )

Double-checking if ‘educ.ba’ was cleaned

unique(sweden_data_lr$educ.ba)
## [1] "No BA"      "BA or more"

Creating linear regression model for Sweden

sweden_pop_model <- lm(populist ~ educ.ba, data = sweden_data_lr) # populism = dependent (outcome), educ.ba = independent (explanatory)

Printing intercept and coefficient

coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)
##  (Intercept) educ.baNo BA 
##    43.490559     7.868121

Interpretation

Intercept (43.5):

  • This value represents the average score on the populist scale for the reference group in Sweden, which is the ‘having a BA’ group. The expected average populist scale score would be 43.5 for the people with a BA, since it is the reference category. Scores higher than 43.5 suggest stronger populist attitudes.

educ.baNo BA Coefficient (7.9):

  • The positive coefficient associated with the “No BA” group indicates that, on average, individuals who have no BA in Sweden have populist scale scores that are 7.9 points higher than those who have a BA.

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

Interpretation

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.

Task 3

Now using the Italian dataset and the authoritarian values scale (as we did in the tutorial), graph the average by survey year. Interpret.

Coding for authoritarian scale

italy_data <- read.fst("italy_data.fst")
## Creates an authoritarian values scale based on human modules items (can also do libertarian value scale)
italy_data <- italy_data |>
  mutate(behave = ipbhprp,    # renaming columns for organization
         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 ))

# Now you can calculate 'schwartzauth' after the NA recoding - creating new column for authoritarian values and adding values in columns, scaling values to 0 and 100
italy_data$auth <- scales::rescale(italy_data$behave + 
                      italy_data$secure + 
                      italy_data$safety + 
                      italy_data$tradition + 
                      italy_data$rules, to=c(0,100), na.rm=TRUE)


italy_data <- italy_data %>% filter(!is.na(auth)) #any rows in whole dataframe where auth = NA value - filter it
table(italy_data$secure)
## 
##    1    2    3    4    5    6 
##   54  139  626 1856 2876 2881

Double-checking rescaling

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

Creating new column for year

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)

Interpretation

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.

Task 4

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.

Creating generation column

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

Running linear regression with ‘auth’ and ‘gen’

model_italy <- lm(auth ~ gen, data = model_data_italy)

Modelsummary

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

Fixed modelsummary

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 with re-named coefficients

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')
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

Interpretation of co-efficients and r-squared

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.

Task 5

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.

Reading in Greece Data

greece_data <- read.fst("greece_data.fst")

Coding Populist Attitudes Scale

# 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))

Cleaning and Re-coding ‘gndr’, ‘lrscale’ and ‘gen’ variables

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"))
  )

Double-checking variables

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

Cleaning NA’s

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

Creating linear regressions

model1 <- lm(populist ~ ID, data = model_data)
model2 <- lm(populist ~ gndr, data = model_data)
model3 <- lm(populist ~ gen, data = model_data)

Grouping models

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

Visualizing Coefficient Estimates and 95% Confidence Intervals

modelplot(models, coef_omit = 'Interc') # omit intercept

Interpretation

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.