Markiewicz_Alexandra_Homework_5
2024-03-25
# 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")
netherlands_data <- read_fst("netherlands_data.fst")germany_data$weight <- germany_data$dweight * germany_data$pweight
survey_design <- svydesign(ids = ~1, data = germany_data, weights = ~weight)
netherlands_data$weight <- netherlands_data$dweight * netherlands_data$pweight
survey_design <- svydesign(ids = ~1, data = netherlands_data, weights = ~weight)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 <- germany_data %>%
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 ))
germany_data$auth <- scales::rescale(germany_data$behave +
germany_data$secure +
germany_data$safety +
germany_data$tradition +
germany_data$rules, to=c(0,100), na.rm=TRUE)
germany_data <- germany_data %>% filter(!is.na(auth))germany_data <- germany_data %>%
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 = germany_data, weights = weight)
model4 <- lm(auth ~ polID + religious + polID*religious, data = germany_data, 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 |
A coefficient of 4.6 for ‘polIDModerate’ suggests that, when comparing individuals with the same religiousness, the average outcome for individuals identifying as politically moderate is 4.6 units higher compared to the reference category (politically left) on the authoritarian values scale.
A coefficient of 7.7 for ‘polIDRight’ suggests that, when comparing individuals with the same religiousness, the average outcome for individuals identifying as politically right is 7.7 units higher compared to the reference category (politically left) on the authoritarian values scale.
A coefficient of 0.9 for ‘religious’ suggests that, when holding ideological self-placement constant, every one-unit increase in religiousness is expected to result in an average increase of 0.9 units on the authoritarian values scale.
All coefficients have ***, indicating a p < 0.001. This suggests that there is less than a 0.1% probability that the observed result is due to chance, leading to an incorrect rejection of the null.
The intercept of 55.3 represents the expected value of the outcome variable when all the predictor variables are set to their reference levels. For continuous variables, this means setting them to zero, and for categorical variables, it refers to the level that is excluded from the model.
The R-squared and adjusted R-squared show the smallest increase in the second model, while the AIC and BIC show the smallest decrease. This suggests that the additional interaction term (polID x religious) explains only a marginal proportion of the variance in the outcome outside of what is already explained by the main effects (polID and “religious”), implying a marginal improvement in model fit.
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.34211The x-axis for the interaction plot is “Religiousness,” the y-axis is “Authoritarian attitudes scale,” and the grouping variable is “polID.” The “Left,” “Moderate,” and “Right” political ideologies would not intersect, meaning the lines are parallel and there is no interaction effect. This suggests the interaction between political ideology and religiousness on authoritarianism is not stable across all political ideologies, implying the impact of religiousness on authoritarianism differs depending on an individual’s 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 <- netherlands_data %>%
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 ))
netherlands_data$auth <- scales::rescale(netherlands_data$behave +
netherlands_data$secure +
netherlands_data$safety +
netherlands_data$tradition +
netherlands_data$rules, to=c(0,100), na.rm=TRUE)
netherlands_data <- netherlands_data %>% filter(!is.na(auth))netherlands_data <- netherlands_data %>%
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(netherlands_data$gen)##
## Interwar Baby Boomers Gen X Millennials
## 3771 6147 4548 2730netherlands_data <- netherlands_data %>%
mutate(religion = case_when(
rlgblg == 2 ~ "No",
rlgblg == 1 ~ "Yes",
rlgblg %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(rlgblg)
))
table(netherlands_data$religion)##
## No Yes
## 10913 6750netherlands_data <- netherlands_data %>%
mutate(ID = case_when(
lrscale >= 0 & lrscale <= 4 ~ "Left",
lrscale >= 6 & lrscale <= 10 ~ "Right",
lrscale > 10 ~ NA_character_,
TRUE ~ NA_character_
))
table(netherlands_data$ID)##
## Left Right
## 5449 7010model5 <- lm(auth ~ religion + ID + gen, data = netherlands_data, weights = weight)
model6 <- lm(auth ~ religion + ID + gen + religion*gen, data = netherlands_data, 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 |
A coefficient of 7.2 for ‘religionYes’ suggests that, when holding other variables constant, individuals who identify with a religion have an average increase of 7.2 units on the authoritarian values scale compared to those who do not identify with any religion, implying religiousness is positively associated with higher levels of authoritarian values in the Netherlands.
A coefficient of 4.2 for ‘IDRight’ suggests that, when holding other variables constant, individuals who identify as politically right have an average increase of 4.2 units on the authoritarian values scale compared to those who identify as politically left, implying being politically right is positively associated with higher levels of authoritarian values in the Netherlands.
A coefficient of -5.2 for ‘genBaby Boomers’ suggests that, when holding other variables constant, Baby Boomers have an average decrease of 5.2 units on the authoritarian values scale compared to the Interwar generation, implying Baby Boomers have lower levels of authoritarian values compared to the Interwar generation in the Netherlands.
A coefficient of -7.7 for ‘genGen X’ suggests that, when holding other variables constant, Gen X has an average decrease of 7.7 units on the authoritarian values scale compared to the Interwar generation, implying Gen X has lower levels of authoritarian values compared to the Interwar generation in the Netherlands.
A coefficient of -8.4 for ‘genMillennials’ suggests that, when holding other variables constant, Millennials have an average decrease of 8.4 units on the authoritarian values scale compared to the Interwar generation, implying Millennials have lower levels of authoritarian values compared to the Interwar generation in the Netherlands.
All coefficients have ***, indicating a p < 0.001. This suggests that there is less than a 0.1% probability that the observed result is due to chance, leading to an incorrect rejection of the null.
The R-squared, adjusted R-squared and BIC show the smallest increase in the second model, while the AIC shows the smallest decrease. This suggests that the additional interaction term (religionYes x gen) explains only a marginal proportion of the variance in the outcome outside of what is already explained by the main effects (“religionYes” and gen), implying a marginal improvement in model fit.
Task 4
Produce the model7 interaction plot from the tutorial, but again using the Netherlands data. Interpret.
netherlands_data <- netherlands_data %>%
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_
)
)model7 <- lm(auth ~ + cohort + polID + cohort*polID, data = netherlands_data, weights = weight)
interact_plot(model7, pred = cohort, modx = polID, jnplot = TRUE)
The x-axis for the interaction plot is “cohort,” the y-axis is “auth,” and the grouping variable is “polID.” The “Moderate” and “Right” political ideologies intersect, meaning the lines are not parallel, and there is an interaction effect. This suggests the effect of cohort on authoritarianism varies depending on the political ideology, implying the relationship between cohort and authoritarianism differs between individuals who are politically moderate and right. However, the “Left” political ideology does not intersect with either. This suggests there is no interaction effect between cohort and political ideology for individuals who are politically left, implying the relationship between cohort and authoritarianism is consistent across different cohorts for individuals who are politically left.
Task 5
Produce the model8 interaction plot from the tutorial, but again using the Netherlands data. Interpret.
netherlands_data <- netherlands_data %>%
mutate(
religious = case_when(
rlgdgr %in% c(77, 88, 99) ~ NA_real_,
TRUE ~ rlgdgr
)
)model8 <- lm(auth ~ religious + ID + religious*ID, data = netherlands_data, weights = weight)
interact_plot(model8, pred = religious, modx = ID, jnplot = TRUE)
The x-axis for the interaction plot is “religious,” the y-axis is “auth,” and the grouping variable is “ID.” The “Left” and “Right” political ideologies do not intersect, meaning the lines are parallel, and there is no interaction effect. This suggests the interaction between political ideology and religiousness on authoritarianism is not stable across either political ideology, implying the impact of religiousness on authoritarianism differs depending on an individual’s political ideology.