# List of packages
packages <- c("tidyverse", "broom", "infer", "fst", "modelsummary", "effects", "survey", "interactions") # 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 'interactions' 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] "broom"     "lubridate" "forcats"   "stringr"   "dplyr"     "purrr"    
##  [7] "readr"     "tidyr"     "tibble"    "ggplot2"   "tidyverse" "stats"    
## [13] "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[3]]
##  [1] "infer"     "broom"     "lubridate" "forcats"   "stringr"   "dplyr"    
##  [7] "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"   "tidyverse"
## [13] "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"  
## [19] "base"     
## 
## [[4]]
##  [1] "fst"       "infer"     "broom"     "lubridate" "forcats"   "stringr"  
##  [7] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [13] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [19] "methods"   "base"     
## 
## [[5]]
##  [1] "modelsummary" "fst"          "infer"        "broom"        "lubridate"   
##  [6] "forcats"      "stringr"      "dplyr"        "purrr"        "readr"       
## [11] "tidyr"        "tibble"       "ggplot2"      "tidyverse"    "stats"       
## [16] "graphics"     "grDevices"    "utils"        "datasets"     "methods"     
## [21] "base"        
## 
## [[6]]
##  [1] "effects"      "carData"      "modelsummary" "fst"          "infer"       
##  [6] "broom"        "lubridate"    "forcats"      "stringr"      "dplyr"       
## [11] "purrr"        "readr"        "tidyr"        "tibble"       "ggplot2"     
## [16] "tidyverse"    "stats"        "graphics"     "grDevices"    "utils"       
## [21] "datasets"     "methods"      "base"        
## 
## [[7]]
##  [1] "survey"       "survival"     "Matrix"       "grid"         "effects"     
##  [6] "carData"      "modelsummary" "fst"          "infer"        "broom"       
## [11] "lubridate"    "forcats"      "stringr"      "dplyr"        "purrr"       
## [16] "readr"        "tidyr"        "tibble"       "ggplot2"      "tidyverse"   
## [21] "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
## [26] "methods"      "base"        
## 
## [[8]]
##  [1] "interactions" "survey"       "survival"     "Matrix"       "grid"        
##  [6] "effects"      "carData"      "modelsummary" "fst"          "infer"       
## [11] "broom"        "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?

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

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

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

sweden_data <- sweden_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_                           
    )
  )
table(sweden_data$ID)
## 
##  Left Right 
##  6188  7406
sweden_pop_model <- lm(populist ~ ID, data = sweden_data)
coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)
## (Intercept)     IDRight 
##  48.3188406  -0.7149274

According to the histogram, Sweden is said to showcase right-wing populist.

Intercept (67.4):

IDRight Coefficient (-0.71):

Task 2

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.

sweden_data <- sweden_data %>%
  mutate(
    educ.ba = case_when(
      essround < 5 & edulvla == 5 ~ "BA or more",
      essround >= 5 & edulvlb > 600 ~ "BA or more",
      TRUE ~ "No BA"                         
    )

  )
table(sweden_data$ID)
## 
##  Left Right 
##  6188  7406
tidy(sweden_pop_model) |>
knitr::kable(digits = 3)
term estimate std.error statistic p.value
(Intercept) 48.319 0.265 182.562 0.000
IDRight -0.715 0.356 -2.011 0.044 
summary(sweden_pop_model)
## 
## Call:
## lm(formula = populist ~ ID, data = sweden_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.319 -14.271  -1.652  12.396  52.396 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  48.3188     0.2647 182.562   <2e-16 ***
## IDRight      -0.7149     0.3555  -2.011   0.0444 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.25 on 11864 degrees of freedom
##   (6350 observations deleted due to missingness)
## Multiple R-squared:  0.0003407,  Adjusted R-squared:  0.0002565 
## F-statistic: 4.044 on 1 and 11864 DF,  p-value: 0.04436

This value of intercept expected the average populist scale score would be 48.32 for the left, since it is the reference category. This higher scores suggest stronger populist attitudes. The negative coefficient associated with the “Right” group indicates that, on average, individuals who identify as ideologically right in Sweden have populist scale scores that are 0.71 points lower than those who identify as left. The p-value is 0.

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.

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))) %>%
  mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))

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

Based on the graph, the data suggest that the measure of authoritarian attitudes has remained relatively stable over the years, with no significant increases or decreases. Indicates that the average level of authoritarian attitudes is above the midpoint in 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.

df <- read.fst("italy_data.fst")
model_data <- df %>%
  mutate(
    gndr = case_when(
      gndr == 1 ~ "Male",
      gndr == 2 ~ "Female",
      gndr == 9 ~ NA_character_,
      TRUE ~ as.character(gndr)
    ),

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

    age = ifelse(agea == 999, NA, agea),

    cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),

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

    gen = factor(gen,
                 levels = c("1", "2", "3", "4"),
                 labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
  )
table(model_data$gen)
## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         1470         3179         2636         2031
models <- list(
   "Model 1" = lm(auth ~ gen, data = model_data),
   "Model 2" = lm(auth ~ gndr, data = model_data),
   "Model 3" = lm(auth ~ educ.ba, data = model_data))
modelsummary(models, fmt = 1,
  estimate  = c("{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_rename = c("genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials", "gndrMale" = "Male", "educ.baBA or more" = "BA or More"),
  title = 'Table 2. Regression models for authoritarian attitudes')
Table 2. Regression models for authoritarian attitudes
Model 1  Model 2  Model 3
(Intercept) 75.3 (0.5)*** 74.4 (0.2)*** 73.3 (0.2)***
Baby Boomers −1.1 (0.6)*
Gen X −2.1 (0.6)***
Millennials −3.6 (0.6)***
Male −2.6 (0.3)***
BA or More −1.1 (0.4)*
Num.Obs. 7631 8432 8432
R2 0.006 0.007 0.001
R2 Adj. 0.006 0.007 0.001
AIC 62934.0 69523.1 69579.5
BIC 62968.7 69544.3 69600.6
Log.Lik. −31462.009 −34758.566 −34786.739
RMSE 14.94 14.93 14.98

Coefficients Interpretation:

Baby Boomers exhibit scores that are, on average, 1.1 points lower on the authoritarian values scale. Generation X shows an even more pronounced decrease, with scores 2.1 points lower than the Interwar generation on the authoritarian values scale. Millennials exhibit the most considerable decrease in authoritarian values, with their scores being 3.6 points lower than those of the Interwar generation. Millennials diverge most strongly from the Interwar generation in comparing their authoritarian values. The adjusted R-squared value of 0.006.

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.

greece_data <- read.fst("greece_data.fst")
model_data <- df %>%
  mutate(
    gndr = case_when(
      gndr == 1 ~ "Male",
      gndr == 2 ~ "Female",
      gndr == 9 ~ NA_character_,
      TRUE ~ as.character(gndr)
    ),

    ID = case_when(
      lrscale %in% 0:4 ~ "Left",                     
      lrscale %in% 6:10 ~ "Right",                   
      lrscale %in% c(77, 88, 99) ~ NA_character_,    
      TRUE ~ NA_character_                           
    ),
    
    age = ifelse(agea == 999, NA, agea),

    cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),

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

    gen = factor(gen,
                 levels = c("1", "2", "3", "4"),
                 labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
  )
table(model_data$gen)
## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         1470         3179         2636         2031
models <- list(
   "Model 1" = lm(auth ~ lrscale, data = model_data),
   "Model 2" = lm(auth ~ gndr, data = model_data),
   "Model 3" = lm(auth ~ 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) 73.4 (0.2)*** 74.4 (0.2)*** 75.3 (0.5)***
lrscale 0.0 (0.0)+
gndrMale −2.6 (0.3)***
genBaby Boomers −1.1 (0.6)*
genGen X −2.1 (0.6)***
genMillennials −3.6 (0.6)***
Num.Obs. 8432 8432 7631
R2 0.000 0.007 0.006
R2 Adj. 0.000 0.007 0.006
AIC 69581.4 69523.1 62934.0
BIC 69602.5 69544.3 62968.7
Log.Lik. −34787.677 −34758.566 −31462.009
RMSE 14.98 14.93 14.94 
modelsummary(
  (models),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_omit = "Intercept")
Model 1  Model 2  Model 3
lrscale 0.0 (0.0)+
gndrMale −2.6 (0.3)***
genBaby Boomers −1.1 (0.6)*
genGen X −2.1 (0.6)***
genMillennials −3.6 (0.6)***
Num.Obs. 8432 8432 7631
R2 0.000 0.007 0.006
R2 Adj. 0.000 0.007 0.006
AIC 69581.4 69523.1 62934.0
BIC 69602.5 69544.3 62968.7
Log.Lik. −34787.677 −34758.566 −31462.009
RMSE 14.98 14.93 14.94 
modelplot(models, coef_omit = 'Interc')

Based on the data, the negative coefficients for generational variables suggest that these cohorts exhibit lower levels of the outcome variable compared to the reference group. Also suggests that right-wing political orientation is associated with a decrease in the dependent variable in the specific context of Model 1.