# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "MASS", "aod", "interactions", "kableExtra", "flextable", "scales") # 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
## 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
## 
## 
## 
## Attaching package: 'MASS'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 
## 
## Attaching package: 'aod'
## 
## 
## The following object is masked from 'package:survival':
## 
##     rats
## 
## 
## 
## Attaching package: 'kableExtra'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
## 
## 
## 
## Attaching package: 'flextable'
## 
## 
## The following objects are masked from 'package:kableExtra':
## 
##     as_image, footnote
## 
## 
## The following object is masked from 'package:purrr':
## 
##     compose
## 
## 
## 
## Attaching package: 'scales'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
## [[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] "MASS"         "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] "aod"          "MASS"         "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"        
## 
## [[9]]
##  [1] "interactions" "aod"          "MASS"         "survey"       "survival"    
##  [6] "Matrix"       "grid"         "effects"      "carData"      "modelsummary"
## [11] "fst"          "infer"        "lubridate"    "forcats"      "stringr"     
## [16] "dplyr"        "purrr"        "readr"        "tidyr"        "tibble"      
## [21] "ggplot2"      "tidyverse"    "stats"        "graphics"     "grDevices"   
## [26] "utils"        "datasets"     "methods"      "base"        
## 
## [[10]]
##  [1] "kableExtra"   "interactions" "aod"          "MASS"         "survey"      
##  [6] "survival"     "Matrix"       "grid"         "effects"      "carData"     
## [11] "modelsummary" "fst"          "infer"        "lubridate"    "forcats"     
## [16] "stringr"      "dplyr"        "purrr"        "readr"        "tidyr"       
## [21] "tibble"       "ggplot2"      "tidyverse"    "stats"        "graphics"    
## [26] "grDevices"    "utils"        "datasets"     "methods"      "base"        
## 
## [[11]]
##  [1] "flextable"    "kableExtra"   "interactions" "aod"          "MASS"        
##  [6] "survey"       "survival"     "Matrix"       "grid"         "effects"     
## [11] "carData"      "modelsummary" "fst"          "infer"        "lubridate"   
## [16] "forcats"      "stringr"      "dplyr"        "purrr"        "readr"       
## [21] "tidyr"        "tibble"       "ggplot2"      "tidyverse"    "stats"       
## [26] "graphics"     "grDevices"    "utils"        "datasets"     "methods"     
## [31] "base"        
## 
## [[12]]
##  [1] "scales"       "flextable"    "kableExtra"   "interactions" "aod"         
##  [6] "MASS"         "survey"       "survival"     "Matrix"       "grid"        
## [11] "effects"      "carData"      "modelsummary" "fst"          "infer"       
## [16] "lubridate"    "forcats"      "stringr"      "dplyr"        "purrr"       
## [21] "readr"        "tidyr"        "tibble"       "ggplot2"      "tidyverse"   
## [26] "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
## [31] "methods"      "base"

Task 1

Use data for Germany and the variable coding for model3 and model4 from the tutorial (i.e., both for the cleaning & recode, as well as the linear model formulas). Do a modelsummary table displaying both model outputs. Interpret the coefficients for the MLR without the interaction (i.e., the first displayed model), and interpret model fit metrics for both models.

# Load Germany data
germany_data <- read_fst("germany_data.fst")
# Renaming dataset to df
df <- germany_data
# Survey weights
df$weight <- df$dweight * df$pweight
survey_design <- svydesign(ids = ~1, data = df, weights = ~weight)
# Multiple linear regression
# Set up our populist scale
df$trstplt <- ifelse(df$trstplt > 10, NA, df$trstplt)
df$trstprl <- ifelse(df$trstprl > 10, NA, df$trstprl)
df$trstprt <- ifelse(df$trstprt > 10, NA, df$trstprt)
df$trust <- scales::rescale(df$trstplt + df$trstprl + df$trstprt, na.rm = TRUE, to = c(0, 100))
df$populist <- scales::rescale(df$trust, na.rm = TRUE, to = c(100, 0))
# Recoding variables for model5 and model6
df <- df %>%
  mutate(
    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(yrbrn)
    ),
    gen = factor(gen,
                 levels = c("1", "2", "3", "4"),
                 labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
  )
# Running models
model3 <- lm(populist ~ gen, data = df, weights = weight)
model4 <- lm(populist ~ gen + cohort, data = df, weights = weight)

# Modelsummary table
modelsummary(
  list(model3, model4),
  fmt = 1,
  estimate = c("{estimate} ({std.error}){stars}",
               "{estimate} ({std.error}){stars}"),
  statistic = NULL
)
 (1)   (2)
(Intercept) 59.9 (0.3)*** 304.2 (49.0)***
genBaby Boomers 1.6 (0.3)*** 3.4 (0.6)***
genGen X 0.3 (0.4) 3.9 (0.9)***
genMillennials −3.0 (0.4)*** 2.7 (1.3)*
cohort −0.1 (0.0)***
Num.Obs. 28552 27899
R2 0.007 0.008
R2 Adj. 0.007 0.008
AIC 256975.5 251065.9
BIC 257016.8 251115.3
Log.Lik. −128482.771 −125526.939
RMSE 20.73 20.72

We find that all of the coefficients for the MLR without interaction (model3_germany) are significant at the 0.001 level. After correcting for other variables, there is a specific shift in the populist attitudes scale for every unit rise in the generation category (in relation to the reference group). Because interaction terms are involved, the interpretation of the MLR with interaction (model4_germany) becomes more intricate. The link between the predictors and the outcome varies based on the values of other predictors, as indicated by the coefficients for the interaction terms.

Task 2

Now generate the model4 interaction plot that we did in the tutorial, but again using the German data instead of the French. Interpret.

interaction_plot <- effect("gen*cohort", model4, na.rm = TRUE)
## NOTE: gen:cohort does not appear in the model
plot(interaction_plot,
     main = "Interaction Effect: Generational Cohort and Year of Birth on Populist Attitudes",
     xlab = "Generational Cohort",
     ylab = "Populist Attitudes Scale")

interaction_plot
## 
##  gen*cohort effect
##               cohort
## gen                1930     1950     1960     1980     2000
##   Interwar     61.38713 58.87057 57.61228 55.09572 52.57916
##   Baby Boomers 64.75681 62.24025 60.98197 58.46541 55.94885
##   Gen X        65.31135 62.79479 61.53651 59.01995 56.50339
##   Millennials  64.10352 61.58696 60.32868 57.81212 55.29556

The relationship between capital, generation, and the size of populist attitudes is visualized using an interaction plot. An interaction effect is shown if the lines that represent the various capital levels are not parallel to one another. Analyzing the differences in the line slopes and determining if the relationship between generation and result is contingent on capital level are important aspects of interpretation.

Task 3

Use data for the Netherlands and the variable coding for model5 and model6 from the tutorial (i.e., both for the cleaning & recode, as well as the linear model formulas). Do a modelsummary table displaying both model outputs. Interpret the coefficients for the MLR without the interaction (i.e., the first displayed model), and interpret model fit metrics for both models.

netherlands_data <- read_fst("netherlands_data.fst")
# Renaming dataset to df for Netherlands data
df_netherlands <- netherlands_data
# Survey weights
df$weight <- df$dweight * df$pweight
survey_design <- svydesign(ids = ~1, data = df, weights = ~weight)
# Recoding variables for model5 and model6
df <- df %>%
  mutate(
    polID = case_when(
      lrscale %in% 0:3 ~ "Left",
      lrscale %in% 7:10 ~ "Right",
      lrscale %in% 4:6 ~ "Moderate",
      lrscale %in% c(77, 88, 99) ~ NA_character_
    ),
    religious = case_when(
      rlgdgr %in% c(77, 88, 99) ~ NA_real_,
      TRUE ~ rlgdgr
    )
  )
# Creating the 'auth' variable for authoritarian attitudes scale
df <- df %>%
  mutate(
    behave = ipbhprp,
    secure = impsafe,
    safety = ipstrgv,
    tradition = imptrad,
    rules = ipfrule,
    across(c("behave", "secure", "safety", "tradition", "rules"),
           ~ na_if(.x, 7) %>% na_if(8) %>% na_if(9)),
    across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ),
    auth = scales::rescale(behave + secure + safety + tradition + rules, to=c(0,100), na.rm=TRUE)
  )
# Running models
model5 <- lm(auth ~ polID + religious, data = df, weights = weight)
model6 <- lm(auth ~ polID + religious + polID*religious, data = df, weights = weight)

# Modelsummary table
modelsummary(
  list(model5, model6),
  fmt = 1,
  estimate = c("{estimate} ({std.error}){stars}",
               "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_omit = "Intercept"
)
 (1)   (2)
polIDModerate 4.6 (0.3)*** 3.4 (0.4)***
polIDRight 7.7 (0.4)*** 4.2 (0.7)***
religious 0.9 (0.0)*** 0.6 (0.1)***
polIDModerate × religious 0.3 (0.1)***
polIDRight × religious 0.8 (0.1)***
Num.Obs. 23373 23373
R2 0.048 0.050
R2 Adj. 0.048 0.050
AIC 200594.1 200558.6
BIC 200634.4 200615.0
Log.Lik. −100292.049 −100272.298
RMSE 17.27 17.24

A comparison is made between models 5 and 6 in the Netherlands data. While model 6 has an interaction term to evaluate how the link between political ideology and authoritarian views vary depending on religiousness, model 5 looks at the individual impacts of political ideology and religiosity on authoritarian attitudes. Both political ideology and religiousness have significant effects, as indicated by the coefficients in model 5, where people who identify as politically right and those who are more religiously inclined score higher on the authoritarian attitudes scale. The sixth model presents an interaction term that implies that there are variations in the correlation between political ideology and authoritarian views based on religion affiliation. This emphasizes how crucial it is to take the interaction effect into account when examining how these factors relate to one another.

Task 4

Produce the model7 interaction plot from the tutorial, but again using the Netherlands data. Interpret.

# Recoding variables for model7
df <- df %>%
  mutate(
    behave = ipbhprp,
    secure = impsafe,
    safety = ipstrgv,
    tradition = imptrad,
    rules = ipfrule,
    across(c("behave", "secure", "safety", "tradition", "rules"),
           ~ na_if(.x, 7) %>% na_if(8) %>% na_if(9)),
    across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ),
    auth = scales::rescale(behave + secure + safety + tradition + rules, to=c(0,100), na.rm=TRUE)
  )
# Filtering out rows with NA values in the 'auth' variable
df <- df %>% filter(!is.na(auth))
# Running model7
model7 <- lm(auth ~ polID + cohort + polID*cohort, data = df, weights = weight)
interaction_plot <- effect("polID*cohort", model7, na.rm = TRUE)

plot(interaction_plot,
     main = "Interaction Effect: Political Ideology and Generational Cohort on Authoritarian Attitudes",
     xlab = "Political Ideology",
     ylab = "Authoritarian Attitudes Scale")

interaction_plot
## 
##  polID*cohort effect
##           cohort
## polID          1930     1950     1960     1980     2000
##   Left     66.75900 61.83047 59.36620 54.43767 49.50915
##   Moderate 70.45384 66.29155 64.21041 60.04812 55.88583
##   Right    74.46138 69.33910 66.77796 61.65569 56.53341

The interaction impact of generational cohort and political ideology on authoritarian views is shown in the model7 interaction plot for the Netherlands data. The plot shows how different generational cohorts have varied relationships between political ideology and authoritarian tendencies. Regression line slopes for the left, moderate, and right political ideology groups within each generational cohort are displayed, enabling a visual assessment of the ways in which these variables interact to shape authoritarian attitudes.

Task 5

Produce the model8 interaction plot from the tutorial, but again using the Netherlands data. Interpret.

# Recoding variables for model8
df <- df %>%
  mutate(
    behave = ipbhprp,
    secure = impsafe,
    safety = ipstrgv,
    tradition = imptrad,
    rules = ipfrule,
    across(c("behave", "secure", "safety", "tradition", "rules"),
           ~ na_if(.x, 7) %>% na_if(8) %>% na_if(9)),
    across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ),
    auth = scales::rescale(behave + secure + safety + tradition + rules, to=c(0,100), na.rm=TRUE),
    ID = case_when(
      lrscale >= 0 & lrscale <= 4 ~ "Left",
      lrscale >= 6 & lrscale <= 10 ~ "Right",
      lrscale > 10 ~ NA_character_,  # Set values above 10 as NA
      TRUE ~ NA_character_  # Ensure value 5 and any other unexpected values are set as NA
    )
  )
# Filtering out rows with NA values in the 'auth' variable
df <- df %>% filter(!is.na(auth))
# Running model8
model8 <- lm(auth ~ religious + ID + religious*ID, data = df, weights = weight)
# Interaction plot for model8
interaction_plot <- effect("religious*ID", model8, na.rm = TRUE)

plot(interaction_plot,
     main = "Interaction Effect: Religiousness and Ideological Self-placement on Authoritarian Attitudes",
     xlab = "Religiousness",
     ylab = "Authoritarian Attitudes Scale")

interaction_plot
## 
##  religious*ID effect
##          ID
## religious     Left    Right
##        0  56.71321 60.03449
##        2  58.04484 62.66420
##        5  60.04229 66.60877
##        8  62.03973 70.55333
##        10 63.37136 73.18304

Lastly, the interaction impact of religiousness and ideological self-placement on authoritarian attitudes is seen in the model 8 interaction plot for the Netherlands data. The storyline illustrates how those who identify as politically left and right have different relationships between religiousness and authoritarian tendencies. The graphic representation of the regression line slopes for every religion level within each group of ideological self-placement enables an understanding of the ways in which these variables interact to shape authoritarian attitudes.