Khurana_Ridhi_Homework_4
2024-03-18
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.
#Homework 4
Step 1: Set up enviorment
Insert all packages
# List of packages
packages <- c("tidyverse", "infer", "broom", "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.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## [[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] "broom" "infer" "lubridate" "forcats" "stringr" "dplyr"
## [7] "purrr" "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [13] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [19] "base"
##
## [[4]]
## [1] "fst" "broom" "infer" "lubridate" "forcats" "stringr"
## [7] "dplyr" "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [13] "tidyverse" "stats" "graphics" "grDevices" "utils" "datasets"
## [19] "methods" "base"
##
## [[5]]
## [1] "modelsummary" "fst" "broom" "infer" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[6]]
## [1] "effects" "carData" "modelsummary" "fst" "broom"
## [6] "infer" "lubridate" "forcats" "stringr" "dplyr"
## [11] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [16] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [21] "datasets" "methods" "base"
##
## [[7]]
## [1] "survey" "survival" "Matrix" "grid" "effects"
## [6] "carData" "modelsummary" "fst" "broom" "infer"
## [11] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [16] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"
##
## [[8]]
## [1] "interactions" "survey" "survival" "Matrix" "grid"
## [6] "effects" "carData" "modelsummary" "fst" "broom"
## [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?
##New Skill (scaling and rescaling)
Isolate Sweden data from ESS data set
sweden_data <- read.fst("sweden_data.fst")##Constructing our populist scale based on the Norris and Ingelhart (2019) conceptualization, which effectively uses trust items and reverses them
# 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))Distribution of values in Sweden data
ggplot(sweden_data, aes(x = populist)) +
geom_histogram(bins = 30, fill = "pink", color = "purple") +
theme_minimal() +
labs(title = "Distribution of Populist Scale for Hungary",
x = "Populist Scale",
y = "Count")## Warning: Removed 2503 rows containing non-finite values (`stat_bin()`).
What do I note: This graph shows a unimodal distribution. It is a type of distribution in a probability that illustrates one singular high point around the concentrated values. In this histogram, values occur mostly between 47.5 and 50 percent. This indictaes that there are more people who are interested in populism then those that are less interested in populism.
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.
df <- read.fst("sweden_data.fst")
model_data <- df %>%
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")),
ID = case_when(
lrscale %in% 0:4 ~ "Left",
lrscale %in% 6:10 ~ "Right",
lrscale %in% c(77, 88, 99) ~ NA_character_,
TRUE ~ NA_character_
),
)
model_data$trstplt <- ifelse(model_data$trstplt > 10, NA, model_data$trstplt)
model_data$trstprl <- ifelse(model_data$trstprl > 10, NA, model_data$trstprl)
model_data$trstprt <- ifelse(model_data$trstprt > 10, NA, model_data$trstprt)
model_data$trust <- scales::rescale(model_data$trstplt + model_data$trstprl + model_data$trstprt, na.rm = TRUE, to = c(0, 100))
model_data$populist <- scales::rescale(model_data$trust, na.rm = TRUE, to=c(100,0))
table(model_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 244table(model_data$educ.ba)##
## No BA BA or more
## 13477 4739Running a linear regression with populist attitudes
model1 <- lm(populist ~ educ.ba, data = model_data)Printing only the intercept and coefficients and then interpreting
coefficients_sweden <- coef(model1)
print(coefficients_sweden)## (Intercept) educ.baBA or more
## 51.358680 -7.868121Interpretation: The interception value is 51.358680, which is about 51.4. This number represents the average score on the populist scale for the “no BA” group.
-7.9 means that people who have a BA or more on average have populist scale scores that are 7.9 points LOWER than those who identify as “No BA”. This means that in the particular state, Sweden - people who do not have a BA have greater stronger populists attitudes than those who have a BA.
Do a tidy model table displaying the p-value
tidy(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 |
Interpenetration: p-value is 0: this tells the viewer that having no BA is statistically associated with stronger populist attitudes.
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.
Loading data and creating the authoritarian value scale
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))) %>%
mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))
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 297df$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)
table(auth_avg)## auth_avg
## year 72.6634225466342 72.780300115429 73.2129963898917 75.1807780320366
## 2012 0 0 0 1
## 2016 1 0 0 0
## 2018 0 1 0 0
## 2020 0 0 1 0Interpretation:
The is a minuscule chance - indicated by the signifinatly small decrease (75.2 for 2012, 72.6 for 2016) in the authoritarian scale average in Italy.
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.
LOAD DATA FOR ITALY
df <- read.fst("italy_data.fst")
model_data <- df %>%
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
## 1470 3179 2636 2031model_data <- model_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 ))
model_data$auth <- scales::rescale(model_data$behave +
model_data$secure +
model_data$safety +
model_data$tradition +
model_data$rules, to=c(0,100), na.rm=TRUE)
model_data <- model_data %>% filter(!is.na(auth))
table(model_data$secure)##
## 1 2 3 4 5 6
## 54 139 626 1856 2876 2881Running a linear regression model (using the model summary package)
model2 <- lm(auth ~ gen, data = model_data)
modelsummary(
list(model2),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_omit = "Intercept",
coef_rename = c("genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials")
)| (1) | |
|---|---|
| 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 |
tidy(model2)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 75.3 0.467 162. 0
## 2 genBaby Boomers -1.12 0.552 -2.04 4.14e- 2
## 3 genGen X -2.08 0.564 -3.69 2.29e- 4
## 4 genMillennials -3.59 0.584 -6.15 8.01e-10Coefficients Interpretation:
Baby boomers (-1.1, p<0.05), gen x (-2.1, p<0.001), millennials (-3.6, p<0.001), are all in comparisoon with the interwal generation which serves as the reference group in this case. Baby boomers are on average 1.1 lower on the authoritation values scale compared to the interwar generation, this means that on average baby boomers tend to have lower scores on these values comparatively to the reference group. Generation X show a significant decrease compared to interwar generation. Specifically, the scores of Generation X are reported to be 2.1 points lower, on average, on the authoritarian values scale compared to the Interwar generation. Lastly, millennials show the most decrease in authoritarian values compared to the Interwar generation. Specifically, the scores of Millennials are reported to be 3.6 points lower.
Adjusted R-Squared interpretation (R2 Adj. = 0.006): In this case the value is 0.006, this means about 0.6% of the variabilty in authortian values is actually explained by the generational categories in the regression model.
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.
Loading data for Greece
df <- read.fst("greece_data.fst")
model_data <- df %>%
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"))
)
model_data$trstplt <- ifelse(model_data$trstplt > 10, NA, model_data$trstplt)
model_data$trstprl <- ifelse(model_data$trstprl > 10, NA, model_data$trstprl)
model_data$trstprt <- ifelse(model_data$trstprt > 10, NA, model_data$trstprt)
model_data$trust <- scales::rescale(model_data$trstplt + model_data$trstprl + model_data$trstprt, na.rm = TRUE, to = c(0, 100))
model_data$populist <- scales::rescale(model_data$trust, na.rm = TRUE, to=c(100,0))
table(model_data$ID)##
## Left Right
## 2688 3522table(model_data$gndr)##
## Female Male
## 6929 5629table(model_data$gen)##
## Interwar Baby Boomers Gen X Millennials
## 2978 3484 3498 2330Creating the 3 models
model3 <- lm(populist ~ ID, data = model_data)
model4 <- lm(populist ~ gndr, data = model_data)
model5 <- lm(populist ~ gen, data = model_data)
models2 <- list(
"Model 1" = lm(populist ~ ID, data = model_data),
"Model 2" = lm(populist ~ gndr, data = model_data),
"Model 3" = lm(populist ~ gen, data = model_data))
modelsummary(models2, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
coef_rename = c("genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials"),
statistic = NULL,
coef_omit = "Intercept")| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| IDRight | −8.8 (0.6)*** | ||
| gndrMale | −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 |
tidy(model3)## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 72.0 0.478 151. 0
## 2 IDRight -8.80 0.643 -13.7 6.33e-42tidy(model4)## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 70.0 0.305 230. 0
## 2 gndrMale -0.967 0.452 -2.14 0.0324tidy(model5)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 65.8 0.508 130. 0
## 2 genBaby Boomers 3.55 0.661 5.37 7.97e- 8
## 3 genGen X 5.58 0.659 8.47 2.85e-17
## 4 genMillennials 5.79 0.706 8.21 2.51e-16Now visualize, using the modelplot() function from the modelsummary package, coefficient estimates and 95% confidence intervals
modelplot(models2, coef_omit = 'Interc', coef_rename = c("genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials"))
Interpretations:
Baby Boomers (3.6): The variance observed in Baby Boomers is moderately lower than that of Gen X and Millennials.
Millennials (5.8): Among the four generations examined, Millennials demonstrate the most significant rise in populist values, surpassing the interwar generation by 5.8 points. Consequently, they exhibit the strongest populist attitudes.
Male-Gender (-1.0): Greek males exhibit populist scores one point lower than females, indicating a slightly weaker tendency toward populist attitudes on average.
Gen X (5.6): Gen X also witnesses a notable surge in populist values, exceeding those of the interwar generation, albeit to a lesser extent than Millennials.
Political leaning (-8.8): In Greece, individuals with right-leaning political inclinations score 8.8 points lower on the populist values scale compared to their left-leaning counterparts. This indicates a weaker inclination towards populist attitudes among the politically right-leaning individuals.
confint(model3, level=0.95)## 2.5 % 97.5 %
## (Intercept) 71.0653 72.938927
## IDRight -10.0591 -7.539754confint(model4, level=0.95)## 2.5 % 97.5 %
## (Intercept) 69.416151 70.60995391
## gndrMale -1.852559 -0.08087043confint(model5, level=0.95)## 2.5 % 97.5 %
## (Intercept) 64.777052 66.768087
## genBaby Boomers 2.254559 4.845308
## genGen X 4.288926 6.872488
## genMillennials 4.409253 7.175453Interpretation:
The confidence interval spans from 61.0 to 65.4 on the populist scores scale, indicating a 95% confidence level that the true population parameter lies within this range, based on random samples.
For males, the confidence interval ranges from 67.5 (69.4-1.9) to 70.5 (70.6-0.1) on the populist scores scale, with the expectation that the true population parameter will be situated within this range across 95% of random samples.
Regarding Baby Boomers, the confidence interval extends from 67.0 (64.8+2.2) to 71.6 (66.8+4.8) on the populist scores scale, suggesting where the true population parameter is expected to lie with a 95% confidence level across random samples.
In the case of Gen X, the confidence interval spans from 69.1 (64.8+4.3) to 73.7 (66.8+6.9) on the populist scores scale, indicating where the true population parameter is anticipated to fall within 95% of random samples.
Similarly, for Millennials, the confidence interval ranges from 69.2 (64.8+4.4) to 74.0 (66.8+7.2) on the populist scores scale, with the expectation that the true population parameter will be within this range across 95% of random samples