Drzymala_Gabriela_Homework

# 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?

sweden_data <- read.fst("sweden_data.fst")

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()`).

Based on this histogram, I can note that it illustrates a normal distribution. Most of the data falls within the middle percentile, with a large peak at around 50% on the populist scale. Therefore, on average Swedish people idenitfy as populists to a moderate degree. ## 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.

table(sweden_data$edulvlb)

## 
##    0  113  213  313  321  323  413  423  520  610  620  710  720  800 5555 7777 
##  144  812 1214  944 1144 1144  325  893 1008  746  573  472  834  137   73    8 
## 8888 9999 
##   14   27

sweden_data <- sweden_data %>%
  mutate(
    ID = case_when(
      edulvla %in% 0:4 ~ "BA or more",                    
      edulvla %in% 5:10 ~ "BA or more",                   
      edulvla %in% c(55, 77, 88, 5555, 7777, 8888, 9999) ~ NA_character_ ,    
      TRUE ~  "No BA"                          
    )
  )


table(sweden_data$ID)

## 
## BA or more      No BA 
##       7679      10524

sweden_pop_model <- lm(populist ~ ID, data = sweden_data)

coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)

## (Intercept)     IDNo BA 
##    51.35224    -3.23235

Intercept (51.35):

This value represents the average score on the populist scale for the reference group in Sweden. The expected average populist scale score would be 51.35224 for BA or more, since it is the reference category. This score indicates a certain level of populism, where higher scores suggest stronger populist attitudes.

IDRight Coefficient (-3.23235):

The negative coefficient associated with the “No BA” group indicates that, on average, individuals who do not have a BA in Sweden have populist scale scores that are 3.23235 points lower than those who have a BA or more.This suggests that individuals with a BA or more exhibit stronger populist attitudes than their No BA counterparts.

tidy(sweden_pop_model)

## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    51.4      0.266    193.   0       
## 2 IDNo BA        -3.23     0.330     -9.80 1.29e-22

tidy(sweden_pop_model) |>
knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value
(Intercept)	51.352	0.266	192.854	0
IDNo BA	-3.232	0.330	-9.802	0

summary(sweden_pop_model)

## 
## Call:
## lm(formula = populist ~ ID, data = sweden_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -51.352 -14.787  -1.453  11.981  51.880 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  51.3522     0.2663 192.854   <2e-16 ***
## IDNo BA      -3.2323     0.3298  -9.802   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.68 on 15700 degrees of freedom
##   (2514 observations deleted due to missingness)
## Multiple R-squared:  0.006082,   Adjusted R-squared:  0.006019 
## F-statistic: 96.07 on 1 and 15700 DF,  p-value: < 2.2e-16

In both cases, p-value indicates it is statistically significant and that the null hypothesis is rejected.

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 <- read.fst("italy_data.fst")
## Creates an authoritarian values scale based on human modules items (can also do libertarian value scale)
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))) %>%
  # Apply the reverse coding
  mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))

# Now you can calculate 'schwartzauth' after the NA recoding
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 2881

table(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 297

ggplot(italy_data, 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")

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 = "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 yearly avergae for authoritarian values of Italians has remained consistent at around 75%. This indicates that on average, Italians have authoritarian tendencies given that they continue to score close to 75% on the authoritarian value scale.

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 gender
    gndr = case_when(
      gndr == 1 ~ "Male",
      gndr == 2 ~ "Female",
      gndr == 9 ~ NA_character_,
      TRUE ~ as.character(gndr)
    ),

 # Recoding education
    educ.ba = case_when(
      essround < 5 & edulvla == 5 ~ "BA or more",
      essround >= 5 & edulvlb > 600 ~ "BA or more",
      TRUE ~ "No BA"
    ),
    
    # Handle NAs for education levels
    edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla),
    edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
    
    # Explicitly making 'No BA' the reference category
    educ.ba = factor(educ.ba, 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(model_data$gen)

## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         1026         2580         2214         1811

model1 <- lm(auth ~ gen, data = model_data)
model2 <- lm(auth ~ gndr, data = model_data)
model3 <- lm(auth ~ educ.ba, 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

Coefficients Interpretation:

genBaby Boomers (-1.1, p<0.001): 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).
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.
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 strongly 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 explain approximately 0.6% of the variance in the authoritarian values scale for Italy, after accounting for the number of predictors. 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. ## 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 <- 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))
ggplot(greece_data, aes(x = populist)) +
  geom_histogram(bins = 30, fill = "blue", color = "black") +
  theme_minimal() +
  labs(title = "Distribution of Populist Scale for Greece",
       x = "Populist Scale",
       y = "Count")

## Warning: Removed 2735 rows containing non-finite values (`stat_bin()`).

table(greece_data$edulvlb)

## 
##    0  113  213  229  313  321  323  421  520  610  620  710  720  800 5555 7777 
##  251  837  576   30 1785  142  255  312   67  167  871   25  164   14    7    9 
## 8888 
##    2

greece_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_                          
    )
  )
table(greece_data$ID)

## 
##  Left Right 
##  2688  3522

greece_pop_model <- lm(populist ~ ID, data = greece_data)

coefficients_greece <- coef(greece_pop_model)
print(coefficients_greece)

## (Intercept)     IDRight 
##   72.002115   -8.799429

tidy(greece_pop_model)

## # 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-42

tidy(greece_pop_model) |>
knitr::kable(digits = 3)

term	estimate	std.error	statistic	p.value
(Intercept)	72.002	0.478	150.677	0
IDRight	-8.799	0.643	-13.695	0

summary(greece_pop_model)

## 
## Call:
## lm(formula = populist ~ ID, data = greece_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -72.002 -16.536   1.331  17.998  36.797 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  72.0021     0.4779   150.7   <2e-16 ***
## IDRight      -8.7994     0.6425   -13.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.44 on 4934 degrees of freedom
##   (7622 observations deleted due to missingness)
## Multiple R-squared:  0.03662,    Adjusted R-squared:  0.03642 
## F-statistic: 187.5 on 1 and 4934 DF,  p-value: < 2.2e-16

model_data <- greece_data %>%
  mutate(
    # Recoding gender
    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)  # 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         2330

model1 <- lm(populist ~ lrscale, data = model_data)
model2 <- lm(populist ~ gndr, data = model_data)
model3 <- lm(populist ~ gen, data = model_data)

modelsummary(
  list(model1),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_omit = "Intercept")

	(1)
lrscale	0.1 (0.0)***
Num.Obs.	9823
R2	0.021
R2 Adj.	0.021
AIC	88667.6
BIC	88689.2
Log.Lik.	−44330.802
RMSE	22.07

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)

	Model 1	Model 2	Model 3
(Intercept)	67.5 (0.3)***	70.0 (0.3)***	65.8 (0.5)***
lrscale	0.1 (0.0)***
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.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')

modelsummary(models, fmt = 1,
  estimate  = c("{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_rename = c("genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials", "gndrMale" = "Male", "educ.ba or more" = "left"),
  title = 'Table 2. Regression models for populist attitudes')

Table 2. Regression models for populist attitudes
	Model 1	Model 2	Model 3
(Intercept)	67.5 (0.3)***	70.0 (0.3)***	65.8 (0.5)***
lrscale	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

Male, relative to female illustrates that on average, men have populist scale scores that are 1.0 points lower than feamles.

The generations model shows no statistical significance.

#End

Drzymala_Gabriela_Homework_4

2024-03-18

Task 1

Task 3

Task 4