RYLEE_BOCCONGELLE_SHAW_HOMEWORK5
Rylee Boccongelle
2024-03-21
# 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"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")germany_data$weight <- germany_data$dweight * germany_data$pweight
survey_design <- svydesign(ids = ~1, data = germany_data, weights = ~weight)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))) %>%
# Apply the reverse coding
mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))
# Now you can calculate 'schwartzauth' after the NA recoding
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 |
PolIDmoderate model 3 - on avwrage, those on the moderate of the political spectrum who are religious score 4.6 points higher on the authoritarian scale than those on the left on the authroiratian scale.
PolIDright model 3 -on avaerage, those on the right of the political spectrum who are religious score 7.7 points higher on the authoritarian scale than those on the left on the authoritarian scale
religious model 3 - for every one unit increase of the religious variable, there is a 0.9 point increase on the authoritarian scale on average.
Comparing Model 3 and 4 - the r squre value of model 3 is 0.048 meaning that the explanatory variable explains 4.8 percent of the variance in the outcome variables included in model 3. Since the r squared value is only 4.8 percent whilst it is statistically signifigant it indicates that there are other variables outside of the ones measured in this model that explain these attitudes. On the other hand, model 4’s r sqaured value of 0.050 insicates that the explanatory variable explains 5 percent of the variace of the outcome variables in model 4. since the r sqaured value is low like i model 3, this tells us that there are other variables to consider that explain these attitudes. Moreover, model 4 provides a slightly more holistic view.
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.34211interaction## function (..., drop = FALSE, sep = ".", lex.order = FALSE)
## {
## args <- list(...)
## narg <- length(args)
## if (narg < 1L)
## stop("No factors specified")
## if (narg == 1L && is.list(args[[1L]])) {
## args <- args[[1L]]
## narg <- length(args)
## }
## for (i in narg:1L) {
## f <- as.factor(args[[i]])[, drop = drop]
## l <- levels(f)
## if1 <- as.integer(f) - 1L
## if (i == narg) {
## ans <- if1
## lvs <- l
## }
## else {
## if (lex.order) {
## ll <- length(lvs)
## ans <- ans + ll * if1
## lvs <- paste(rep(l, each = ll), rep(lvs, length(l)),
## sep = sep)
## }
## else {
## ans <- ans * length(l) + if1
## lvs <- paste(rep(l, length(lvs)), rep(lvs, each = length(l)),
## sep = sep)
## }
## if (anyDuplicated(lvs)) {
## ulvs <- unique(lvs)
## while ((i <- anyDuplicated(flv <- match(lvs,
## ulvs)))) {
## lvs <- lvs[-i]
## ans[ans + 1L == i] <- match(flv[i], flv[1:(i -
## 1)]) - 1L
## ans[ans + 1L > i] <- ans[ans + 1L > i] - 1L
## }
## lvs <- ulvs
## }
## if (drop) {
## olvs <- lvs
## lvs <- lvs[sort(unique(ans + 1L))]
## ans <- match(olvs[ans + 1L], lvs) - 1L
## }
## }
## }
## structure(as.integer(ans + 1L), levels = lvs, class = "factor")
## }
## <bytecode: 0x13fd222b0>
## <environment: namespace:base>- When interpreting the output, it can be seen that all three political views, left, moderate, and right have postive slopes to the variable of religion with different degrees. The left variable has the least postive slope however is still postive slope. In constrast right political views illustrate that people who identiy right on the political spectrum on average, hold higher religious views in comparision than those moderate and left on the political spectrum.
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")netherlands_data$weight <- netherlands_data$dweight * netherlands_data$pweight
survey_design <- svydesign(ids = ~1, data = netherlands_data, weights = ~weight)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))) %>%
# Apply the reverse coding
mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))
# Now you can calculate 'schwartzauth' after the NA recoding
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),
# 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(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)
))
# check
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_, # Set values above 10 as NA
TRUE ~ NA_character_ # Ensure value 5 and any other unexpected values are set as NA
))
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 |
- Religionyes this variable indicates that those who do identify with a religion and are part of the same political identification group, tend to be 7.2 points higher on the authoritairan scale compared to those who do not identify with a religion
IDright - this variable indicates that those who identify with a religion and identify with the right of the political spectrum on average score 4.2 points higher on the authoritarian scale compared to those on the left of the poltiical spectrum.
genBaby Boomers - on average baby boomers who are religious tend to score 5.2 points less points on the authoritarian scale than the interwar generation
genx - on average gen x who are religious score 7.7 points less on the authoritarian scale than the interwar generation
Millennials - on average millenials who are religious score 8.4 points less on the authoritarian scale than the interwar generation
Comparing models - In model 5 the r squared value is 0.124 this illustrates that the explanatory variable explains 12.4 percent of the vairence in the outcome variables. Since the r squared value is only 12.4 percent whilst it is statistically signifigant it indicates that there are other variables outside of the ones measured in this model that explain these attitudes. In comparision, in model 6 the r squared value is 0.125 this illustrates that the explanatory variable explains 12.5 percent of the vairence in the outcome variables. Since the r squared value is only 12.5 percent whilst it is signifigant it indicates that there are other vairables outside of those measured in this model that explain these attitudes. However, when comparing the 2 r sqaure values in model 5 and 6, model 6 is slightly better explains the interactions between the variables than in model 5.
Task 4
Produce the model7 interaction plot from the tutorial, but again using the Netherlands data. Interpret.
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))) %>%
# Apply the reverse coding
mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))
# Now you can calculate 'schwartzauth' after the NA recoding
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(
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_
)
)netherlands_data <- netherlands_data %>%
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(netherlands_data$gen)##
## Interwar Baby Boomers Gen X Millennials
## 3771 6147 4548 2730model7 <- lm(auth ~ + cohort + polID + cohort*polID, data = netherlands_data, weights = weight)
interact_plot(model7, pred = cohort, modx = polID, jnplot = TRUE)
- The first noticiable interpretation from this graph is that
authortarian values have overall decreased with all political
orientations but has decreased to the fewest points on the scale with
those who identify on the right of the political spectrum. When looking
at this graph it is evident that there is a stronger interaction between
authoritarian values and right and moderate political identifications
because the lines outputed are not parrelel. ## 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)
It is certainly evident from this graph that as you increase
religiosity, both those who fall on the left and the right of the
political spectrum increase on the authoritarian scale. Additionally
from this graph it can also be said that those who identify with the
left of the political spectrum have fall lower on the authoritarian
scale than those who are on the right of the political spectrum.