Setting up Environment

# 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
## Warning: 程辑包'infer'是用R版本4.3.3 来建造的
## Warning: 程辑包'effects'是用R版本4.3.3 来建造的
## 载入需要的程辑包:carData
## Warning: 程辑包'carData'是用R版本4.3.3 来建造的
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## Warning: 程辑包'survey'是用R版本4.3.3 来建造的
## 载入需要的程辑包: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
## Warning: 程辑包'aod'是用R版本4.3.3 来建造的
## 
## 载入程辑包:'aod'
## 
## The following object is masked from 'package:survival':
## 
##     rats
## Warning: 程辑包'interactions'是用R版本4.3.3 来建造的
## 
## 载入程辑包:'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
## Warning: 程辑包'flextable'是用R版本4.3.3 来建造的
## 
## 载入程辑包:'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")

df <- germany_data
df$weight <- df$dweight * df$pweight
survey_design <- svydesign(ids = ~1, data = df, weights = ~weight)
# Setting values greater than 10 to NA
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)

# Creating and rescaling the trust variable
df$trust <- scales::rescale(df$trstplt + df$trstprl + df$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
df$populist <- scales::rescale(df$trust, na.rm = TRUE, to=c(100,0))
df <- df %>%
  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(df$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
df <- df %>%
  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_
    )
  )
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))
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
    )
  )
model3 <- lm(auth ~ polID + religious, data = df, weights = weight)
model4 <- lm(auth ~ polID + religious + polID*religious, data = df, 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

Interpretation:

polIDModerate (4.6): This coefficient suggests that, comparing individuals with the same religiousness (i.e., what some term “holding constant”), the average outcome for individuals identifying as politically moderate are 4.6 units higher compared to the reference category (in our case, identifying ideologically as on the ‘left’). The “baseline group” refers to the category that is omitted, and is thus interpreted as part of our intercept.

polIDRight (7.7): This coefficient indicates that, comparing individuals with the same religiousness, individuals identifying as politically right are expected to be, on average, 7.7 units higher than the reference category (i.e., identifying ideologically as on the ‘left’) on the authoritarian values scale.

religious (0.9): For the ‘religious’ predictor, a coefficient of 0.9 means that with every one-unit increase in religiousness, there is an associated expected average increase of 0.9 units on the authoritarian values scale, assuming that we are comparing people with the same ideological categorization (i.e., holding ideological self-placement constant).

All three coefficients have ***, indicating a p < 0.001. Meaning: There is less than a 0.1% probability that the observed result is due to chance, which suggests that we can be very confident in asserting the existence of a relationship between the predictor and the outcome variable, under frequentist assumptions.

For model3 (without interaction terms):

R-squared (0.048): This indicates that 4.8% of the variance in the authoritarian values scale can be explained by the predictors included in the model. In other words, this model captures a modest portion of the variation in authoritarian values.

Adjusted R-squared (0.048): This metric adjusts the R-squared to account for the number of predictors in the model relative to the number of observations. An adjusted R-squared of 0.048 suggests that after adjustment for the number of predictors, the explanatory power is still comparable to the R-squared, which is indicative of a model that is not overly complex for the amount of data.

AIC (200594.1): The Akaike Information Criterion helps to compare models with a trade-off between the goodness of fit and model complexity. A lower AIC value generally indicates a better model. An AIC of 200594.1 for this model provides a baseline for comparison with other models.

For model4 (with interaction terms):

R-squared (0.050): A slightly higher R-squared value than model 3 suggests that the addition of interaction terms offers a marginal increase in explaining the variance in authoritarian values.

Adjusted R-squared (0.050): Despite the slight increase in R-squared, the adjusted R-squared remains the same as that of model 3, indicating that the addition of interaction terms does not significantly improve the explanatory power of the model when considering the penalty for additional predictors.

AIC (200558.6): The AIC is slightly lower than that of model 3, suggesting that the added complexity of the interaction terms might be justifiable in terms of model fit. Typically, we prefer models with lower AIC values, so in this comparison, model 4 might be marginally more preferable.

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

interaction_plot
## 
##  polID*religious effect
##           religious
## polID             0        2        5        8       10
##   Left     56.40706 57.56775 59.30879 61.04983 62.21053
##   Moderate 59.85063 61.66176 64.37844 67.09513 68.90626
##   Right    60.60237 63.35032 67.47224 71.59416 74.34211

Interpretation:

There is an overall positive relationship between religiousness and authoritarian attitudes across all political ideology groups. This means that as religiousness increases, so does the score on the authoritarian attitudes scale, regardless of political ideology.

The slopes of the lines for the Moderate and Right groups appear steeper than for the Left group. This suggests that the increase in authoritarian attitudes with respect to religiousness is more pronounced for individuals who identify as Moderate or Right compared to those who identify as Left.

The difference in slopes between the Moderate and Right groups compared to the Left group implies that the effect of religiousness on authoritarian attitudes is moderated by political ideology. That is, political ideology influences how strongly religiousness predicts authoritarian attitudes.

Takeaway:

While religiousness is associated with higher levels of authoritarian attitudes across all political spectrums, this association is particularly stronger for those who identify as Moderate or Right in their political ideology. It’s important to note that interaction effects are about the differences in slopes and not necessarily the starting or ending points of the lines. Hence, what’s significant here is that the trend is not parallel, indicating an interaction effect where the relationship between religiousness and authoritarian attitudes is different across the categories of political ideology.

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

df <- netherlands_data
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))
df <- df %>%
  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(df$gen)
## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         3771         6147         4548         2730
df <- df %>%
  mutate(religion = case_when(
    rlgblg == 2 ~ "No",
    rlgblg == 1 ~ "Yes",
    rlgblg %in% c(7, 8, 9) ~ NA_character_,
    TRUE ~ as.character(rlgblg)
  ))

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

Interpretation:

ReligionYes (8.2)*: This coefficient suggests that, comparing individuals with the same political ideology (i.e., holding constant), the average authoritarian value for religious individuals is 8.2 units higher compared to non-religious individuals. The reference category here is the non-religious individuals.

IDRight (3.4)*: This coefficient indicates that, comparing individuals with the same level of religiousness, individuals identifying as politically right are expected to be, on average, 3.4 units higher on the authoritarian values scale compared to those identifying ideologically as on the ‘left’, which is the reference category.

genBabyBoomers (-5.8), genGenX (-7.3), and genMillennials (-8.6)*: These coefficients suggest that, holding religiousness and political ideology constant, Baby Boomers, Generation X, and Millennials have, on average, -5.8, -7.3, and -8.6 units lower on the authoritarian values scale, respectively, compared to the reference category, which is the Interwar generation.

All three coefficients have ***, indicating a p-value less than 0.001, meaning there is a very low probability (less than 0.1%) that these results could occur by chance, suggesting strong evidence against the null hypothesis.

For model5 (without interaction terms):

R-squared (0.126): This indicates that 12.6% of the variance in the authoritarian values scale can be explained by the predictors included in the model. In other words, this model captures a modest portion of the variation in authoritarian values.

Adjusted R-squared (0.123): This metric adjusts the R-squared to account for the number of predictors in the model relative to the number of observations. An adjusted R-squared of 0.123 suggests that after adjustment for the number of predictors, the explanatory power is still comparable to the R-squared, which is indicative of a model that is not overly complex for the amount of data.

AIC (9479.3): The Akaike Information Criterion helps to compare models with a trade-off between the goodness of fit and model complexity. A lower AIC value generally indicates a better model. An AIC of 9479.3 for this model provides a baseline for comparison with other models.

For model6 (with interaction terms):

R-squared (0.127): A slightly higher R-squared value than model5 suggests that the addition of interaction terms offers a marginal increase in explaining the variance in authoritarian values.

Adjusted R-squared (0.121): Despite the slight increase in R-squared, the adjusted R-squared is lower than that of model5, indicating that the addition of interaction terms does not improve the explanatory power of the model sufficiently when considering the penalty for additional predictors.

AIC (9484.5): The AIC is higher than that of model5, suggesting that the added complexity of the interaction terms may not be justifiable in terms of model fit. Typically, we prefer models with lower AIC values, so in this comparison, model5 might be more preferable.

Task 4

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

netherlands_data <- read_fst("netherlands_data.fst")

df <- netherlands_data
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))
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
    )
  )
df <- df %>%
  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(df$gen)
## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         3771         6147         4548         2730
model7 <- lm(auth ~ + cohort + polID + cohort*polID, data = df, weights = weight)
interact_plot(model7, pred = cohort, modx = polID, jnplot = TRUE)

Interpretation:

The slopes for the Left, Moderate, and Right political ideologies are not parallel, suggesting that the relationship between cohort and auth scores differs by political ideology.

It seems that for the Left and Moderate groups, as the cohort gets younger (moving towards the right on the x-axis), the auth score decreases, indicating less authority or conservatism. However, the Left group starts with a higher auth score and decreases more sharply.

For the Right group, the slope is relatively flatter, indicating a more stable auth score across generations.

In summary, generational changes in attitudes towards authority or conservatism are more pronounced among the Left and Moderate groups than among the Right group.

Task 5

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

netherlands_data <- read_fst("netherlands_data.fst")

df <- netherlands_data
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))
df <- df %>%
  mutate(
   religious = case_when(
      rlgdgr %in% c(77, 88, 99) ~ NA_real_,
      TRUE ~ rlgdgr
    )
  )
df <- df %>%
  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(df$ID)
## 
##  Left Right 
##  5449  7010
model8 <- lm(auth ~ religious + ID + religious*ID, data = df, weights = weight)
interact_plot(model8, pred = religious, modx = ID, jnplot = TRUE)

Interpretation:

The slopes for the Left and Right are not parallel, suggesting that the effect of religious on auth varies by political orientation. Specifically, it seems that the slope is steeper for the Right, which could imply that increases in religiousness are associated with a greater increase in auth for those on the political right compared to those on the left.