Markiewicz_Alexandra_Homework_4
2024-03-16
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "interactions", "broom")
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_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
## Loading required package: carData
##
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
##
## 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
##
##
## 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"sweden_data <- read.fst("sweden_data.fst")
italy_data <- read.fst("italy_data.fst")
greece_data <- read.fst("greece_data.fst")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$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)
sweden_data$trust <- scales::rescale(sweden_data$trstplt + sweden_data$trstprl + sweden_data$trstprt, na.rm = TRUE, to = c(0, 100))
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 histogram shows the distribution of the populist scale for Sweden. The horizontal axis depicts the populist scale, re-scaled to a 0 to 100 range, while the vertical axis indicates the count of observations within each range. It seems to have a uni-modal distribution, with a peak near the middle of the horizontal axis, suggesting that most observations have a populist score within a particular range near 50. The histogram is right-skewed, demonstrating fewer observations with higher populist scores than lower ones.
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.
model_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"
),
edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = factor(educ.ba, levels = c("No BA", "BA or more"))
)
table(model_data$educ.ba)##
## No BA BA or more
## 13477 4739model1 <- lm(populist ~ educ.ba, data = model_data)models <- list(
"Model 1" = lm(populist ~ educ.ba, data = model_data))
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_rename = c("educ.baBA or more" = "BA or More"),
title = 'Table 1. Regression models for populist attitudes')| Model 1 | |
|---|---|
| (Intercept) | 51.4 (0.2)*** |
| BA or More | −7.9 (0.3)*** |
| Num.Obs. | 15713 |
| R2 | 0.031 |
| R2 Adj. | 0.031 |
| AIC | 137833.2 |
| BIC | 137856.2 |
| Log.Lik. | −68913.611 |
| RMSE | 19.43 |
modelplot(models, coef_omit = 'Interc')
tidy(model1)## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 51.4 0.181 283. 0
## 2 educ.baBA or more -7.87 0.350 -22.5 2.58e-110tidy(model1) |>
knitr::kable(digits = 3)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 51.359 | 0.181 | 283.363 | 0 |
| educ.baBA or more | -7.868 | 0.350 | -22.494 | 0 |
summary(model1)##
## Call:
## lm(formula = populist ~ educ.ba, data = model_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -51.359 -14.692 -1.359 11.975 56.509
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 51.3587 0.1812 283.36 <2e-16 ***
## educ.baBA or more -7.8681 0.3498 -22.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.43 on 15711 degrees of freedom
## (2503 observations deleted due to missingness)
## Multiple R-squared: 0.0312, Adjusted R-squared: 0.03114
## F-statistic: 506 on 1 and 15711 DF, p-value: < 2.2e-16The intercept value of 51.4 represents the average score on the populist scale for the reference group in Sweden, the “BA or More” group. When this value is higher, it implies that individuals within the group tend to show stronger populist attitudes.
The educ.baBA or More coefficient value of -7.9 indicates that, on average, individuals who have a BA or more in Sweden have populist scale scores that are 7.9 points lower than those who have no BA, suggesting that individuals with no BA exhibit stronger populist attitudes than individuals with a BA or more.
The p-value of <2e-16, for both the intercept and the BA or More coefficient, indicates statistical significance and strong evidence against the null hypothesis, suggesting having a BA or more is associated with a decrease in populist attitudes in Sweden.
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.
italy_data <- read.fst("italy_data.fst")
italy_data <- italy_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))) %>%
mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))
italy_data$auth <- scales::rescale(italy_data$behave +
italy_data$secure +
italy_data$safety +
italy_data$tradition +
italy_data$rules, to=c(0,100), na.rm=TRUE)
italy_data <- italy_data %>% filter(!is.na(auth))
table(italy_data$secure)##
## 1 2 3 4 5 6
## 54 139 626 1856 2876 2881table(italy_data$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 297italy_data$year <- NA
replacements <- c(2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020)
for(i in 1:10){
italy_data$year[italy_data$essround == i] <- replacements[i]
}
auth_avg <- italy_data %>%
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)
The graph shows the distribution of the average authoritarian scale in Italy over the years. The horizontal axis depicts the survey years from 2012 to 2020, while the vertical axis indicates the authoritarian scale from 0 to 100. Each point on the graph depicts the average authoritarian scale score for a certain survey year, with the blue dashed line linking them, demonstrating the general trend of the average authoritarian scale score over the years. It displays a peak in authoritarian values in 2012, shown by the highest point on the graph, and an incremental decrease until 2016. However, the trend is shown to reverse slightly in 2016, with a slight increase in authoritarian values until 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 <- italy_data %>%
mutate(
age = ifelse(agea == 999, NA, agea),
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)
),
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)models <- list(
"Model 1" = lm(auth ~ gen, data = model_data))
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_rename = c("genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials"),
title = 'Table 2. Regression models for authoritarian attitudes')| Model 1 | |
|---|---|
| (Intercept) | 75.3 (0.5)*** |
| Baby Boomers | −1.1 (0.6)* |
| Gen X | −2.1 (0.6)*** |
| Millennials | −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 |
modelplot(models, coef_omit = 'Interc')
The genBaby Boomers coefficient value of -1.1 indicates that, on average, Baby Boomers in Italy have authoritarian scale scores that are 1.1 points lower than the Interwar generation, suggesting that the Interwar generation exhibits stronger authoritarian attitudes than Baby Boomers. The ’*’ indicates statistical significance with a p-value below the alpha threshold, 0.05.
The genGen X coefficient value of -2.1 indicates that, on average, Gen X in Italy has authoritarian scale scores that are 2.1 points lower than the Interwar generation, suggesting that the Interwar generation exhibit stronger authoritarian attitudes than Gen X. The ’***’ indicates statistical significance with a p-value significantly below the alpha threshold, 0.05.
The genMillennials coefficient value of -3.6 is the lowest of the three generations, indicating that, on average, Millennials in Italy have authoritarian scale scores that are 3.6 points lower than the Interwar generation, suggesting that the Interwar generation exhibits stronger authoritarian attitudes than Millennials. The ’***’ indicates statistical significance with a p-value significantly below the alpha threshold, 0.05.
The adjusted R-squared value of 0.006 means that the ‘gen’ variable explains approximately 0.6% of the variance in the authoritarian values scale for Italy. Even though it is statistically significant, it suggests that generational differences alone do not account for much variation, implying other important factors influence authoritarian values.
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$trstplt <- ifelse(greece_data$trstplt > 10, NA, greece_data$trstplt)
greece_data$trstprl <- ifelse(greece_data$trstprl > 10, NA, greece_data$trstprl)
greece_data$trstprt <- ifelse(greece_data$trstprt > 10, NA, greece_data$trstprt)
greece_data$trust <- scales::rescale(greece_data$trstplt + greece_data$trstprl + greece_data$trstprt, na.rm = TRUE, to = c(0, 100))
greece_data$populist <- scales::rescale(greece_data$trust, na.rm = TRUE, to=c(100,0))model_data <- greece_data %>%
mutate(
ID = case_when(
lrscale %in% 0:4 ~ "Left",
lrscale %in% 6:10 ~ "Right",
lrscale %in% c(77, 88, 99) ~ NA_character_,
TRUE ~ NA_character_
),
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
age = ifelse(agea == 999, NA, agea),
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)
),
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
## 2978 3484 3498 2330model1 <- lm(populist ~ lrscale, data = model_data)
model2 <- lm(populist ~ gndr, data = model_data)
model3 <- lm(populist ~ gen, data = model_data)models <- list(
"Model 1" = lm(populist ~ lrscale, data = model_data),
"Model 2" = lm(populist ~ gndr, data = model_data),
"Model 3" = lm(populist ~ gen, data = model_data))
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_rename = c("lrscale" = "Right", "gndrMale" = "Male", "genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials"),
title = 'Table 3. Regression models for populist attitudes')| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 67.5 (0.3)*** | 70.0 (0.3)*** | 65.8 (0.5)*** |
| Right | 0.1 (0.0)*** | ||
| Male | −1.0 (0.5)* | ||
| Baby Boomers | 3.5 (0.7)*** | ||
| Gen X | 5.6 (0.7)*** | ||
| Millennials | 5.8 (0.7)*** | ||
| Num.Obs. | 9823 | 9823 | 9562 |
| R2 | 0.021 | 0.000 | 0.009 |
| R2 Adj. | 0.021 | 0.000 | 0.009 |
| AIC | 88667.6 | 88873.2 | 86473.9 |
| BIC | 88689.2 | 88894.8 | 86509.7 |
| Log.Lik. | −44330.802 | −44433.590 | −43231.944 |
| RMSE | 22.07 | 22.30 | 22.25 |
modelplot(models, coef_omit = 'Interc')
The lrscale coefficient value of 0.1 indicates that, on average, individuals who identify as ideologically right in Greece have populist scale scores that are 0.1 points higher than those who identify as left, suggesting that individuals who identify as right exhibit stronger populist attitudes than individuals who identify as left.
The gndrMale coefficient value of -1.0 indicates that, on average, men in Greece have populist scale scores that are 1.0 points lower than women, suggesting that women exhibit stronger populist attitudes than men.
The genBaby Boomers coefficient value of 3.5 indicates that, on average, Baby Boomers in Greece have populist scale scores that are 3.5 points higher than the Interwar generation, suggesting that Baby Boomers exhibit stronger populist attitudes than the Interwar generation.
The genGen X coefficient value of 5.6 indicates that, on average, Gen X in Greece has populist scale scores that are 5.6 points higher than the Interwar generation, suggesting that Gen X exhibits stronger populist attitudes than the Interwar generation.
The genMillennials coefficient value of 5.8 is the highest of the three generations, indicating that, on average, Millennials in Greece have populist scale scores that are 5.8 points higher than the Interwar generation, suggesting that Millennials exhibit stronger populist attitudes than the Interwar generation.