packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "MASS", "aod", "interactions", "kableExtra", "flextable", "scales") 

new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_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
## 
## 
## 
## 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.

germany_data <- read_fst("germany_data.fst")

task1 <- germany_data
task1 <- task1 %>%
  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 ))

task1$auth <- scales::rescale(task1$behave + 
                      task1$secure + 
                      task1$safety + 
                      task1$tradition + 
                      task1$rules, to=c(0,100), na.rm=TRUE)

task1 <- task1 %>% filter(!is.na(auth))
task1 <- task1 %>%
  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
    ))
model3.1 <- lm(auth ~ polID + religious, data = task1, weights = weight)
model4.1 <- lm(auth ~ polID + religious + polID*religious, data = task1, weights = weight)
modelsummary(
  list(model3.1, model4.1),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL)
 (1)   (2)
(Intercept) 57.5 (0.7)*** 59.6 (0.9)***
polIDModerate 5.3 (0.7)*** 2.6 (1.1)*
polIDRight 9.3 (1.1)*** 3.4 (1.9)+
religious 0.3 (0.1)** −0.3 (0.2)+
polIDModerate × religious 0.8 (0.2)**
polIDRight × religious 1.5 (0.4)***
Num.Obs. 2855 2855
R2 0.035 0.041
R2 Adj. 0.034 0.039
AIC 24680.0 24665.5
BIC 24709.7 24707.2
Log.Lik. −12334.978 −12325.764
RMSE 17.13 17.13

A coefficient of 5.3 for polIDModerate in model 1 means that in comparison to individuals who are left-leaning, those who score in the middle on the left-right political scale will score approximately 5.3 units higher on the authoritarian scale than left-leaning individuals.

A coefficient of 9.3 for polIDRight represents indicates that right-leaning individuals score an average of 9.3 units higher on the authoritarian scale than left-leaning individuals.

polIDModerate x religious has a coefficient of 0.8, indicating that on average, individuals who score moderately on the left-right scale and are more religious score 0.8 units higher than left-leaning religious individuals on the authoritarian scale. However, the asterisks signify that this coefficient is statistically significant, but does not have the lowest p-value.

polIDRight x religious indicates that religious individuals who are right-leaning score an average of 1.5 units higher than left-leaning religious individuals on the authoritarian scale.

polIDModerate, polIDRight, and polIDRight x religious have the most statistically significant coefficients, meaning they have a very small p-value and the null hypothesis can be rejected. Along with this, an r-squared adjusted value of 0.034 and 0.039 or 3.4% and 3.9% means that the intercepts explain only 3.4% and 3.9% of the variance in the models.

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_plot2 <- effect("polID*religious", model4.1, na.rm=TRUE)

plot(interaction_plot2,
     main="Interaction effect",
     xlab="Religiousness",
     ylab="Authoritarian attitudes scale")

interaction_plot2
## 
##  polID*religious effect
##           religious
## polID             0        2        5        8       10
##   Left     59.61010 58.93398 57.91980 56.90562 56.22950
##   Moderate 62.21299 63.14317 64.53843 65.93370 66.86388
##   Right    63.00456 65.26702 68.66072 72.05442 74.31689

The slope of these lines are not parallel to one another and indicate that there is an interaction effect between an individual’s authoritarian score religiousness and their political standing.

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.

nld_data <- read_fst("netherlands_data.fst")

task3 <- nld_data
task3 <- task3 %>%
  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 ))

task3$auth <- scales::rescale(task3$behave + 
                      task3$secure + 
                      task3$safety + 
                      task3$tradition + 
                      task3$rules, to=c(0,100), na.rm=TRUE)


task3 <- task3 %>% filter(!is.na(auth))
task3 <- task3 %>%
  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_      
    ))

task3 <- task3 %>%
   mutate(religion = case_when(
    rlgblg == 2 ~ "No",
    rlgblg == 1 ~ "Yes",
    rlgblg %in% c(7, 8, 9) ~ NA_character_,
    TRUE ~ as.character(rlgblg)
  ))

task3 <- task3 %>%
  mutate(ID = case_when(
    lrscale >= 0 & lrscale <= 4 ~ "Left",
    lrscale >= 6 & lrscale <= 10 ~ "Right",
    lrscale > 10 ~ NA_character_, 
    TRUE ~ NA_character_  
  ))
task3 <- task3 %>%
  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"))
  )
model5.3 <- lm(auth ~ religion + ID + gen, data = task3, weights = weight)
model6.3 <- lm(auth ~ religion + ID + gen + religion*gen, data = task3, weights = weight)
modelsummary(
  list(model5.3, model6.3),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL)
 (1)   (2)
(Intercept) 63.9 (1.3)*** 64.2 (1.6)***
religionYes 8.2 (0.9)*** 7.5 (2.1)***
IDRight 3.4 (0.9)*** 3.4 (0.9)***
genBaby Boomers −5.8 (1.3)*** −6.5 (1.7)***
genGen X −7.3 (1.4)*** −7.6 (1.8)***
genMillennials −8.6 (1.5)*** −8.6 (1.9)***
religionYes × genBaby Boomers 1.6 (2.6)
religionYes × genGen X 0.5 (2.8)
religionYes × genMillennials −0.6 (3.2)
Num.Obs. 1147 1147
R2 0.126 0.127
R2 Adj. 0.123 0.121
AIC 9479.3 9484.5
BIC 9514.7 9535.0
Log.Lik. −4732.669 −4732.251
RMSE 14.94 14.94

On average, individuals who identify as belonging to a particular religion score an average of 8.2 units higher on the authoritarian scale than individuals who do not belong to a religion.

A coefficient of 3.4 for right-leaning individuals represents an average authoritarian score 3.4 units higher than left-leaning individuals.

A coefficient of -5.8 for genBaby Boomers indicates that in comparison to the interwar generation, they score an average of 5.8 units lower on the authoritarian scale.

A coefficient of -7.3 for genGenX means that Generation X on average 7.3 units lower on the authoritarian scale than the interwar generation.

On average, millenials scored 8.6 units lower on the authoritarian than the interwar generation.

The lack of asterisks on the coefficients associated with generation and a person being religious in comparison to non-religious interwar individuals indicate that these values are not statistically significant while all other coefficients are statistically significant.

Task 4

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

model7.4 <- lm(auth ~ + cohort + polID + cohort*polID, data = task3, weights = weight)
interact_plot(model7.4, pred = cohort, modx = polID, jnplot = TRUE)

This model indicates that there is an interaction effect on a persons authoritarian scale and their generation and political stance. The lines have varying slopes and intersect on the graph rather than having parallel slopes, indicating that there is an interaction effect.

Task 5

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

task5 <- task3 %>%
  mutate(
    religious = case_when(
      rlgdgr %in% c(77, 88, 99) ~ NA_real_,
      TRUE ~ rlgdgr))
model8.5 <- lm(auth ~ religious + ID + religious*ID, data = task5, weights = weight)
interact_plot(model8.5, pred = religious, modx = ID, jnplot = TRUE)

While the lines appear to be parallel with no visible intersections on this graph, the slopes of the line are not the same so they are parallel which indicates that there is an interaction effect.

End