Homework 5
Luna Marin
2024-03-25
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.
rm(list=ls()); gc()## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 526472 28.2 1169497 62.5 NA 669417 35.8
## Vcells 970375 7.5 8388608 64.0 16384 1851612 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"# Loading Germany's data in
germany_data <- read_fst("germany_data.fst")# Renaming the dataset to df
df <- germany_data# Survey weights
df$weight <- df$dweight * df$pweight
survey_design <- svydesign(ids = ~1, data = df, weights = ~weight)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 |
Interpreting Coefficients (No Interaction):
A coefficient of 4.6 for polIDModerate suggests that, all else being equal in terms of religiousness, individuals identifying as politically moderate have an average outcome 4.6 units higher compared to those identifying as left-wing. Similarly, a coefficient of 7.7 for polIDRight indicates that, controlling for religiousness, individuals identifying as politically right have an average outcome 7.7 units higher compared to the left-wing reference category. Regarding the “religious” predictor, a coefficient of 0.9 implies that for every one-unit increase in religiousness, there is an expected average increase of 0.9 units on the authoritarian values scale, assuming individuals are categorized similarly in terms of ideology. All three coefficients have p-values < 0.001, indicating high statistical significance.
Interpreting Model Fit Metrics:
An adjusted R-Squared value of 0.048 for the first model indicates that political ideology among individuals with the same level of religiousness explains approximately 4.8% of the variance in the authoritarian values scale for Germany, after considering the predictors. While statistically significant, this suggests that the model only accounts for a small portion of the variability in authoritarian values. In comparison, the second model’s adjusted R-Squared value is similar at 0.050 (~5%), but it explains slightly more variability than the first model. AIC, BIC, and Log.Lik are utilized to gauge the goodness of fit while penalizing for complexity. The second model exhibits slightly lower AIC and BIC values, suggesting improved performance compared to the first model. Additionally, the Log.Lik is higher (less negative) for the second model, indicating a better fit. Therefore, overall, the second model is preferable for use and evaluation.
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.34211Interpretation: Upon observing the slopes of each line, it becomes evident that they are not perfectly parallel. While the lines corresponding to polID = Moderate and polID = Right appear somewhat more parallel to each other compared to the line for polID = Left, it is clear that all three lines possess distinct slopes. This discrepancy suggests the presence of an interaction effect among the variables of religiousness, political affiliation, and authoritarian value scale score. In essence, both religiousness and political affiliation exert influence on the authoritarian value scale score. This interpretation aligns with the findings from model fit metrics, where the second model (with interaction accounted for) demonstrates superior goodness of fit compared to the first model (without interaction). Notably, in Germany, higher levels of religiousness and a right-wing political ideology seem to correlate with the highest authoritarian value scale scores (75), while lower levels of religiousness and a left-wing political ideology appear to correlate with lower authoritarian value scale scores (>60).
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")
nf <- netherlands_datanf$weight <- nf$dweight * nf$pweight
survey_design <- svydesign(ids = ~1, data = nf, weights = ~weight)nf <- nf %>%
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
nf$auth <- scales::rescale(nf$behave +
nf$secure +
nf$safety +
nf$tradition +
nf$rules, to=c(0,100), na.rm=TRUE)
nf <- nf %>% filter(!is.na(auth))nf <- nf %>%
mutate(religion = case_when(
rlgblg == 2 ~ "No",
rlgblg == 1 ~ "Yes",
rlgblg %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(rlgblg)
))
# check
table(nf$religion)##
## No Yes
## 10913 6750nf <- nf %>%
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(nf$ID)##
## Left Right
## 5449 7010nf <- nf %>%
mutate(
age = ifelse(agea == 999, NA, agea),
# Recoding cohort variable, setting years before 1930 and after 2000 to NA
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
# Recoding generational cohorts based on year of birth
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(cohort) # Keeping other values as character if they do not fit the ranges
),
# Converting cohort to a factor with labels
gen = factor(gen,
levels = c("1", "2", "3", "4"),
labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
)
table(nf$gen)##
## Interwar Baby Boomers Gen X Millennials
## 3771 6147 4548 2730model5 <- lm(auth ~ religion + ID + gen, data = nf, weights = weight)
model6 <- lm(auth ~ religion + ID + gen + religion*gen, data = nf, 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 |
Interpreting Coefficients (No Interaction):
In the absence of interaction effects, the coefficient for ReligionYes stands at 7.2, indicating that, on average, religious individuals score 7.2 units higher on the authoritarian values scale compared to their non-religious counterparts. Similarly, IDRight exhibits a coefficient of 4.2, suggesting that, on average, individuals leaning towards right-wing political ideology score 4.2 units higher on the authoritarian values scale than those leaning towards left-wing political ideology. Regarding generational differences, Baby Boomers score approximately 5.2 units lower, Gen Xers score around 7.7 units lower, and Millennials score roughly 8.4 units lower on the authoritarian values scale compared to the Interwar generation (the reference category). These findings imply an overall decline in authoritarian attitudes across generations in the Netherlands. However, identifying as politically right-wing or religious is associated with higher scores on the authoritarian values scale. All coefficients exhibit three asterisks, indicating high statistical significance (p < 0.001), with a less than 0.1% probability of the observed results being due to chance.
Interpreting Model Fit Metrics:
The adjusted R-Squared value for the first model stands at 0.124, indicating that approximately 12.4% of the variance in authoritarian value scale scores is explained by the variables included (generation, religion, and political ideology). In comparison, the adjusted R-Squared value for the second model is 0.125, suggesting that about 12.5% of the variance in authoritarian value scale scores is explained. Although the difference between these percentages is minimal, the second model demonstrates slightly better goodness of fit. The AIC is slightly lower for the second model, while the BIC is slightly higher. Moreover, the Log.Lik value is slightly higher (less negative) for the second model. Overall, these metrics indicate that the second model is preferable for evaluation, albeit the difference may stem from the complexity disparity between the models.
Task 4
Produce the model7 interaction plot from the tutorial, but again using the Netherlands data. Interpret.
netherlands_data <- read_fst("netherlands_data.fst")
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.34211model7 <- lm(auth ~ + cohort + ID + cohort*ID, data = nf, weights = weight)
interact_plot(model7, pred = cohort, modx = ID, jnplot = TRUE)
Interpretation: The inclinations of the two lines diverge, with the slope for ID = Right (~30 degrees) displaying a steeper inclination compared to ID = Left (~20 degrees). This suggests the presence of an interaction effect among the variables of political affiliation, generation (cohort), and authoritarian value scale score, wherein both cohort and political affiliation contribute to influencing the authoritarian value scale score. In the context of the Netherlands, it appears that later generations (e.g., millennials) identifying with left-wing ideology yield the lowest authoritarian value scale scores (~55), while earlier generations (boomers) endorsing right-wing ideology yield the highest authoritarian value scale scores (~70+).
Task 5
Produce the model8 interaction plot from the tutorial, but again using the Netherlands data. Interpret.
nf <- nf %>%
mutate(
religious = case_when(
rlgdgr %in% c(77, 88, 99) ~ NA_real_,
TRUE ~ rlgdgr
)
)
model8 <- lm(auth ~ religious + ID + religious*ID, data = nf, weights = weight)
interact_plot(model8, pred = religious, modx = ID, jnplot = TRUE)
Interpretation: Upon observation of the graph, the inclinations of the two lines corresponding to ID = Left and ID = Right appear largely parallel. Although ID = Right exhibits a slightly steeper angle compared to ID = Left, their similarities are notable. This suggests that in the Netherlands, there appears to be no discernible interaction effect between religiousness, political affiliation, and authoritarian value scale score. On average, individuals identifying as politically right-wing exhibit slightly higher authoritarian value scale scores (~60-75) compared to those identifying as left-wing (~55-65). However, the presence of religiousness does not seem to significantly impact authoritarian value scale scores. This distinction contrasts with findings from Germany, where religiousness exhibited an interaction effect on authoritarian value scale scores.