Homework 4
Mingdi Zhou
2024-03-17
Setting up environment
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "interactions") # 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: 程辑包'infer'是用R版本4.3.3 来建造的## Warning: 程辑包'effects'是用R版本4.3.3 来建造的## 载入需要的程辑包：carData## Warning: 程辑包'carData'是用R版本4.3.3 来建造的## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Warning: 程辑包'survey'是用R版本4.3.3 来建造的## 载入需要的程辑包：grid
## 载入需要的程辑包：Matrix
##
## 载入程辑包：'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## 载入需要的程辑包：survival
##
## 载入程辑包：'survey'
##
## The following object is masked from 'package:graphics':
##
## dotchart## Warning: 程辑包'interactions'是用R版本4.3.3 来建造的## [[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"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")# Populist scale (Norris/Inglehart use low levels of trust in parliaments, politicians and parties as an indicator of populism)
# Note: Schafer does it 2 ways, we will focus now on the most straightforward
# Setting values greater than 10 to NA
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()`).
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
)
)
table(sweden_data$ID)##
## Left Right
## 6188 7406sweden_pop_model <- lm(populist ~ ID, data = sweden_data)coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)## (Intercept) IDRight
## 48.3188406 -0.7149274Answer:
Intercept (48.3188046):
This value represents the average score on the populist scale for the reference group in Sweden, which is the “Left” group. In simpler terms, the expected average populist scale score would be 48.32 for the left, since it is the reference category. This score indicates a certain level of populism, where higher scores suggest stronger populist attitudes.
IDRight Coefficient (-0.7149274):
The negative coefficient associated with the “Right” group indicates that, on average, individuals who identify as ideologically right in Sweden have populist scale scores that are 0.71 points lower than those who identify as left (note: based on our recode).
Based on what we learnt for interpreting coefficients, we might say: “Comparing right-leaning individuals to left-leaning ones in Sweden, we expect to see the average populist scale score to be 0.71 points lower for the right-leaning group.” This suggests that left-leaning individuals exhibit slightly stronger populist attitudes than their right-leaning counterparts.
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.
# Load libraries
library(dplyr)
library(broom)## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtoolssweden_data <- sweden_data %>%
mutate(
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
)
)%>%
# Ensure all NAs are handled uniformly
mutate(
edulvla = ifelse(edulvla %in% c(77, 88), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba) # Force educ.ba to be treated as character
)
# create linear regression model
sweden_pop_model <- lm(populist ~ educ.ba, data = sweden_data)
# Print the coefficients
coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)## (Intercept) educ.baNo BA
## 43.490559 7.868121# Create a tidy model table with p-value
tidy_model <- tidy(sweden_pop_model) %>%
mutate(p.value = format.pval(p.value, digits = 3))
# Print the tidy model table
print(tidy_model)## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 (Intercept) 43.5 0.299 145. <2e-16
## 2 educ.baNo BA 7.87 0.350 22.5 <2e-16Interpretation:
The positive coefficient for “No BA” suggests that lower educational attainment (not having a BA degree) is associated with higher scores on the populist scale. This could be interpreted as a higher tendency towards populist attitudes among those with less formal education, as measured in your study.
The small p-values for both the intercept and the educ.ba No BA term strongly suggest that these findings are not due to random chance, and that there is a statistically significant relationship between educational attainment and populist attitudes in the data.
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")
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 <- df %>% filter(!is.na(auth))
table(df$secure)##
## 1 2 3 4 5 6
## 54 139 626 1856 2876 2881table(df$auth)##
## 0 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84
## 7 4 4 9 13 20 27 57 113 155 163 323 367 596 647 789 823 897 992 733
## 88 92 96 100
## 632 465 299 297ggplot(df, aes(x = auth)) +
geom_histogram(bins = 30, fill = "blue", color = "black") +
theme_minimal() +
labs(title = "Distribution of Authoritarian Scale for Italy",
x = "Authoritarian Scale",
y = "Count")
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)
Interpretation:
The output displays the average scores on an authoritarian scale over several survey years.
The blue dashed line in the graph indicates the average score for each year. The graph shows a relatively stable trend over the years, with the average score staying around the same level, indicating that the authoritarian values measured by the scale have not changed significantly among the survey participants in Italy over the observed period.
There’s no significant upward or downward trend observable, suggesting that the Italian population’s authoritarian values, as captured by this particular scale, have remained consistent over the time span from 2002 to 2020.
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.
model_data <- df %>%
mutate(
# 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"))
)
table(model_data$gen)##
## Interwar Baby Boomers Gen X Millennials
## 1026 2580 2214 1811model1 <- lm(auth ~ gen, data = model_data)modelsummary(
list(model1),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_omit = "Intercept")| (1) | |
|---|---|
| genBaby Boomers | -1.1 (0.6)* |
| genGen X | -2.1 (0.6)*** |
| genMillennials | -3.6 (0.6)*** |
| Num.Obs. | 7631 |
| R2 | 0.006 |
| R2 Adj. | 0.006 |
| AIC | 62934.0 |
| BIC | 62968.7 |
| Log.Lik. | -31462.009 |
| RMSE | 14.94 |
Interpretation:
Coefficients Interpretation:
genBaby Boomers (-1.1, p<0.05): Compared to the Interwar generation (reference category that is omitted), Baby Boomers exhibit scores that are, on average, 1.1 points lower on the authoritarian values scale (0-100). The single asterisk (*) denotes statistical significance at a p-value between 0.01 and 0.05 (under our alpha threshold of 0.05).
genGen X (-2.1, p<0.001): Generation X shows an even more pronounced decrease, with scores 2.1 points lower than the Interwar generation on the authoritarian values scale. The triple asterisks (***) denote statistical significance at a very low p-value (under our alpha threshold of 0.05).
genMillennials (-3.6, p<0.001): 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 significantly 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 (with our four coded generations) explains approximately 0.6% of the variance in the authoritarian values scale for Italy, after accounting for the number of predictors (three in this case). While the generation coefficients are statistically significant, this suggests that the model explains a very 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 <- read.fst("greece_data.fst")greece_data <- greece_data model_data <- df %>%
mutate(
# Recoding gender
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
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
),
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
# 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),
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"))
)model1 <- lm(auth ~ lrscale, data = model_data)
model2 <- lm(auth ~ gndr, data = model_data)
model3 <- lm(auth ~ gen, data = model_data)modelsummary(
list(model1),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_omit = "Intercept")| (1) | |
|---|---|
| lrscale | 0.0 (0.0)+ |
| Num.Obs. | 8432 |
| R2 | 0.000 |
| R2 Adj. | 0.000 |
| AIC | 69581.4 |
| BIC | 69602.5 |
| Log.Lik. | -34787.677 |
| RMSE | 14.98 |
models <- list(
model1 <- lm(auth ~ ID, data = model_data),
model2 <- lm(auth ~ gndr, data = model_data),
model3 <- lm(auth ~ gen, data = model_data))
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL)| (1) | (2) | (3) | |
|---|---|---|---|
| (Intercept) | 71.7 (0.3)*** | 74.4 (0.2)*** | 75.3 (0.5)*** |
| IDRight | 2.9 (0.4)*** | ||
| gndrMale | -2.6 (0.3)*** | ||
| genBaby Boomers | -1.1 (0.6)* | ||
| genGen X | -2.1 (0.6)*** | ||
| genMillennials | -3.6 (0.6)*** | ||
| Num.Obs. | 4604 | 8432 | 7631 |
| R2 | 0.010 | 0.007 | 0.006 |
| R2 Adj. | 0.009 | 0.007 | 0.006 |
| AIC | 37859.6 | 69523.1 | 62934.0 |
| BIC | 37878.9 | 69544.3 | 62968.7 |
| Log.Lik. | -18926.783 | -34758.566 | -31462.009 |
| RMSE | 14.76 | 14.93 | 14.94 |
modelplot(models, coef_omit = 'Interc')
Answer:
Model 1: This model uses ‘lrscale’ as the predictor. The coefficient for IDRight is positive (2.9) and statistically significant, meaning right-leaning individuals score approximately 2.9 points higher on the populist attitudes scale than left-leaning individuals. The R-squared value is 0.010, suggesting that political orientation only explains about 1% of the variance in populist attitudes in Greece, indicating other factors also play a significant role.
Model 2: Here, ‘gndr’ is used as the predictor. The negative coefficient for gndrMale (-2.6) indicates that males have lower scores on the populist attitudes scale compared to females, and this result is statistically significant. The R-squared value is 0.007, which is quite low, signifying that gender alone does not explain much of the variation in populist attitudes.
Model 3: The ‘gen’ variable is used as the predictor, with coefficients for each generation compared to the omitted category (presumably the oldest generation not shown). The negative coefficients for genBaby Boomers, genGen X, and genMillennials suggest that these generations score lower on the populist attitudes scale than the reference generation. The R-squared value is 0.006, meaning only 0.6% of the variance in populist attitudes is accounted for by generational differences.
And the plot visualizes the coefficient estimates from the models and their 95% confidence intervals.
For Model 1 (IDRight), the plot shows a positive coefficient with a 95% confidence interval that does not include zero, confirming its statistical significance.
For Model 2 (gndrMale), the negative coefficient and its 95% confidence interval also do not include zero, supporting the statistical significance of males having lower populist attitudes.
For Model 3, all generational categories have negative coefficients, with 95% confidence intervals that do not include zero, indicating that these are statistically significant findings.
The AIC and BIC values across all models are quite high, suggesting that the models may not be the best fit for the data.
In conclusion, while each model has identified statistically significant predictors of populist attitudes in Greece, the low R-squared values in all models indicate that the predictors—political orientation (lrscale), gender (gndr), and generation (gen)—each explain a very small portion of the variability in populist attitudes. This suggests that more complex models including additional variables may be needed to more fully understand the determinants of populist attitudes.