Petsche_Lauren_Homework

Homework 4 (5%): due by next lecture on March 19th

packages <- c("tidyverse", "broom",  "infer", "fst", "modelsummary", "effects", "survey", "interactions") 

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.5
## ✔ 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] "broom"     "lubridate" "forcats"   "stringr"   "dplyr"     "purrr"    
##  [7] "readr"     "tidyr"     "tibble"    "ggplot2"   "tidyverse" "stats"    
## [13] "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[3]]
##  [1] "infer"     "broom"     "lubridate" "forcats"   "stringr"   "dplyr"    
##  [7] "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"   "tidyverse"
## [13] "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"  
## [19] "base"     
## 
## [[4]]
##  [1] "fst"       "infer"     "broom"     "lubridate" "forcats"   "stringr"  
##  [7] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [13] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [19] "methods"   "base"     
## 
## [[5]]
##  [1] "modelsummary" "fst"          "infer"        "broom"        "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"          "infer"       
##  [6] "broom"        "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"          "infer"        "broom"       
## [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"          "infer"       
## [11] "broom"        "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()`).

It can be noted that the populist scale for sweden is highest at just below 2000 at approximately just under 50.

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")

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 
##  396 2123 2852 3673 3702 1582

model_data <- df %>%
  mutate(
    trstplt = case_when(
      trstplt %in% 0:4 ~ "no trust",                     
      trstplt %in% 6:10 ~ "trust",                   
      trstplt %in% c(77, 88, 99) 
      ~ NA_character_,    
     TRUE ~ as.character(trstplt)  
    ),

    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$gen)

## < table of extent 0 >

model1 <- lm(auth ~ trstplt, data = model_data)
model2 <- lm(auth ~ educ.ba, data = model_data)

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

coefficients_sweden <- coef(sweden_pop_model)
print(coefficients_sweden)

## (Intercept)     trstplt 
##    88.87824    -8.59462

The expected average populus score of 88.8 suggest higher populist attitudes. -8.6 indicate the average individual identification as right have populist score -8.6 point lower than those that identify as left.

tidy(sweden_pop_model)

## # A tibble: 2 × 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)    88.9     0.137       648.       0
## 2 trstplt        -8.59    0.0270     -319.       0

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

term	estimate	std.error	statistic	p.value
(Intercept)	88.878	0.137	648.404	0
trstplt	-8.595	0.027	-318.685	0

A P value close to zero indicates that it is statistically significant.

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.

do not flip scale

dff <- read.fst("italy_data.fst")

dff <- dff %>%
  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 ))

dff$auth <- scales::rescale(dff$behave + 
                      dff$secure + 
                      dff$safety + 
                      dff$tradition + 
                      dff$rules, to=c(0,100), na.rm=TRUE)


dff <- dff %>% filter(!is.na(auth))
table(dff$secure)

## 
##    1    2    3    4    5    6 
##   54  139  626 1856 2876 2881

dff$year <- NA
replacements <- c(2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020)
for(i in 1:10){
  dff$year[dff$essround == i] <- replacements[i]
}
auth_avg <- dff %>%
  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)

There is very minimal fluctuation of the authoritarian scale average throughout the years for 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.

coiffenent rename function use. interpret coifficients

model_data2 <- dff %>%
  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_data2$gen)

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

model3 <- lm(auth ~ gen, data = model_data2)

modelsummary(
  list(model3),
  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

modelplot(model3, coef_omit = 'Interc')

modelsummary(model3, 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')

Table 2. Regression models 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

genBaby Boomers (-1,1) Compared to the Interwar generation
genGen X (-2.1): Generation X showed a decrease, with scores 2.1 points lower than the Interwar generation on the authoritarian values scale. The triple asterisks (***) indicates statistical significance at a very low p-value (under alpha threshold of 0.05)
genMillennials (-2.5): Millennials shows the most decrease in authoritarian values, at 2.5 points lower than the other generations. This is statistically significant, indicating that Millennials contrast the other generations in comparing their authoritarian values. The triple asterisks (***) denote statistical significance at a very low p-value (under our alpha threshold of 0.05)

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.

dx <- read.fst("greece_data.fst")

dx <- dx %>%
  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 ))

dx$auth <- scales::rescale(dff$behave + 
                      dx$secure + 
                      dx$safety + 
                      dx$tradition + 
                      dx$rules, to=c(0,100), na.rm=TRUE)

## Warning in dff$behave + dx$secure: longer object length is not a multiple of
## shorter object length

dx <- dx %>% filter(!is.na(auth))
table(dx$secure)

## 
##    1    2    3    4    5    6 
##   44  217  551 1347 3931 6240

model_data3 <- dx %>%
  mutate(
    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),

    
    yrbrn = 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)  
    ),
    lrscale = case_when(
      lrscale %in% 0:4 ~ "Left",                     
      lrscale %in% 6:10 ~ "Right",                   
      lrscale %in% c(77, 88, 99) ~ NA_character_,    
      TRUE ~ NA_character_                           
    ),

    yrbrn = factor(yrbrn,
                 levels = c("1", "2", "3", "4"),
                 labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
  )
table(model_data3$yrbrn)

## 
##     Interwar Baby Boomers        Gen X  Millennials 
##         2919         3425         3440         2285

model11 <- lm(auth ~ lrscale, data = model_data)
model22 <- lm(auth ~ gndr, data = model_data)
model33 <- lm(auth ~ yrbrn, data = model_data)

models <- list(
   "Model 11" = lm(auth ~ lrscale, data = model_data),
   "Model 22" = lm(auth ~ gndr, data = model_data),
   "Model 33" = lm(auth ~ yrbrn, data = model_data))
modelsummary(models, fmt = 1,
  estimate  = c("{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL)

	Model 11	Model 22	Model 33
(Intercept)	56.3 (0.2)***	54.8 (0.5)***	62.7 (2.0)***
lrscale	0.1 (0.0)***
gndr		1.4 (0.3)***
yrbrn			0.0 (0.0)**
Num.Obs.	14328	14328	14328
R2	0.003	0.002	0.001
R2 Adj.	0.003	0.001	0.001
AIC	122722.2	122737.5	122751.3
BIC	122744.9	122760.2	122774.0
Log.Lik.	−61358.091	−61365.729	−61372.661
RMSE	17.52	17.53	17.54

genBaby Boomers (-1,1) Compared to the Interwar generation
genGen X (-2.1): Generation X showed a decrease, with scores 2.1 points lower than the Interwar generation on the authoritarian values scale. The triple asterisks (***) indicates statistical significance at a very low p-value (under alpha threshold of 0.05)
genMillennials (-3.6): Millennials shows the most decrease in authoritarian values, at 3.6 points lower than the other generations. This is statistically significant, indicating that Millennials contrast the other generations in comparing their authoritarian values. this is statistical significance at a very low p-value
gndr is statistically significant

Petsche_Lauren_Homework_4

2024-03-12