# 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.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ 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
## Version 2.0.0 of `modelsummary`, to be released soon, will introduce a
##   breaking change: The default table-drawing package will be `tinytable`
##   instead of `kableExtra`. All currently supported table-drawing packages
##   will continue to be supported for the foreseeable future, including
##   `kableExtra`, `gt`, `huxtable`, `flextable, and `DT`.
##   
##   You can always call the `config_modelsummary()` function to change the
##   default table-drawing package in persistent fashion. To try `tinytable`
##   now:
##   
##   config_modelsummary(factory_default = 'tinytable')
##   
##   To set the default back to `kableExtra`:
##   
##   config_modelsummary(factory_default = 'kableExtra')
## 
## 载入需要的程辑包：carData
## 
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## 
## 载入需要的程辑包：grid
## 
## 载入需要的程辑包：Matrix
## 
## 
## 载入程辑包：'Matrix'
## 
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## 
## 载入需要的程辑包：survival
## 
## 
## 载入程辑包：'survey'
## 
## 
## The following object is masked from 'package:graphics':
## 
##     dotchart
## 
## 
## 
## 载入程辑包：'MASS'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 
## 
## 载入程辑包：'aod'
## 
## 
## The following object is masked from 'package:survival':
## 
##     rats
## 
## 
## 
## 载入程辑包：'kableExtra'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
## 
## 
## 
## 载入程辑包：'flextable'
## 
## 
## The following objects are masked from 'package:kableExtra':
## 
##     as_image, footnote
## 
## 
## The following object is masked from 'package:purrr':
## 
##     compose
## 
## 
## 
## 载入程辑包：'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")
first <- germany_data
first <- first %>%
  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 ))

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

first <- first %>% filter(!is.na(auth))
first <- first %>%
  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 <- lm(auth ~ polID + religious, data = first, weights = weight)
model4 <- lm(auth ~ polID + religious + polID*religious, data = first, weights = weight)
modelsummary(
  list(model3, model4),
  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

pollDRight(9.3) is represent The level of religiosity is 9.3 points higher for those on the political right compared to those on the “left” ideology. pollDModerate(5.3) is In groups with similar levels of religiosity, politically moderate people scored on average 5.3 points higher compared to people with leftist ideologies. religious(0.9) is represent for each unit increase in religiosity, the mean expected score on the Authoritarian Values Scale increases by 0.9 units.

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

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

pollDLeft images indicate a lower level of religiosity and pollDModerated images indicate a more subdued level of religiosity.pollDRight images indicate high religiosity.

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")
third <- nld_data
third$weight <- third$dweight * third$pweight
survey_design2 <- svydesign(ids = ~1, data = third, weights = ~weight)
third$trstplt <- ifelse(third$trstplt > 10, NA, third$trstplt)
third$trstprl <- ifelse(third$trstprl > 10, NA, third$trstprl)
third$trstprt <- ifelse(third$trstprt > 10, NA, third$trstprt)

third$trust <- scales::rescale(third$trstplt + third$trstprl + third$trstprt, na.rm = TRUE, to = c(0, 100))

third$populist <- scales::rescale(third$trust, na.rm = TRUE, to=c(100,0))
third <- third %>%
  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"))
  )
table(third$gen)
## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         3967         6359         4684         2803
third <- third %>%
  mutate(religion = case_when(
    rlgblg == 2 ~ "No",
    rlgblg == 1 ~ "Yes",
    rlgblg %in% c(7, 8, 9) ~ NA_character_,
    TRUE ~ as.character(rlgblg)
  ))

# check
table(third$religion)
## 
##    No   Yes 
## 11289  7008
third <- third %>%
  mutate(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
  ))
table(third$ID)
## 
##  Left Right 
##  5605  7247
third <- third %>%
  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 ))

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


third <- third %>% filter(!is.na(auth))
model5 <- lm(auth ~ religion + ID + gen, data = third, weights = weight)
model6 <- lm(auth ~ religion + ID + gen + religion*gen, data = third, weights = weight)
modelsummary(
  list(model5, model6),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_omit = "Intercept")
 (1)   (2)
religionYes 7.2 (0.3)*** 5.7 (0.7)***
IDRight 4.2 (0.3)*** 4.3 (0.3)***
genBaby Boomers -5.2 (0.4)*** -5.9 (0.6)***
genGen X -7.7 (0.4)*** -8.8 (0.6)***
genMillennials -8.4 (0.5)*** -9.7 (0.6)***
religionYes × genBaby Boomers 1.1 (0.8)
religionYes × genGen X 2.3 (0.9)**
religionYes × genMillennials 2.9 (1.0)**
Num.Obs. 12116 12116
R2 0.124 0.125
R2 Adj. 0.124 0.125
AIC 101219.3 101213.5
BIC 101271.2 101287.5
Log.Lik. -50602.672 -50596.757
RMSE 15.27 15.26

ReligionYes 7.2 means joining a religion for those with no religion adds 7.2 units.

IDRight 4.2 When all other factors are held constant, a one-unit rise in IDRight results in a 4.2-unit rise in the variable.

genBaby Boomers 5.2 The baby boomer generation would be reduced by 5.2 variables.

genGen X(7.7) genGen X generation would be reduced by 7.7 variables.

genMillennials(8.4) The variable for those categorized as Millennials decreases by 8.4.

R-Squared model 0.124 is represent Models with no interaction accounted for 12.4% of the change in the dependent variable.

Task 4

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

netherlands_data <- read_fst("netherlands_data.fst")
dn <- netherlands_data
dn$weight <- dn$dweight * dn$pweight
survey_design <- svydesign(ids = ~1, data = dn, weights = ~weight)
dn <- dn %>%
  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 ))


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


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

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


dn <- dn %>% filter(!is.na(auth))
dn <- dn %>%
  mutate(
    cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
    # Recoding generational cohorts based on the year of birth (yrbrn).
    # The year of birth is categorized into different generational cohorts.
    # Interwar (1900-1945), Baby Boomers (1946-1964), Gen X (1965-1979), Millennials (1980-1996).
    # The 'TRUE' line is a catch-all that keeps the original year of birth for those not in these ranges.
    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)  
    ),
    # After recoding, the gen variable is converted into a factor with labels for clearer interpretation.
    # Factors are used in R to handle categorical variables.
    gen = factor(gen,
                 levels = c("1", "2", "3", "4"),
                 labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
  )
table(dn$gen)
## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         3771         6147         4548         2730
model7 <- lm(auth ~ + cohort + polID + cohort*polID, data = dn, weights = weight)
interact_plot(model7, pred = cohort, modx = polID, jnplot = TRUE)

This let image shows the connection between rightist political formations, leftist political formations and moderate political formations. It can be seen in the graph that in the year 2000, the lines for the moderate and rightist political formations intersected, indicating that the data for the rightists and moderates were the same in that year. Since 1950, all three political formations have shown a downward trend. The leftist political pattern is the lowest, and the gap between the rightist and moderate political patterns gets smaller with each passing year.

Task 5

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

dn1 <- netherlands_data
dn1$weight <- dn1$dweight * dn1$pweight
survey_design <- svydesign(ids = ~1, data = dn1, weights = ~weight)
dn1 <- dn1 %>%
  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 ))


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


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

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


dn1 <- dn1 %>% filter(!is.na(auth))
dn1 <- dn1 %>%
  mutate(
    cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
    # Recoding generational cohorts based on the year of birth (yrbrn).
    # The year of birth is categorized into different generational cohorts.
    # Interwar (1900-1945), Baby Boomers (1946-1964), Gen X (1965-1979), Millennials (1980-1996).
    # The 'TRUE' line is a catch-all that keeps the original year of birth for those not in these ranges.
    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)  
    ),
    # After recoding, the gen variable is converted into a factor with labels for clearer interpretation.
    # Factors are used in R to handle categorical variables.
    gen = factor(gen,
                 levels = c("1", "2", "3", "4"),
                 labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
  )
table(dn1$gen)
## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         3771         6147         4548         2730
dn1 <- dn1 %>%
  mutate(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
  ))
table(dn1$ID)
## 
##  Left Right 
##  5449  7010
model8 <- lm(auth ~ religious + ID + religious*ID, data = dn1, weights = weight)
interact_plot(model8, pred = religious, modx = ID, jnplot = TRUE)

This image shows that the left and the right political formations are two unrelated. Both show an upward trend. Rightist political formations are always higher than leftist political formations.