Kaur_Gurjyot_Homework_5
Gurjyot Kaur
2024-03-26
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"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_datadf$weight <- df$dweight * df$pweight
survey_design <- svydesign(ids = ~1, data = df, weights = ~weight)## Setting up our populist scale
# Setting values greater than 10 to NA
df$trstplt <- ifelse(df$trstplt > 10, NA, df$trstplt) # anything above 10 is NA, if not keep it as is
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)) # scale to 100, add all the columns (trst ones)
# 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)) # scale is flipped, trust is inverse of populistdf <- df %>% # auth variable
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 %>% #ID
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) # this one is the interactionmodelsummary(
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, the average outcome for individuals identifying as politically moderate is 4.6 units higher compared to the ideologically left reference category (omitted category)
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 compared to the ideologically left reference category (omitted category)
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, given that such individuals have the same ideological categorization (i.e., holding ideological self-placement constant).
All three coefficients have three asterisks relating to p values less than 0.001, suggesting sufficient evidence to reject the null hypothesis (less than a 0.1% probability that the outcome is due to pure chance).
In terms of the model fit metrics for both models, there are not any significant changes. For instance, the R-squared increased by 0.002, from 0.048 to 0.050. However, the BIC the AIC and BIC are slightly less for the second model, highlighting that the increase in variance explained may be better justified (there are less potential underlying issues), such as the Akaike Information Criterion decreased from 200594.1 to 200558.6.
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.
germany_data <- read_fst("germany_data.fst")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
)
)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:
There seems to be a slight interaction effect between authoritarian attitudes and religiousness in relation to political ideology because the slope is slightly increasing as it goes from politically left-leaning to politically right-leaning. For instance, the graph demonstrates that those who are politically right-leaning and very religious, tend to have stronger authoritarian attitudes compared to those who are politically left-leaning and very religious.
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")df2 <- netherlands_datadf2 <- df2 %>% # first variable auth
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
df2$auth <- scales::rescale(df2$behave +
df2$secure +
df2$safety +
df2$tradition +
df2$rules, to=c(0,100), na.rm=TRUE)
df2 <- df2 %>% filter(!is.na(auth))df2 <- df2 %>% # second variable
mutate(religion = case_when(
rlgblg == 2 ~ "No",
rlgblg == 1 ~ "Yes",
rlgblg %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(rlgblg)
))
# check
table(df2$religion)##
## No Yes
## 10913 6750df2 <- df2 %>% # third variable
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
))
# check
table(df2$ID)##
## Left Right
## 5449 7010df2 <- df2 %>%
mutate(
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(df2$gen)##
## Interwar Baby Boomers Gen X Millennials
## 3771 6147 4548 2730model5 <- lm(auth ~ religion + ID + gen, data = df2, weights = weight)
model6 <- lm(auth ~ religion + ID + gen + religion*gen, data = df2, 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
religiousYES (8.2): This coefficient suggests that, comparing individuals with the same political ideology and generation, the average outcome for individuals who are religious is 8.2 units higher compared to non-religious individuals (omitted category).
IDRight (3.4): This coefficient indicates that, comparing individuals with the same religiousness level and generation, individuals identifying as politically right are expected to be, on average, 3.4 units higher compared to the ideologically left reference category (omitted category).
genBaby Boomers (-5.8): A coefficient of -5.8 means that with every one-unit increase in the generation category of Baby Boomers, there is an associated expected average decrease of 5.8 units on the authoritarian values scale, assuming that we are comparing people with the same ideological position and religiousness level (i.e., holding ideological position and religiousness constant).
genGen X (-7.3): A coefficient of -7.3 means that with every one-unit increase in the generation category of Gen X, there is an associated expected average decrease of 7.3 units on the authoritarian values scale, assuming that we are comparing people with the same ideological position and religiousness level.
genMillennials (-8.6): A coefficient of -8.6 means that with every one-unit increase in the generation category of Gen X, there is an associated expected average decrease of 8.6 units on the authoritarian values scale, assuming that we are comparing people with the same ideological position and religiousness level.
All three coefficients have three asterisks relating to p values less than 0.001, suggesting sufficient evidence to reject the null hypothesis.
In terms of the model fit metrics for both models, there are not any significant changes. For instance, the R-squared increased by 0.001, from 0.126 to 0.127. However, the AIC and BIC are slightly higher for the second model, highlighting that the increase in variance explained may not justified (there is more complexity or potential underlying issues).
Task 4
Produce the model7 interaction plot from the tutorial, but again using the Netherlands data. Interpret.
netherlands_data <- read_fst("netherlands_data.fst")df3 <- netherlands_datadf3 <- df3 %>% # first variable auth
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
df3$auth <- scales::rescale(df3$behave +
df3$secure +
df3$safety +
df3$tradition +
df3$rules, to=c(0,100), na.rm=TRUE)
df3 <- df3 %>% filter(!is.na(auth))df3 <- df3 %>%
mutate(
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn)
)
table(df3$cohort)##
## 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945
## 117 143 132 136 142 170 173 148 195 204 221 209 223 267 269 247
## 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961
## 365 339 320 285 309 269 316 284 383 304 364 328 287 328 312 314
## 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977
## 347 351 342 305 333 350 356 373 362 341 329 276 266 242 267 279
## 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993
## 242 227 215 194 213 200 199 194 203 188 174 145 140 140 119 115
## 1994 1995 1996 1997 1998 1999 2000
## 100 103 88 90 76 72 61df3 <- df3 %>%
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_
))
table(df3$polID)##
## Left Moderate Right
## 3420 8646 4747model7 <- lm(auth ~ + cohort + polID + cohort*polID, data = df3, weights = weight)
interact_plot(model7, pred = cohort, modx = polID, jnplot = TRUE)
Interpretation
There is an indication of an interaction effect because one of the slopes are not parallel, specifically polID Right. This means the slope does vary across subgroups.
Task 5
Produce the model8 interaction plot from the tutorial, but again using the Netherlands data. Interpret.
netherlands_data <- read_fst("netherlands_data.fst")df4 <- netherlands_datadf4 <- df4 %>% # first variable auth
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
df4$auth <- scales::rescale(df4$behave +
df4$secure +
df4$safety +
df4$tradition +
df4$rules, to=c(0,100), na.rm=TRUE)
df4 <- df4 %>% filter(!is.na(auth))df4 <- df4 %>% # second variable
mutate(religious = case_when(
rlgdgr %in% c(77, 88, 99) ~ NA_real_,
TRUE ~ rlgdgr))
# check
table(df4$religious)##
## 0 1 2 3 4 5 6 7 8 9 10
## 3348 988 1094 1129 923 1860 2168 2692 1996 718 728df4 <- df4 %>%
mutate(
ID = 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_
))
table(df4$ID)##
## Left Moderate Right
## 3420 8646 4747model8 <- lm(auth ~ religious + ID + religious*ID, data = df4, weights = weight)
interact_plot(model8, pred = religious, modx = ID, jnplot = TRUE)
Interpretation
There is no indication of an interaction effect because the slopes are parallel, meaning the slopes do not vary across subgroups.