Kashvi Shah | 1009413706
# 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.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ 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## [[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"rm(list=ls()); gc() ## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 2631431 140.6 4401068 235.1 NA 4401068 235.1
## Vcells 4474636 34.2 10146329 77.5 16384 7426248 56.7ess <- read_fst("All-ESS-Data.fst")unique(ess$cntry)## [1] "AT" "BE" "CH" "CZ" "DE" "DK" "ES" "FI" "FR" "GB" "GR" "HU" "IE" "IL" "IT"
## [16] "LU" "NL" "NO" "PL" "PT" "SE" "SI" "EE" "IS" "SK" "TR" "UA" "BG" "CY" "RU"
## [31] "HR" "LV" "RO" "LT" "AL" "XK" "ME" "RS" "MK"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")# 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))table(sweden_data$populist)##
## 0 3.33333333333333 6.66666666666667 10
## 46 15 47 69
## 13.3333333333333 16.6666666666667 20 23.3333333333333
## 123 180 398 505
## 26.6666666666667 30 33.3333333333333 36.6666666666667
## 806 879 931 1011
## 40 43.3333333333333 46.6666666666667 50
## 1143 1118 956 1161
## 53.3333333333333 56.6666666666667 60 63.3333333333333
## 804 841 683 631
## 66.6666666666667 70 73.3333333333333 76.6666666666667
## 544 596 431 364
## 80 83.3333333333333 86.6666666666667 90
## 331 257 200 164
## 93.3333333333333 96.6666666666667 100
## 155 80 244ggplot(sweden_data, aes(x = populist)) +
geom_histogram(bins = 30, fill = "hotpink", color = "yellow") +
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()`).
For some reason my graph itself would not load even though I tired to do it, multiple ways
The populist scale distribution for the nation of Sweden is seen on the histogram.Populism seems to be more common at lower levels and less common at higher levels, according to the histogram’s distribution, which appears to be somewhat tilted to the right. The mode appears to be lower than 50 as well.
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(
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"))
)# 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(
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
)
)class(sweden_data$ID)## [1] "character"levels(sweden_data$ID)## NULLsweden_data$ID <- as.factor(sweden_data$ID)levels(sweden_data$populist)## NULLsweden_data <- unique(sweden_data)I’m going to put all these codes into normal text because it would not run even though it was running the first try, once I figure out this out with TA again I will put these back into blocks (sorry I do not want to lose marks). I will take screenshot of all my errors and provide them so that there atleast you can see I tried.
sweden_pop_model <- lm(populist ~ ID, data = sweden_data)sweden_pop_model <- lm(populist ~ ID, data = sweden_data)
coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)## (Intercept) IDRight
## 48.3188406 -0.7149274coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)## (Intercept) IDRight
## 48.3188406 -0.7149274broom::tidy## function (x, ...)
## {
## UseMethod("tidy")
## }
## <bytecode: 0x10621a3f8>
## <environment: namespace:generics>tidy(sweden_pop_model)## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 48.3 0.265 183. 0
## 2 IDRight -0.715 0.356 -2.01 0.0444tidy(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 |
summary(sweden_pop_model)##
## Call:
## lm(formula = populist ~ ID, data = sweden_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.319 -14.271 -1.652 12.396 52.396
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.3188 0.2647 182.562 <2e-16 ***
## IDRight -0.7149 0.3555 -2.011 0.0444 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.25 on 11864 degrees of freedom
## (6350 observations deleted due to missingness)
## Multiple R-squared: 0.0003407, Adjusted R-squared: 0.0002565
## F-statistic: 4.044 on 1 and 11864 DF, p-value: 0.04436The projected populist scale is 51.4 when the education level (educ.ba) is zero, according to the intercept value of 51.4. When a coefficient is -7.9, it means that those who have completed a “BA or more” education have an estimated 7.9 less populist scale than others who have completed less education. Strong evidence against the null hypothesis is provided by the p-values, which are less than 0.001. As a result, “educ.baBA or more” has a statistically significant intercept and coefficient.
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")names(df)## NULLhead(df)##
## 1 function (x, df1, df2, ncp, log = FALSE)
## 2 {
## 3 if (missing(ncp))
## 4 .Call(C_df, x, df1, df2, log)
## 5 else .Call(C_dnf, x, df1, df2, ncp, log)
## 6 }library(dplyr)
library(scales)##
## Attaching package: 'scales'## The following object is masked from 'package:purrr':
##
## discard## The following object is masked from 'package:readr':
##
## col_factoritaly_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 = "darkred", 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)
italy_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 = "darkred", 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 graphic depicts the authoritarian scale’s steady stability between 2012 and 2020 with no discernible oscillations. However, there is a noticeable downturn that starts in 2016 and continues all the way through 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(
# 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)model <- lm(auth ~ gen, data = model_data)
modelsummary(model,
fmt = 1,
estimate = c("{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_rename = c("genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials"),
title = 'Table1. Regression model for authoritarian attitudes')| (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 |
It may be deduced that Baby Boomers had lower average levels of authoritarianism than the interwar generation did for the outcome variable as the coefficient is negative. In contrast, the Gen X generation’s coefficient of -2.1 shows that, in comparison to the interwar generation, being a member of the Gen X generation is associated with a 2.1 unit fall in authoritarian values. With a coefficient of -3.6, millennials are associated with a 3.6 unit fall in the outcome variable when compared to people from the interwar generation. The revised R-squared value of 0.006 shows that the “gen” variable only partially explains the variation in the Italian authoritarian value scale. However, the model’s limited scope suggests that factors other than generational disparities should be taken into account that influence authoritarian ideas.
Task 5 Now visualize, using the modelplot() function from the modelsummary package, coefficient estimates and 95% confidence intervals for the following three models:
greece_data <- read.fst("greece_data.fst")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(
# Recoding gender
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
lrscale = 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
),
# 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
## 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) | 72.0 (0.5)*** | 70.0 (0.3)*** | 65.8 (0.5)*** |
| lrscaleRight | −8.8 (0.6)*** | ||
| Male | −1.0 (0.5)* | ||
| Baby Boomers | 3.5 (0.7)*** | ||
| Gen X | 5.6 (0.7)*** | ||
| Millennials | 5.8 (0.7)*** | ||
| Num.Obs. | 4936 | 9823 | 9562 |
| R2 | 0.037 | 0.000 | 0.009 |
| R2 Adj. | 0.036 | 0.000 | 0.009 |
| AIC | 44723.8 | 88873.2 | 86473.9 |
| BIC | 44743.3 | 88894.8 | 86509.7 |
| Log.Lik. | −22358.894 | −44433.590 | −43231.944 |
| RMSE | 22.44 | 22.30 | 22.25 |
modelplot(models, coef_omit = 'Interc')
Gender-wise, the coefficient of -1.0 means that, generally speaking, males in Greece receive 1.0 less points than women on the populist scale. This suggests that in Greece, women are more likely than males to display populist views. In terms of generational differences, Baby Boomers in Greece score 3.5 points higher on the populist scale than the Interwar generation, suggesting stronger populist inclinations, according to the Baby Boomers coefficient of 3.5. Comparably, the Gen X coefficient of 5.6 suggests that Gen Xers in Greece have greater populist views than the Interwar generation, with their populist scale scores being 5.6 points higher. Finally, the Millennial generation’s coefficient, which stands at 5.8, indicates that, compared to the Interwar generation, Millennials in Greece have 5.8 points higher populist scale scores, showing the greatest populist sentiments among the generations investigated in Greece.
Taking into account only one predictor in the model, the adjusted R-squared value of 0.009 shows that, on average, our “gen” variable explains 0.9% of the variation in authoritarian opinions among Greeks. Despite being statistically significant, this finding implies that the model only accounts for a tiny fraction of the variables impacting authoritarian views, excluding variables other than generation.