Su_Yizhan_homework4
Yizhan Su
2024-03-17
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "interactions","broom") # 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## Warning: package 'interactions' was built under R version 4.3.3## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools## [[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] "interactions" "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] "broom" "interactions" "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"Task 1 Graph the histogram distribution using the populist scale (as we did in the tutorial) for the country of Sweden. What do you note?
sweden_data <- read.fst("sweden_data.fst")sweden_data$trstplt <- ifelse(sweden_data$trstplt > 10, NA, sweden_data$trstplt)
sweden_data$trstprl <- ifelse(sweden_data$trstprl > 10, NA, sweden_data$trstprl)
sweden_data$trstprt <- ifelse(sweden_data$trstprt > 10, NA, sweden_data$trstprt)
# Creating and rescaling the trust variable
sweden_data$trust <- scales::rescale(sweden_data$trstplt + sweden_data$trstprl + sweden_data$trstprt, na.rm = TRUE, to = c(0, 100))
# 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
sweden_data$populist <- scales::rescale(sweden_data$trust, na.rm = TRUE, to=c(100,0))ggplot(sweden_data, aes(x = populist)) +
geom_histogram(bins = 30, fill = "blue", color = "black") +
theme_minimal() +
labs(title = "Distribution of Populist Scale for Sweden",
x = "Populist Scale",
y = "Count")## Warning: Removed 2503 rows containing non-finite values (`stat_bin()`).
The distribution is unimodal and skewed to the left, with a sharp peak
around the 25–30 range on the Populist Scale. This indicates that a
majority of the data points fall within this range, suggesting that
within the sample, the most common level of populist support or attitude
is relatively low to moderate.
Task 2 Run a linear regression with populist attitudes as the outcome, educational attainment (recoded as a binary “BA or more” and “No BA”), and using the Swedish data. Print the intercept and coefficients only and interpret. Then do a tidy model table displaying the p-value (with digits = 3) and interpret.
sweden_data <- sweden_data %>%
mutate(
ID = case_when(
lrscale %in% 0:4 ~ "Left", # This time: Values 0 to 3 in 'lrscale' are categorized as "Left"
lrscale %in% 6:10 ~ "Right", # This time: Values 6 to 10 in 'lrscale' are categorized as "Right"
lrscale %in% c(77, 88, 99) ~ NA_character_, # Values 77, 88, 99 in 'lrscale' are set as NA
TRUE ~ NA_character_ # Any other values not explicitly matched above
),
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
)
)sweden_pop_model <- lm(populist ~ ID, data = sweden_data)
coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)## (Intercept) IDRight
## 48.3188406 -0.7149274tidy(sweden_pop_model) |>
knitr::kable(digits = 3)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 48.319 | 0.265 | 182.562 | 0.000 |
| IDRight | -0.715 | 0.356 | -2.011 | 0.044 |
intercept (48.319): the expected average populist scale score would be 48.319 for the left, since it is the reference category. IDRight Coefficient (-7.5): The negative coefficient associated with the “Right” group indicates that, on average, individuals who identify as ideologically right in have populist scale scores that are 7.5 points lower than those who identify as left.
The p-value for IDRight is 0.044, which is less than the typical alpha level of 0.05, indicating that the relationship between educational attainment and populist attitudes is statistically significant.There is a statistically significant decrease in populist scale scores associated with having higher educational attainment
Task 3 Now using the Italian dataset and the authoritarian values scale (as we did in the tutorial), graph the average by survey year. Interpret.
df <- read.fst("italy_data.fst")
## Creates an authoritarian values scale based on human modules items (can also do libertarian value scale)
df <- df %>%
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$year <- NA
replacements <- c(2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020)
for(i in 1:10){
df$year[df$essround == i] <- replacements[i]
}
auth_avg <- df %>%
group_by(year) %>%
summarize(auth_avg = mean(auth, na.rm = TRUE))
plot_auth <- ggplot(auth_avg, aes(x = year, y = auth_avg)) +
geom_point(aes(color = auth_avg), alpha = 1.0) +
geom_line(aes(group = 1), color = "blue", linetype = "dashed") +
labs(title = "Authoritarian Scale Average by Year (Italy)",
x = "Survey Year",
y = "Authoritarian Scale Average") +
theme_minimal() +
theme(legend.position = "none") +
scale_y_continuous(limits = c(0, 100))
print(plot_auth)## Warning: Removed 1 rows containing missing values (`geom_point()`).## Warning: Removed 1 row containing missing values (`geom_line()`).
There is small fluctuations,this could imply that there has been no
significant shift in authoritarian attitudes among the surveyed
population in Italy during these years.
Task 4 Now run a linear regression model, using the modelsummary package, with authoritarian values as the outcome, and generations as the predictor/potential explanatory variable, and using the Italian data again. Rename the coefficients using the coef_rename function. Interpret the coefficients, whether they are statistically significant, and the adjusted R-squared.
italy_data <- df %>%
mutate(
# Recoding gender
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
# Recoding education
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
),
# Handle NAs for education levels
edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
# Explicitly making 'No BA' the reference category
educ.ba = factor(educ.ba, levels = c("No BA", "BA or more")),
# Recoding age, setting 999 to NA
age = ifelse(agea == 999, NA, agea),
# Recoding cohort variable, setting years before 1930 and after 2000 to NA
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
# Recoding generational cohorts based on year of birth
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(cohort) # Keeping other values as character if they do not fit the ranges
),
# Converting cohort to a factor with labels
gen = factor(gen,
levels = c("1", "2", "3", "4"),
labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
)models <- list(
"Model 1" = lm(auth ~ gen, data = italy_data),
"Model 2" = lm(auth ~ gndr, data = italy_data),
"Model 3" = lm(auth ~ educ.ba, data = italy_data))modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_rename = c("genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials", "gndrMale" = "Male", "educ.baBA or more" = "BA or More"),
title = 'Regression models for authoritarian attitudes')| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 75.3 (0.5)*** | 74.4 (0.2)*** | 73.3 (0.2)*** |
| Baby Boomers | −1.1 (0.6)* | ||
| Gen X | −2.1 (0.6)*** | ||
| Millennials | −3.6 (0.6)*** | ||
| Male | −2.6 (0.3)*** | ||
| BA or More | −1.1 (0.4)* | ||
| Num.Obs. | 7631 | 8432 | 8432 |
| R2 | 0.006 | 0.007 | 0.001 |
| R2 Adj. | 0.006 | 0.007 | 0.001 |
| AIC | 62934.0 | 69523.1 | 69579.5 |
| BIC | 62968.7 | 69544.3 | 69600.6 |
| Log.Lik. | −31462.009 | −34758.566 | −34786.739 |
| RMSE | 14.94 | 14.93 | 14.98 |
genBaby Boomers (-1.1): Baby Boomers exhibit scores that are, on average, 1.1 points lower on the authoritarian values scale (0-100). genGen X (-2.1 ): Generation X shows an even more pronounced decrease, with scores 2.1 points lower than the Interwar generation on the authoritarian values scale. genMillennials (-3.6): Millennials exhibit the most considerable decrease in authoritarian values, with their scores being 3.6 points lower than those of the Interwar generation. This marked difference is statistically significant, indicating that Millennials diverge most strongly from the Interwar generation in comparing their authoritarian values.
Adjusted R-Squared (Adj. R² = 0.006):
The adjusted R-squared value of 0.006 means that the ‘gen’ variable explain approximately 0.6% of the variance in the authoritarian values scale for italy. While statistically significant, this suggests that the model explains only a small portion of the variability in authoritarian values, pointing to other factors outside generational differences that may play critical roles in shaping these attitudes.
Task 5 Now visualize, using the modelplot() function from the modelsummary package, coefficient estimates and 95% confidence intervals for the following three models:
Model 1 = populist attitudes scale as the outcome; left/right as the single predictor (omitting moderates as we did in the tutorial), and using data for Greece.
Model 2 = populist attitudes scale as the outcome; gender as the single predictor, and using data for Greece.
Model 3 = populist attitudes scale as the outcome; generations as the single predictor (using the four generational category recode we did in the tutorial), and using data for Greece. Interpret.
greece_data <- df %>%
mutate(
# Recoding gender
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
# Recoding education
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
),
# Handle NAs for education levels
edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
# Explicitly making 'No BA' the reference category
educ.ba = factor(educ.ba, levels = c("No BA", "BA or more")),
# Recoding age, setting 999 to NA
age = ifelse(agea == 999, NA, agea),
# Recoding cohort variable, setting years before 1930 and after 2000 to NA
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
# Recoding generational cohorts based on year of birth
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(cohort) # Keeping other values as character if they do not fit the ranges
),
# Converting cohort to a factor with labels
gen = factor(gen,
levels = c("1", "2", "3", "4"),
labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
)
models <- list(
"Model 1" = lm(auth ~ lrscale, data = greece_data),
"Model 2" = lm(auth ~ gndr, data = greece_data),
"Model 3" = lm(auth ~ gen, data = greece_data))
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 73.4 (0.2)*** | 74.4 (0.2)*** | 75.3 (0.5)*** |
| lrscale | 0.0 (0.0)+ | ||
| gndrMale | −2.6 (0.3)*** | ||
| genBaby Boomers | −1.1 (0.6)* | ||
| genGen X | −2.1 (0.6)*** | ||
| genMillennials | −3.6 (0.6)*** | ||
| Num.Obs. | 8432 | 8432 | 7631 |
| R2 | 0.000 | 0.007 | 0.006 |
| R2 Adj. | 0.000 | 0.007 | 0.006 |
| AIC | 69581.4 | 69523.1 | 62934.0 |
| BIC | 69602.5 | 69544.3 | 62968.7 |
| Log.Lik. | −34787.677 | −34758.566 | −31462.009 |
| RMSE | 14.98 | 14.93 | 14.94 |
modelplot(models, coef_omit = 'Interc')
right-leaning individuals and Baby Boomers are likely to have lower
populist attitudes, while gender (being male) also appears to be
significantly associated with lower populist attitudes in Greece.