Homework 5
Christian Dizon
2024-4-25
Initial set-up
# 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## Warning: package 'infer' was built under R version 4.3.3## Warning: package 'effects' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Warning: package 'survey' was built under R version 4.3.3## 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## Warning: package 'aod' was built under R version 4.3.3##
## Attaching package: 'aod'
##
## The following object is masked from 'package:survival':
##
## rats## Warning: package 'interactions' was built under R version 4.3.3##
## 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"Homework 5 (5%): due by next lecture on March 26th
Instructions: Start a new R markdown for the homework and call it “Yourlastname_Firstname_Homework_5”.
Copy everything below from Task 1 to Task 5. Keep the task prompt and questions, and provide your code and answer underneath.
Remember: you need all the steps for your code to work, including loading your data – otherwise it will not knit.
To generate a new code box, click on the +C sign above. Underneath your code, provide your answer to the task question.
When you are done, click on “Knit” above, then “Knit to Html”. Wait for everything to compile. If you get an error like “Execution halted”, it means there are issues with your code you must fix. When all issues are fixed, it will prompt a new window. Then click on “Publish” in the top right, and then Rpubs (the first option) and follow the instructions to create your Rpubs account and get your Rpubs link for your document (i.e., html link as I provide for the tutorial).
Note: Make sure to provide both your markdown file and R pubs link. If you do not submit both, you will be penalized 2 pts. out of the 5 pts. total.
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 <- germany_datagermany_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 (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 3.9 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.
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 and we incorrectly reject the null, under frequentist assumptions.
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.34211There is a positive correlation of religion in all graphs. The polID Left however, does not seem to follow the same trends as polID Moderate and 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")
netherlands_data <- netherlands_datanetherlands_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$trstplt <- ifelse(netherlands_data$trstplt > 10, NA, netherlands_data$trstplt)
netherlands_data$trstprl <- ifelse(netherlands_data$trstprl > 10, NA, netherlands_data$trstprl)
netherlands_data$trstprt <- ifelse(netherlands_data$trstprt > 10, NA, netherlands_data$trstprt)
netherlands_data$trust <- scales::rescale(netherlands_data$trstplt + netherlands_data$trstprl + netherlands_data$trstprt, na.rm = TRUE, to = c(0, 100))
netherlands_data$populist <- scales::rescale(netherlands_data$trust, na.rm = TRUE, to=c(100,0))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_
),
religious = case_when(
rlgdgr %in% c(77, 88, 99) ~ NA_real_,
TRUE ~ rlgdgr
)
)netherlands_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)
))netherlands_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 7010netherlands_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 2730model5 <- 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 |
religiousYES: On average, religious individuals score 7.2 points higher on the dependent variable than other categories.
IDright: People that are politically RIGHT, scored 4.2 points higher than those within the three generational categories.
genBaby Boomers (-5.2), genGen X (-7.7), genMillenials (-8.4): All three generational categories scored the lowest, with Millenials scoring the lowest of -8.4 on the dependent variable.
Task 4
Produce the model7 interaction plot from the tutorial, but again using the Netherlands data. Interpret.
model7 <- lm(auth ~ + cohort + polID + cohort*polID, data = netherlands_data, weights = weight)
interact_plot(model7, pred = cohort, modx = polID, jnplot = TRUE)
## Left (Blue Line): Individual that were identified as politically LEFT
has a lower, downward trend in authoritarian attitudes over generations.
This suggests that over time, younger generations identifying as
politically LEFT has lower authoritarian score compared the older
generation
Task 5
Produce the model8 interaction plot from the tutorial, but again using the Netherlands data. Interpret.
model8 <- lm(auth ~ religious + ID + religious*ID, data = netherlands_data, weights = weight)
interact_plot(model8, pred = religious, modx = ID, jnplot = TRUE)
## People that identify as politicaly left: There seems to be a positive
relationship between religion and authoritarian perspective, with both
correlating in an upward pattern.
End
Appendix
If you want some examples of ESS papers and how they talk about different elements like data and findings, here are some links:
Food Comes First, Then Morals: Redistribution Preferences, Parochial Altruism, an Immigration in Western Europe https://www.journals.uchicago.edu/doi/pdf/10.1086/694201?casa_token=eUhZ0HbtlD4AAAAA:SGZ0609PlTu7pmc-8cE-ztOIaCUfD9pbkUWlMrTUlYZP-32wRLHQQfRA7sT8_WcwCp39mE6C2BnI
Cultural Backlash? How (Not) to Explain the Rise of Authoritarian Populism https://www.cambridge.org/core/services/aop-cambridge-core/content/view/FFE9742798D8CC4BF6ED325FDBAFA251/S0007123421000363a.pdf/div-class-title-cultural-backlash-how-not-to-explain-the-rise-of-authoritarian-populism-div.pdf
Party membership and closeness and the development of trust in political institutions: An analysis of the European Social Survey, 2002–2010 https://journals.sagepub.com/doi/abs/10.1177/1354068813509519?casa_token=cWvN2GK3tjYAAAAA%3A88ipxKnceH2FJ0suJsd7Q5RzK3IfsCOsqARBVonjnQnUuIj7UZR0hWDebi4fxbBdtFddpiLHFCebjw&journalCode=ppqa
Deservingness perceptions, welfare state support and vote choice in Western Europe https://www-tandfonline-com.myaccess.library.utoronto.ca/doi/full/10.1080/01402382.2020.1715704?casa_token=c0kCwK30vRAAAAAA%3AQgU2EcX6gggkMbguBnSNUdb-YeLTu89ath6GuWppbBGZfIHEGfAl7zG269mPl-aIA6FjpQ1E8FgCXg
Increasingly unequal? The economic crisis, social inequalities and trust in the European Parliament in 20 European countries https://ejpr.onlinelibrary.wiley.com/doi/pdf/10.1111/1475-6765.12126?casa_token=H3MiTpwFHA0AAAAA:yP3OQdvt26J46rojOyB_kx5Y4x4Fsg6vI56Z7P6CKm4-8APdmnmLxY5wgNcwimB97CbvnNQOP_ueUx4