rm(list=ls()); gc() 
##          used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 527659 28.2    1172837 62.7         NA   669520 35.8
## Vcells 973558  7.5    8388608 64.0      16384  1851694 14.2
# 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.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"
germany_data <- read_fst("germany_data.fst")
dg <- germany_data

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.

dg$weight <- dg$dweight * dg$pweight
survey_design <- svydesign(ids = ~1, data = dg, weights = ~weight)
# Setting values greater than 10 to NA
dg$trstplt <- ifelse(dg$trstplt > 10, NA, dg$trstplt)
dg$trstprl <- ifelse(dg$trstprl > 10, NA, dg$trstprl)
dg$trstprt <- ifelse(dg$trstprt > 10, NA, dg$trstprt)

# Creating and rescaling the trust variable
dg$trust <- scales::rescale(dg$trstplt + dg$trstprl + dg$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
dg$populist <- scales::rescale(dg$trust, na.rm = TRUE, to=c(100,0))
dg <- dg %>%
  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(dg$gen)
## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         6373        11474         7648         6704
# The dataframe is further updated to recode trust and political interest # note: you could have done it one, integrating above, but I broke it down
dg <- dg %>%
  mutate(
    # Recoding the 'ppltrst' variable, which seems to measure trust on a scale from 0 to 10.
    # Special codes like 77, 88, 99 are probably used for missing data and are set to NA.
    # Trust levels are categorized as 'Low', 'Mid', and 'High'.
    Capital = case_when(
      ppltrst %in% c(77, 88, 99) ~ NA_character_,
      ppltrst >= 0 & ppltrst <= 3 ~ "Low",
      ppltrst >= 4 & ppltrst <= 6 ~ "Mid",
      ppltrst >= 7 & ppltrst <= 10 ~ "High",
      TRUE ~ as.character(ppltrst)
    ),
    # Recoding 'polintr' variable, which seems to measure political interest.
    # Values of 1 and 2 indicate 'Interested', 3 and 4 'Not Interested'.
    # Special codes like 7, 8, 9 (likely missing data) are set to NA.
    Politically = case_when(
      polintr %in% c(1, 2) ~ "Interested",
      polintr %in% c(3, 4) ~ "Not Interested",
      polintr %in% c(7, 8, 9) ~ NA_character_
    )
  )
model1 <- lm(populist ~ gen, data = dg, weights = weight)
model2 <- lm(populist ~ Capital + Politically, data = dg, weights = weight)
modelsummary(
  list(model1, model2), # you can add as many models as you have generated by putting them in this list. But you MUST refer to them to their stored name. So if you named your model say Regression1, you would have to put that instead of model1. By default, you can use the naming convention we rely upon of model1, model2, ...
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}", # what we are specifying here is what we want to show in our model summary table
                  "{estimate} ({std.error}){stars}"),
  statistic = NULL)
 (1)   (2)
(Intercept) 59.9 (0.3)*** 49.8 (0.2)***
genBaby Boomers 1.6 (0.3)***
genGen X 0.3 (0.4)
genMillennials −3.0 (0.4)***
CapitalLow 17.3 (0.3)***
CapitalMid 8.3 (0.3)***
PoliticallyNot Interested 4.6 (0.2)***
Num.Obs. 28552 30565
R2 0.007 0.121
R2 Adj. 0.007 0.121
AIC 256975.5 271785.2
BIC 257016.8 271826.9
Log.Lik. −128482.771 −135887.607
RMSE 20.73 19.51 
dg <- dg %>%
  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
dg$auth <- scales::rescale(dg$behave + 
                      dg$secure + 
                      dg$safety + 
                      dg$tradition + 
                      dg$rules, to=c(0,100), na.rm=TRUE)


dg <- dg %>% filter(!is.na(auth))
dg <- dg %>%
  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 = dg, weights = weight)
model4 <- lm(auth ~ polID + religious + polID*religious, data = dg, weights = weight)
modelsummary(
  list(model3, model4),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL)
 (1)   (2)
(Intercept) 55.3 (0.3)*** 56.4 (0.3)***
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

Model 1: Intercept (55.3): This coefficient suggests that the expected outcome value, assuming all predictors are at their reference level or zero (for continuous variables like ‘religious’), is 55.3. The stars indicate that this result is highly significant, with a p-value < 0.001.

polIDModerate (4.6): Individuals identifying as politically moderate, holding all else constant, score on average 4.6 units higher on the outcome variable compared to the baseline group (assumed to be ‘left’ political identification), with high significance.

polIDRight (7.7):Similarly, those identifying as politically right are expected to score, on average, 7.7 units higher than the baseline political identification group, again with high significance.

religious (0.9):For each unit increase in religiousness, there is an associated expected increase of 0.9 units in the outcome variable, holding other variables constant. This is highly significant.

Model 2: Intercept (56.4):The expected value of the outcome variable at the reference levels of predictors is slightly higher in this model than in the first, indicating a base shift. polIDModerate (3.4):The effect of being politically moderate has decreased compared to model (1), suggesting a 3.4 units increase over the baseline, with significant certainty.

polIDRight (4.2):The effect of identifying as politically right also sees a reduction to 4.2 units higher than the baseline, which is significant but less pronounced than in model (1).

religious (0.6):The effect of religiousness decreases to 0.6 units increase for every unit increase in religiousness, showing high significance but a smaller effect size than in model (1).

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

The interaction plot illustrates the relationship between religiousness and authoritarian attitudes across three political identifications: Left, Moderate, and Right. The slopes represent how religiousness predicts authoritarian attitudes for each political group. For the Left, there’s a positive but relatively gentle slope, indicating a slight increase in authoritarian attitudes with increased religiousness. The Moderate slope is steeper, suggesting a more substantial increase in authoritarian attitudes as religiousness rises. The Right shows the steepest slope of all, implying that increases in religiousness are associated with a more pronounced increase in authoritarian attitudes compared to the other groups. Notably, the lines are not parallel, indicating interaction effects: the influence of religiousness on authoritarian attitudes varies by political identification. This non-parallel trend signifies that the relationship between religiousness and authoritarian attitudes is not consistent across political groups; it intensifies from Left to Right.

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")
dn <- netherlands_data
dn$weight <- dn$dweight * dn$pweight
survey_design <- svydesign(ids = ~1, data = dn, weights = ~weight)
# Setting values greater than 10 to NA
dn$trstplt <- ifelse(dn$trstplt > 10, NA, dn$trstplt)
dn$trstprl <- ifelse(dn$trstprl > 10, NA, dn$trstprl)
dn$trstprt <- ifelse(dn$trstprt > 10, NA, dn$trstprt)

# Creating and rescaling the trust variable
dn$trust <- scales::rescale(dn$trstplt + dn$trstprl + dn$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
dn$populist <- scales::rescale(dn$trust, na.rm = TRUE, to=c(100,0))
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 
##         3967         6359         4684         2803
# The dataframe is further updated to recode trust and political interest # note: you could have done it one, integrating above, but I broke it down
dn <- dn %>%
  mutate(
    # Recoding the 'ppltrst' variable, which seems to measure trust on a scale from 0 to 10.
    # Special codes like 77, 88, 99 are probably used for missing data and are set to NA.
    # Trust levels are categorized as 'Low', 'Mid', and 'High'.
    Capital = case_when(
      ppltrst %in% c(77, 88, 99) ~ NA_character_,
      ppltrst >= 0 & ppltrst <= 3 ~ "Low",
      ppltrst >= 4 & ppltrst <= 6 ~ "Mid",
      ppltrst >= 7 & ppltrst <= 10 ~ "High",
      TRUE ~ as.character(ppltrst)
    ),
    # Recoding 'polintr' variable, which seems to measure political interest.
    # Values of 1 and 2 indicate 'Interested', 3 and 4 'Not Interested'.
    # Special codes like 7, 8, 9 (likely missing data) are set to NA.
    Politically = case_when(
      polintr %in% c(1, 2) ~ "Interested",
      polintr %in% c(3, 4) ~ "Not Interested",
      polintr %in% c(7, 8, 9) ~ NA_character_
    )
  )
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))) %>%
  # Apply the reverse coding
  mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))

# Now you can calculate 'schwartzauth' after the NA recoding
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(religion = case_when(
    rlgblg == 2 ~ "No",
    rlgblg == 1 ~ "Yes",
    rlgblg %in% c(7, 8, 9) ~ NA_character_,
    TRUE ~ as.character(rlgblg)
  ))

# check
table(dn$religion)
## 
##    No   Yes 
## 10913  6750
dn <- dn %>%
  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(dn$ID)
## 
##  Left Right 
##  5449  7010
model5 <- lm(auth ~ religion + ID + gen, data = dn, weights = weight)
model6 <- lm(auth ~ religion + ID + gen + religion*gen, data = dn, 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

Model 1: religionYes (7.2):Individuals who identify as religious have an expected 7.2 units higher score on the outcome variable, likely authoritarian attitudes, compared to those who are not religious, with a high level of statistical significance.

IDRight (4.2):Being politically right is associated with an average increase of 4.2 units on the outcome variable compared to the baseline political identification, again with strong statistical significance.

genBaby Boomers (−5.2): Baby Boomers score, on average, 5.2 units lower on the outcome variable compared to the reference generational group (presumably the Silent Generation), with a highly significant p-value.

genGen X (−7.7):Generation X scores 7.7 units lower on the outcome variable than the reference generational group, with high significance.

genMillennials (−8.4): Millennials show an 8.4 unit decrease in the outcome variable score compared to the reference generational group, which is also highly significant.

Model 2: religionYes × genBaby Boomers (1.1):he interaction between being religious and being a Baby Boomer is positive but not statistically significant.

religionYes × genGen X (2.3):Being religious and from Generation X is associated with an additional 2.3 unit increase in the outcome variable, significant at the 0.05 level.

religionYes × genMillennials (2.9):Being religious and a Millennial is associated with an additional 2.9 unit increase in the outcome variable, which is significant at the 0.05 level.

Task 4

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

model7 <- lm(auth ~ religion + ID + gen, data = dn, weights = weight)
interaction_plot <- effect("religion*gen", model7, na.rm=TRUE)
## NOTE: religion:gen does not appear in the model
plot(interaction_plot,
     main="Interaction effect",
     xlab="Generations",
     ylab="Authoritarian attitudes scale")

interaction_plot
## 
##  religion*gen effect
##         gen
## religion Interwar Baby Boomers    Gen X Millennials
##      No  65.54806     60.34304 57.89188    57.10745
##      Yes 72.78179     67.57677 65.12561    64.34118

The interaction plot displays a downward trend in authoritarian attitudes from the Interwar generation to Millennials, with non-religious individuals showing a more pronounced decrease than their religious counterparts. Although starting higher, the authoritarian attitudes of non-religious individuals steeply decline across generations, suggesting a strong generational effect on these attitudes. For religious individuals, the decline is less steep, indicating a more moderate generational effect. Interestingly, the gap in authoritarian attitudes between religious and non-religious individuals narrows over time, particularly among Millennials, suggesting that the influence of religion on authoritarian attitudes diminishes with younger generations. This convergence implies a potential shift in the way generational and religious factors interact to shape political attitudes.

Task 5

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

model8 <- lm(auth ~ religion + ID + gen + religion*gen, data = dn, weights = weight)
interaction_plot <- effect("religion*gen", model8, na.rm=TRUE)

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

interaction_plot
## 
##  religion*gen effect
##         gen
## religion Interwar Baby Boomers    Gen X Millennials
##      No  66.39342     60.50327 57.62650    56.73392
##      Yes 72.13454     67.31264 65.70738    65.34644

The interaction plot depicts how generational cohorts’ authoritarian attitudes vary by religious affiliation. Non-religious individuals show a pronounced downward trend from the Interwar generation to Millennials, suggesting a significant generational shift towards less authoritarian attitudes over time. In contrast, religious individuals start with lower authoritarian attitudes in the older Interwar cohort than their non-religious counterparts and exhibit a more gradual decline, with the youngest Millennial generation converging towards similar attitudes regardless of religious affiliation. The trends imply that while generational shifts in authoritarian attitudes are evident in both religious and non-religious groups, the magnitude of this shift is more pronounced among the non-religious, but both groups appear to converge in the Millennial cohort. The error bars suggest some variability in the estimates but do not overlap between generations, reinforcing the generational effect.