Hakim_Rashad_Homework_4
2024-03-16
# 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.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] "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"ess <- read_fst("All-ESS-Data.fst")sweden_data <- read.fst("sweden_data.fst")Homework 4 (5%): due by next lecture on March 19th
Instructions: Start a new R markdown for the homework and call it “Yourlastname_Firstname_Homework_4”.
Copy everything below from Task 1 to Task 5. Keep the task prompt and questions, and provide your code and answer underneath.
Remember: you need all the steps for your code to work, including loading your data – otherwise it will not knit.
To generate a new code box, click on the +C sign above. Underneath your code, provide your answer to the task question.
When you are done, click on “Knit” above, then “Knit to Html”. Wait for everything to compile. If you get an error like “Execution halted”, it means there are issues with your code you must fix. When all issues are fixed, it will prompt a new window. Then click on “Publish” in the top right, and then Rpubs (the first option) and follow the instructions to create your Rpubs account and get your Rpubs link for your document (i.e., html link as I provide for the tutorial).
Note: Make sure to provide both your markdown file and R pubs link. If you do not submit both, you will be penalized 2 pts. out of the 5 pts. total.
Task 1
##ANSWER: The Populist scale on the left tend to be around the 40-50 range which indicates that there is around half or even amount of populist attitudes.
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()`).
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.
##ANSWER: Sweden data : The intercept at the number 49.7 means that the populist attitudes are almost even on the left side of the populist scale. On the other hand, the ID “No-BA” of 5.37 means that those who have No BA have populist scores 5.37 highier than those who are “BA or More”. P Value of 0 means that it is statistically significant and that “No BA” is has a slight increase in populist attitudes sompared to those with “BA or more”.
sweden_data <- sweden_data %>%
mutate(
ID = case_when(
edulvla %in% 0:2 ~ "No BA",
edulvla %in% 3:5 ~ "BA or more",
lrscale %in% c(55,77, 88, 99) ~ NA_character_,
TRUE ~ NA_character_
)
)
table(sweden_data$ID)##
## BA or more No BA
## 5211 2468sweden_pop_model <- lm(populist ~ ID, data = sweden_data)coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)## (Intercept) IDNo BA
## 49.727249 5.375657broom::tidy(sweden_pop_model)## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 49.7 0.313 159. 0
## 2 IDNo BA 5.38 0.569 9.45 4.80e-21broom::tidy(sweden_pop_model) |>
knitr::kable(digits = 3)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 49.727 | 0.313 | 159.032 | 0 |
| IDNo BA | 5.376 | 0.569 | 9.452 | 0 |
summary(sweden_pop_model)##
## Call:
## lm(formula = populist ~ ID, data = sweden_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55.10 -13.06 -1.77 13.61 50.27
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.7272 0.3127 159.032 <2e-16 ***
## IDNo BA 5.3757 0.5687 9.452 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.31 on 5463 degrees of freedom
## (12751 observations deleted due to missingness)
## Multiple R-squared: 0.01609, Adjusted R-squared: 0.01591
## F-statistic: 89.34 on 1 and 5463 DF, p-value: < 2.2e-16Task 3
Now using the Italian dataset and the authoritarian values scale (as we did in the tutorial), graph the average by survey year. Interpret.
##ANSWER: Based off the authoritarian scale per year for Italy plotted below the 5 categories behave, secure, safety, tradition and rules for Italy under the authoritarian scale stays the same at around the 70-75 range. Which means that majority of people living in Italy are content with how Italy is run based of the 5 categories mentioned above.
italy_data <- read.fst("italy_data.fst")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)
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.
##ANSWER: Baby Boomers -1.1 and P Value of 0.6 means that the generation of Baby Boomers have a scores thats 1.1 points lower on the authoritarian values scale based on the 0-100 scale (0-100). The P value of 0.6 and the one asterik (*) means that statistical significance is little highier than the threshold which is 0.5.
##ANSWER: Gen X -2.1 and P Value of 0.006 means that Gen X have scores thats 2.1 points lower on our authoritarian values scale based from (0-100). The P value of 0.6 and the three asterik (***) means that statistical significance is also little highier than the threshold which is 0.5.
##ANSWER: Millennials -3.6 and P Value of 0.6 means that the Millennials have 3.6 points lower on our authoritarian values scale which is based on a scale of (0-100). The P value of 0.6 and the three asterik (***) means that statistical significance is highier than the other three generations.
##ANSWER Adjusted R-Squared (Adj. R² = 0.006): For the adjusted r squared value which is 0.006, the meaning behind that is that 0.6 percent of the variance is explained (approximately) for the authoritarian values for Italy.
The adjusted R-squared value of 0.031 means that the ‘gen’ variable (with our four coded generations) explain approximately 3.1% of the variance in the authoritarian values scale for Portugal, after accounting for the number of predictors (just one in this case). 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.
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 2881df <- df %>%
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
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")),
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(df$gen)##
## Interwar Baby Boomers Gen X Millennials
## 1026 2580 2214 1811model1 <- lm(auth ~ gen, data = df)
model2 <- lm(auth ~ gndr, data = df)
model3 <- lm(auth ~ educ.ba, data = df)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 |
modelsummary(model1, fmt = 1,
estimate = c("{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 = 'Table 2. Regression models for authoritarian attitudes for Italy')| (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 |
Task 5
Now visualize, using the modelplot() function from the modelsummary package, coefficient estimates and 95% confidence intervals for the following three models:
##ANSWER: Model 1: The intercept at the number 70.3 means that our populist attitudes are highier on the left side of the populist scale. The ID “BA or more” of -5.6 means that those who have No BA have populist scores -5.6 lower than people who are “No BA”
##ANSWER: Model 2: Gender based on males have -1.0 lower populist attitudes in comparison to females on the Populist scale.
##ANSWER: Model 3:
##ANSWER: Baby Boomers 3.5 and P Value of 0.7 means that the generation of Baby Boomers have a scores thats 3.5 points higher on the populist attitudes scale based on the 0-100 scale (0-100). The P value of 0.7 and the three asteriks (***) means that statistical significance is highier than the threshold which is 0.5.
##ANSWER: Gen X 5.6 and P Value of 0.7 means that Gen X have scores thats 5.6 points higher on populist attitudes scale based from (0-100). The P value of 0.7 and the three asterik (***) means that statistical significance is highier than the threshold which is 0.7.
##ANSWER: Millennials 5.8 and P Value of 0.7 means that the Millennials have 5.8 points higher on our populist attitudes scale which is based on a scale of (0-100). The P value of 0.7 and the three asterik (***) means that statistical significance is high.
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)
# Creating and rescaling the trust variable
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))greece_data <- greece_data %>%
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
),
ltrscale = case_when(
essround < 5 & lrscale == 5 ~ "BA or more",
essround >= 5 & lrscale > 5 ~ "No BA",
TRUE ~ "No BA"
),
# Handle NAs for education levels
lrscale = ifelse(lrscale %in% c(77, 88, 99), NA_integer_, lrscale),
# Explicitly making 'No BA' the reference category
ltrscale = factor(ltrscale, 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"))
)
table(greece_data$gen)##
## Interwar Baby Boomers Gen X Millennials
## 2978 3484 3498 2330modell1 <- lm(populist ~ ltrscale, data = greece_data)
modell2 <- lm(populist ~ gndr, data = greece_data)
modell3 <- lm(populist ~ gen, data = greece_data)modelss <- list(
"Modell 1" = lm(populist ~ ltrscale, data = greece_data),
"Modell 2" = lm(populist ~ gndr, data = greece_data),
"Modell 3" = lm(populist ~ gen, data = greece_data))
modelsummary(modelss, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL)| Modell 1 | Modell 2 | Modell 3 | |
|---|---|---|---|
| (Intercept) | 70.3 (0.2)*** | 70.0 (0.3)*** | 65.8 (0.5)*** |
| ltrscaleBA or more | −5.6 (0.7)*** | ||
| gndrMale | −1.0 (0.5)* | ||
| genBaby Boomers | 3.5 (0.7)*** | ||
| genGen X | 5.6 (0.7)*** | ||
| genMillennials | 5.8 (0.7)*** | ||
| Num.Obs. | 9823 | 9823 | 9562 |
| R2 | 0.007 | 0.000 | 0.009 |
| R2 Adj. | 0.007 | 0.000 | 0.009 |
| AIC | 88803.9 | 88873.2 | 86473.9 |
| BIC | 88825.4 | 88894.8 | 86509.7 |
| Log.Lik. | −44398.934 | −44433.590 | −43231.944 |
| RMSE | 22.22 | 22.30 | 22.25 |
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.