rm(list=ls()); gc() 
##          used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 526709 28.2    1170174 62.5         NA   669417 35.8
## Vcells 970510  7.5    8388608 64.0      16384  1851678 14.2
# 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"
ess <- read_fst("All-ESS-Data.fst")

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

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

The height of the bar located at the range of value at around 45 of the Populist scale is significantly higher than the others, indicating that a greater number or proportion of observations within that range.

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.

spain_data <- read.fst("spain_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
spain_data$trstplt <- ifelse(spain_data$trstplt > 10, NA, spain_data$trstplt)
spain_data$trstprl <- ifelse(spain_data$trstprl > 10, NA, spain_data$trstprl)
spain_data$trstprt <- ifelse(spain_data$trstprt > 10, NA, spain_data$trstprt)

# Creating and rescaling the trust variable
spain_data$trust <- scales::rescale(spain_data$trstplt + spain_data$trstprl + spain_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
spain_data$populist <- scales::rescale(spain_data$trust, na.rm = TRUE, to=c(100,0))
ggplot(spain_data, aes(x = populist)) +
  geom_histogram(bins = 30, fill = "blue", color = "black") +
  theme_minimal() +
  labs(title = "Distribution of Populist Scale for Spain",
       x = "Populist Scale",
       y = "Count")
## Warning: Removed 2775 rows containing non-finite values (`stat_bin()`).

norway_data <- read.fst("norway_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
norway_data$trstplt <- ifelse(norway_data$trstplt > 10, NA, norway_data$trstplt)
norway_data$trstprl <- ifelse(norway_data$trstprl > 10, NA, norway_data$trstprl)
norway_data$trstprt <- ifelse(norway_data$trstprt > 10, NA, norway_data$trstprt)

# Creating and rescaling the trust variable
norway_data$trust <- scales::rescale(norway_data$trstplt + norway_data$trstprl + norway_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
norway_data$populist <- scales::rescale(norway_data$trust, na.rm = TRUE, to=c(100,0))
ggplot(norway_data, aes(x = populist)) +
  geom_histogram(bins = 30, fill = "blue", color = "black") +
  theme_minimal() +
  labs(title = "Distribution of Populist Scale for Norway",
       x = "Populist Scale",
       y = "Count")
## Warning: Removed 2198 rows containing non-finite values (`stat_bin()`).

sweden_data <- sweden_data %>%
  mutate(
    ID = case_when(
      lrscale %in% 0:4 ~ "Left",                     # This time: Values 0 to 3 in 'lrscale' are categorized as "Left"
      lrscale %in% 6:10 ~ "Right",                   # This time: Values 6 to 10 in 'lrscale' are categorized as "Right"
      lrscale %in% c(77, 88, 99) ~ NA_character_,    # Values 77, 88, 99 in 'lrscale' are set as NA
      TRUE ~ NA_character_                           # Any other values not explicitly matched above
    )
  )

spain_data <- spain_data %>%
  mutate(
    ID = case_when(
      lrscale %in% 0:4 ~ "Left",                     # Values 0 to 4 in 'lrscale' are categorized as "Left"
      lrscale %in% 6:10 ~ "Right",                   # Values 6 to 10 in 'lrscale' are categorized as "Right"
      lrscale %in% c(77, 88, 99) ~ NA_character_,    # Values 77, 88, 99 in 'lrscale' are set as NA
      TRUE ~ NA_character_                           # Any other values not explicitly matched above
    )
  )
table(sweden_data$ID)
## 
##  Left Right 
##  6188  7406
table(spain_data$ID)
## 
##  Left Right 
##  7618  4113
sweden_pop_model <- lm(populist ~ ID, data = sweden_data)
spain_pop_model <- lm(populist ~ ID, data = spain_data)
coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)
## (Intercept)     IDRight 
##  48.3188406  -0.7149274

Intercept (43.32): The value represents the average score on the populist scale for the reference group in Sweden, the expected average populist scale score would be 43.32. This score indicates a certain level of populism, where higher scores suggest stronger populist attitudes.

IDRight Coefficient (-0.71): The negative coeffcient indicates that, on average, Sweden have populist scale scires that are 0.71 points lower than those who identify as left.

tidy(sweden_pop_model)
## # A tibble: 2 × 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)   48.3       0.265    183.    0     
## 2 IDRight       -0.715     0.356     -2.01  0.0444
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

The linear regression analysis conducted on Sweden’s data indicates that being categorized as “Right” is associated with a statistically significant decrease of approximately 0.715 in the estimated mean value of the “populist” scale compared to the baseline category. However, the model explains only a minute proportion of the variation in the “populist” scale, suggesting that other factors not accounted for in the analysis may significantly influence populist tendencies 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.

italy_data <- read.fst("italy_data.fst")
df <- read.fst("italy_data.fst")
## Creates an authoritarian values scale based on human modules items (can also do libertarian value scale)
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 ))

# Now you can calculate 'schwartzauth' after the NA recoding
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 survey 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)

There is small fluctuation for the authoritarian value scale.

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.

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

 # Recoding education
    educ.ba = case_when(
      essround < 5 & edulvla == 5 ~ "BA or more",
      essround >= 5 & edulvlb > 600 ~ "BA or more",
      TRUE ~ "No BA"
    ),
    
    # Handle 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")),

    # 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_data$gen)
## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         1026         2580         2214         1811
model1 <- lm(auth ~ gen, data = model_data)
model2 <- lm(auth ~ gndr, data = model_data)
model3 <- lm(auth ~ educ.ba, data = model_data)
modelsummary(
  list(model1),
  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

Generational categories significantly impact the outcome variable, with Baby Boomers associated with a decrease of 1.1 units, Gen X with a decrease of 2.1 units, and Millennials with a decrease of 3.6 units compared to the reference category. All coefficients except Baby Boomers are highly statistically significant. However, the model explains only a negligible proportion of the outcome variable’s variance, indicating that factors beyond generational categories likely influence it.

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")
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)
Model 1  Model 2  Model 3
(Intercept) 75.3 (0.5)*** 74.4 (0.2)*** 73.3 (0.2)***
genBaby Boomers −1.1 (0.6)*
genGen X −2.1 (0.6)***
genMillennials −3.6 (0.6)***
gndrMale −2.6 (0.3)***
educ.baBA 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 
modelplot(models, coef_omit = 'Interc')

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

Across models, Baby Boomers, Gen X, and Millennials all exhibit lower authoritarian attitudes compared to a reference category, with males and individuals with higher education levels also showing lower authoritarian attitudes. However, the models have limited explanatory power, as indicated by low adjusted R-squared values, suggesting that additional unaccounted-for factors likely influence authoritarian attitudes.