Basa_Luis_Homework_5
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.5
## ✔ 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## Warning: package 'flextable' was built under R version 4.3.3##
## 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#
germany_data <- read_fst("germany_data.fst")gd <- germany_datagd <- gd %>%
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))
gd$auth <- scales::rescale(gd$behave +
gd$secure +
gd$safety +
gd$tradition +
gd$rules, to=c(0,100), na.rm = TRUE)
gd <- gd %>%
filter(!is.na(auth))gd <- gd %>%
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 = gd, weights = weight)
model4 <- lm(auth ~ polID + religious + polID*religious, data = gd, weights = weight)modelsummary(
list(model3, model4),
fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL)| (1) | (2) | |
|---|---|---|
| (Intercept) | 57.5 (0.7)*** | 59.6 (0.9)*** |
| polIDModerate | 5.3 (0.7)*** | 2.6 (1.1)* |
| polIDRight | 9.3 (1.1)*** | 3.4 (1.9)+ |
| religious | 0.3 (0.1)** | −0.3 (0.2)+ |
| polIDModerate × religious | 0.8 (0.2)** | |
| polIDRight × religious | 1.5 (0.4)*** | |
| Num.Obs. | 2855 | 2855 |
| R2 | 0.035 | 0.041 |
| R2 Adj. | 0.034 | 0.039 |
| AIC | 24680.0 | 24665.5 |
| BIC | 24709.7 | 24707.2 |
| Log.Lik. | −12334.978 | −12325.764 |
| RMSE | 17.13 | 17.13 |
#INTERPRETATION#
MLR:
Firstly, it is important to note that all levels of the MLR model indicate a p-value < .0001, meaning that there is less than a 0.1% chance that the obtained value be due to chance; in other words that the obtained values are statistically significant.
In relation to our Intercept which represents our omitted reference category which in our case is the “Left” political identification. Thus, in terms of its relationship with the “Moderate” political identification, there is a positive and very strong relationship between the two variables. Further. We can further interpret this as indicating that as individuals who identify as politically moderate are 5.3 units higher compared to those who identify as “Left. Political”Right” also have a positive and extremely strong relationship. Additionally, individuals who identify as such are 9.3 units larger than those who identify as “Left”.
Religiosity is also given a coefficient of 0.3, which represents our rate of increase per unit of 1.
Interaction:
Regarding the interaction between PolID x Religious, we see a moderate to substantial relationship as individuals identify more further to the right. As well, there is a negative relationship between both relationships. As well, we see a rate of increase by (0.3) with moderate political identification and (0.5) with political Rightness.
#TASK 2#
interaction_plot <- effect("polID*religious", model4, na.rm = TRUE)
plot(interaction_plot,
main = "Interaction Effect",
xlab = "Religiousness",
ylab = "Authoritarian attitudes scale")
#INTERPRETATION#
Based on the slopes of all interaction plots, they seem to have a point of intersection if continued to infinity, which indicates to us that there exists a relationship between political identification and authoritarian scores. Political “Right” also seem to have the steepest slope which may suggest a strong interaction between political rightness and authoritarian scores. Likewise, political leftness shows a negative correlation and a moderately even slope which indicates an interaction much less and less predictive of authoritarian scales than political rightness.
#TASK 3#
netherlands_data <- read_fst("netherlands_data.fst")nd <- netherlands_dataVariables : auth, religion, ID, gen, religion*gen
#auth
nd <- nd %>%
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 ))
nd$auth <- scales::rescale(nd$behave +
nd$secure +
nd$safety +
nd$tradition +
nd$rules, to=c(0,100), na.rm=TRUE)
nd <- nd %>%
filter(!is.na(auth))#religion
nd <- nd %>%
mutate(religion = case_when(
rlgblg == 2 ~ "No",
rlgblg == 1 ~ "Yes",
rlgblg %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(rlgblg)
))
table(nd$religion)##
## No Yes
## 10913 6750#ID
nd <- nd %>%
mutate(
ID = case_when(
lrscale >= 0 & lrscale <= 4 ~ "Left",
lrscale >= 6 & lrscale <= 10 ~ "Right",
lrscale > 10 ~ NA_character_
) )
table(nd$ID)##
## Left Right
## 5449 7010#gen
nd <- nd %>%
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)
),
gen = factor(gen,
levels = c("1", "2", "3", "4"),
labels = c("Interwar", "Baby Boomers", "Gen X", "Millenials"))
)
table(nd$gen)##
## Interwar Baby Boomers Gen X Millenials
## 3771 6147 4548 2730model5 <- lm(auth ~ religion + ID + gen, data = nd, weights = weight)
model6 <- lm(auth ~ religion + ID + gen + religion*gen, data = nd, 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) *** |
| genMillenials | −8.6 (1.5) *** | −8.6 (1.9) *** |
| religionYes × genBaby Boomers | 1.6 (2.6) | |
| religionYes × genGen X | 0.5 (2.8) | |
| religionYes × genMillenials | −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#
Interpretation MLR w/o Interaction (Model 5)
Firstly, at all levels of our independent variables seem to have statistical findings of p < .0001 (***). This indicates a 0.1% chance of the obtained value to be a result of chance. For the variable of religionYes, we see a positive and extremely strong relationship between that and those who answer as religionNO which is in our reference category and that there is an 8.2 unit average in authoritarian scales.
Individuals who identify as politically right also show an average higher than those in the reference category, which in our case is the politically left identification. With a coefficient of 0.9 which indicates a positive and extremely strong relationship between political identification and authoritarian scales.
Finally, all levels of generation other than the “Interwar” show negative relationships since it was placed in the reference category and we are attempting to compare other generations and their interaciton with the authoratitive scale given that all other variables are held constant. We see again that all generations also show an extremely strong relationship between that and authoritarian scores. Further, comparing to the “Interwar” reference group, it seems that as generation groups increase, there is also an inverse reaction to authoritarian scales which tells us that as generations transition further, individuals become less authoritarian.
Interpretation MLR w/Interaction (Model 6)
We see an inverse reaction when individuals of generations also identify as religious as well. As the generations progress, the relationship seems to lean towards a negative relationship. Which tells us that as the years progress, authoritarian scales seem to go down.
#TASK 4#
Variables: cohort, polID
nd <- nd %>%
mutate(
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
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 = nd, weights = weight)
interact_plot(model7, pred = cohort, modx = polID, jnplot = TRUE)
#INTERPRETATION#
This graph illustrates the relationship between authoritarian scales or individuals born in indicated cohorts based on their political affiliation for the country of Netherlands. Individuals who identify as politically right seem to have the steepest slope amongst the three which indicates that individuals who identify as politically right strongly lower as the years progress. Similarly, for individuals who identify as politically left, have a more gentle slope, which tells us that authoritarian scales for the Dutchmen who identify as politically left gradually decline as the years progress. Finally, for politically moderate individuals, we see a particularly moderate slope but anyways the same type of interaction between authoritarian scales and year for other political identifications, just a much less decline in authoritarian scales than those of the political right but steeper than those of the politically left.
#TASK 5#
netherlands_data <- read_fst("netherlands_data.fst")nd <- netherlands_datand <- nd %>%
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 ))
nd$auth <- scales::rescale(nd$behave +
nd$secure +
nd$safety +
nd$tradition +
nd$rules, to=c(0,100), na.rm=TRUE)
nd <- nd %>%
filter(!is.na(auth))
table(nd$auth)##
## 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60
## 8 10 17 43 42 87 113 200 326 423 529 757 1000 1133 1467 1640
## 64 68 72 76 80 84 88 92 96 100
## 1734 1803 1605 1472 1380 709 514 311 202 165nd <- nd %>%
mutate(
religious = case_when(
rlgdgr %in% c(77, 88, 99) ~ NA_real_,
TRUE ~ rlgdgr
)
)
table(nd$religious)##
## 0 1 2 3 4 5 6 7 8 9 10
## 3348 988 1094 1129 923 1860 2168 2692 1996 718 728nd <- nd %>%
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_
))
table(nd$ID)##
## Left Right
## 5449 7010model8 <- lm(auth ~ religious + ID + religious*ID, data = nd, weights = weight)
interact_plot(model8, pred = religious, modx = ID, jnplot = TRUE)
#INTERPRETATION#
Here, we see the interaction between authoritarian scales of individuals from Netherlands who identify as politically Left or Right while ommiting for moderate political affiliation. Already, we see that both political affiliations seem to follow a similar trend and relationship, indicating that an increase in the religiousity scale also results in an increase of authoritarian scales regardless of political affiliation. However, these lines may seem to be parallel however one seems to be slightly askewed than the other so I would not be quick to conclude a lack of a relationship between the two but instead would conclude that there is a little to almost no relationship between the two.