Shah_Kashvi_Homework5

#Setting Up the Environment

# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "interactions","broom") # 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
## 
## 
## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools

## [[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"        
## 
## [[8]]
##  [1] "broom"        "interactions" "survey"       "survival"     "Matrix"      
##  [6] "grid"         "effects"      "carData"      "modelsummary" "fst"         
## [11] "infer"        "lubridate"    "forcats"      "stringr"      "dplyr"       
## [16] "purrr"        "readr"        "tidyr"        "tibble"       "ggplot2"     
## [21] "tidyverse"    "stats"        "graphics"     "grDevices"    "utils"       
## [26] "datasets"     "methods"      "base"

ess <- 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 Use data for Germany and the variable coding for model3 and model4 from the tutorial (i.e., both for the cleaning & recode, as well as the linear model formulas). Do a modelsummary table displaying both model outputs. Interpret the coefficients for the MLR without the interaction (i.e., the first displayed model), and interpret model fit metrics for both models.

#Here we read the germany data fst
germany_data <- read_fst("germany_data.fst")

# Renaming data set to 'df' , so whenever 'df' is used it actually means 'germany_data'
df <- germany_data

# Preparing weighted data set
df$weight <- df$dweight * df$pweight
survey_design <- svydesign(ids = ~1, data = df, weights = ~weight)

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))) %>%
  # Apply the reverse coding
  mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))

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

df <- df %>%
  mutate(
    polID = case_when(
      lrscale %in% 0:3 ~ "Left",                    
      lrscale %in% 7:10 ~ "Right",                     
      lrscale %in% 4:6 ~ "Moderate",                  
      lrscale %in% c(77, 88, 99) ~ NA_character_      
    ),
   religious = case_when(
      rlgdgr %in% c(77, 88, 99) ~ NA_real_,
      TRUE ~ rlgdgr
    )
  )

# Recoding the varaible
df <- df %>%
  mutate(
    cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
    # Recoding generational cohorts based on the year of birth (yrbrn).
    # The year of birth is categorized into different generational cohorts.
    # Interwar (1900-1945), Baby Boomers (1946-1964), Gen X (1965-1979), Millennials (1980-1996).
    # The 'TRUE' line is a catch-all that keeps the original year of birth for those not in these ranges.
    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(yrbrn)  
    ),
    # After recoding, the gen variable is converted into a factor with labels for clearer interpretation.
    # Factors are used in R to handle categorical variables.
    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 
##         5360         8583         5601         4652

# Inputting into the models
model3 <- lm(auth ~ polID + religious, data = df, weights = weight)
model4 <- lm(auth ~ polID + religious + polID*religious, data = df, weights = weight)

model3 <- lm(auth ~ polID + religious, data = df, weights = weight)
model4 <- lm(auth ~ polID + religious + polID*religious, data = df, weights = weight)

modelsummary(
  list(model3, model4),
  # you can add as many models as you have generated by putting them in this list. But you MUST refer to them to their stored name. So if you named your model say Regression1, you would have to put that instead of model1. By default, you can use the naming convention we rely upon of model1, model2, ...
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}", # what we are specifying here is what we want to show in our model summary table
                  "{estimate} ({std.error}){stars}"),
  statistic = NULL)

	(1)	(2)
(Intercept)	55.3 (0.3)***	56.4 (0.3)***
polIDModerate	4.6 (0.3)***	3.4 (0.4)***
polIDRight	7.7 (0.4)***	4.2 (0.7)***
religious	0.9 (0.0)***	0.6 (0.1)***
polIDModerate × religious		0.3 (0.1)***
polIDRight × religious		0.8 (0.1)***
Num.Obs.	23373	23373
R2	0.048	0.050
R2 Adj.	0.048	0.050
AIC	200594.1	200558.6
BIC	200634.4	200615.0
Log.Lik.	−100292.049	−100272.298
RMSE	17.27	17.24

olIDModerate (4.6): According to this coefficient, the average result for people who identify as politically moderate is 4.6 units higher than the group of reference (in our case, identifying ideologically as on the “left”) when comparing individuals who have the same religiousness (i.e., what some term “holding constant”). The category that is left out is referred to as the “baseline group” and is therefore considered to be a component of our intercept.
polIDRight (7.7): According to this coefficient, people who identify as politically right should, on average, score 7.7 units higher on the authoritarian values scale than the reference category, which is made up of people who identify as ideologically “left,” when comparing people of similar religiousness.
religious (0.9): According to the “religious” predictor, a coefficient of 1.6 indicates that, assuming that individuals are being compared who share the same ideological classification, an increase in religiousness of one unit is expected to result in an average increase of 0.9 units on the authoritarian values scale (i.e., holding ideological self-placement constant).

To sum up, The *** in each of the three coefficients indicates a p < 0.001. Meaning: Under frequentist assumptions, there is less than a 0.1% chance that the observed outcome is the result of chance, and hence we mistakenly reject the null.

Task 2 Now generate the model4 interaction plot that we did in the tutorial, but again using the German data instead of the French. Interpret.

germany_data <- read_fst("germany_data.fst")

interaction_plot <- effect("polID*religious", model4, na.rm=TRUE)

plot(interaction_plot,
     main="Interaction effect",
     xlab="Religiousness",
     ylab="Authoritarian attitudes scale")

# Displaying the interaction plot
interaction_plot <- effect("polID*religious", model4, na.rm=TRUE)

plot(interaction_plot,
     main="Interaction effect",
     xlab="Religiousness",
     ylab="Authoritarian attitudes scale")

interaction_plot

## 
##  polID*religious effect
##           religious
## polID             0        2        5        8       10
##   Left     56.40706 57.56775 59.30879 61.04983 62.21053
##   Moderate 59.85063 61.66176 64.37844 67.09513 68.90626
##   Right    60.60237 63.35032 67.47224 71.59416 74.34211

Similar to individuals who identify as conservative, the interaction plot for moderate political identity exhibits an increasing trend. This suggests that people with moderate political leanings generally rank higher on the authoritarian scale as they grow more devout. The slope of the two figures, however, indicates that for moderates, the influence of religion is less steep than for the Right group.

There is an increasing tendency in the interaction plot for those who identify as left-leaning politically. This suggests that Left-leaning people often score higher on the authoritarian scale as they grow more devout.There is a rising tendency in the interaction plot for those who identify as right-wing politically as well. This suggests that the average authoritarian score of Right-leaning people rises as they grow more devout.

Task 3 Use data for the Netherlands and the variable coding for model5 and model6 from the tutorial (i.e., both for the cleaning & recode, as well as the linear model formulas). Do a modelsummary table displaying both model outputs. Interpret the coefficients for the MLR without the interaction (i.e., the first displayed model), and interpret model fit metrics for both models.

# Reading the netherlands data fst
netherlands_data <- read_fst("netherlands_data.fst")

# Preparing weighted data set
netherlands_data$weight <- netherlands_data$dweight * netherlands_data$pweight
survey_design <- svydesign(ids = ~1, data = netherlands_data, weights = ~weight)

# Assigning the netherlands_data to 'df' so that it can be used instead
df <- netherlands_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
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))

df <- df %>%
  mutate(
    polID = case_when(
      lrscale %in% 0:3 ~ "Left",                    
      lrscale %in% 7:10 ~ "Right",                     
      lrscale %in% 4:6 ~ "Moderate",                  
      lrscale %in% c(77, 88, 99) ~ NA_character_      
    ),
   religious = case_when(
      rlgdgr %in% c(77, 88, 99) ~ NA_real_,
      TRUE ~ rlgdgr
    )
  )

df <- df %>%
  mutate(religion = case_when(
    rlgblg == 2 ~ "No",
    rlgblg == 1 ~ "Yes",
    rlgblg %in% c(7, 8, 9) ~ NA_character_,
    TRUE ~ as.character(rlgblg)
  ))

# check
table(df$religion)

## 
##    No   Yes 
## 10913  6750

# Code to produce a table of left and right 
df <- df %>%
  mutate(ID = case_when(
    lrscale >= 0 & lrscale <= 4 ~ "Left",
    lrscale >= 6 & lrscale <= 10 ~ "Right",
    lrscale > 10 ~ NA_character_,  # Set values above 10 as NA
    TRUE ~ NA_character_  # Ensure value 5 and any other unexpected values are set as NA
  ))
table(df$ID)

## 
##  Left Right 
##  5449  7010

# Code to showcase generations
df <- df %>%
  mutate(
    cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
    # Recoding generational cohorts based on the year of birth (yrbrn).
    # The year of birth is categorized into different generational cohorts.
    # Interwar (1900-1945), Baby Boomers (1946-1964), Gen X (1965-1979), Millennials (1980-1996).
    # The 'TRUE' line is a catch-all that keeps the original year of birth for those not in these ranges.
    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(yrbrn)  
    ),
    # After recoding, the gen variable is converted into a factor with labels for clearer interpretation.
    # Factors are used in R to handle categorical variables.
    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 
##         3771         6147         4548         2730

model5 <- lm(auth ~ religion + ID + gen, data = df, weights = weight)
model6 <- lm(auth ~ religion + ID + gen + religion*gen, data = df, weights = weight)

# Displaying a summary table
modelsummary(
  list(model5, model6),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_omit = "Intercept")

	(1)	(2)
religionYes	7.2 (0.3)***	5.7 (0.7)***
IDRight	4.2 (0.3)***	4.3 (0.3)***
genBaby Boomers	−5.2 (0.4)***	−5.9 (0.6)***
genGen X	−7.7 (0.4)***	−8.8 (0.6)***
genMillennials	−8.4 (0.5)***	−9.7 (0.6)***
religionYes × genBaby Boomers		1.1 (0.8)
religionYes × genGen X		2.3 (0.9)**
religionYes × genMillennials		2.9 (1.0)**
Num.Obs.	12116	12116
R2	0.124	0.125
R2 Adj.	0.124	0.125
AIC	101219.3	101213.5
BIC	101271.2	101287.5
Log.Lik.	−50602.672	−50596.757
RMSE	15.27	15.26

interaction_plot <- effect("religion*gen", model6, na.rm=TRUE)

plot(interaction_plot,
     main="Interaction effect",
     xlab="Generations",
     ylab="Authoritarian attitudes scale")

interaction_plot

## 
##  religion*gen effect
##         gen
## religion Interwar Baby Boomers    Gen X Millennials
##      No  66.39342     60.50327 57.62650    56.73392
##      Yes 72.13454     67.31264 65.70738    65.34644

cat_plot(model5, pred = gen, modx = religion, jnplot = TRUE)

## Warning: gen and religion are not included in an interaction with one another in the
## model.

cat_plot(model6, pred = gen, modx = religion, jnplot = TRUE)

religionYes (7.2): According to this coefficient, there is a 7.2 rise in authitorian values for those who identify as religiously inclined compared to those who do not.
IDRight (4.2): People who identify as belonging to the right wing of politics have a 4.2 in terms of authoritarian ideals.
genBaby Boomers (-5.2): The authoritarian values of this group have decreased by 5.2.
genGen X (-7.7): Uthoritarian values have decreased by 7.7, according to genGen X. -genMillenials (-8.4): genMillenails exhibits an 8.4 decline in liberal ideals.

Model 5: Overall Conclusion The triple asterisks next to each generational category (genBaby Boomers, genGen X, and genMillenials) indicate that the p-value is less than 0.005. This provides evidence to reject the null hypothesis by demonstrating the statistical significance of the explanatory variable.The relegionYes coefficient of 8.2 shows that there is an 8.2-point difference between the average authoritarian score and the identification of oneself as religious. The IDRright coefficient of 3.4 indicates that, on average, those who identify as politically right score 3.4 units higher than our reference group, which is made up of people who identify as left. R-squared corrected is 0.123. It emphasizes that 12.3% of the variation in the data can be explained by the model.

Model 6: Overall Conclusion The relegionYes coefficient of 7.5 indicates a 7.5-point increase in the average authoritarian score for those who identify as religious. The IDRright coefficient of 3.4 indicates that, on average, those who identify as politically right score 3.4 units higher than our reference group, which is made up of people who identify as left. 12.1% of the variance in the data is explained by the model, according to the corrected R-squared of 0.121. Each coefficient’s *** indicates that, under the most prevalent assumptions, the p-value is less than 0.005, indicating that the coefficients are statistically significant.

Task 4 Produce the model7 interaction plot from the tutorial, but again using the Netherlands data. Interpret.

netherlands_data <- netherlands_data %>%
  mutate(
    polID = case_when(
      lrscale %in% 0:3 ~ "Left",                    
      lrscale %in% 7:10 ~ "Right",                     
      lrscale %in% 4:6 ~ "Moderate",                  
      lrscale %in% c(77, 88, 99) ~ NA_character_      
    ),
   religious = case_when(
      rlgdgr %in% c(77, 88, 99) ~ NA_real_,
      TRUE ~ rlgdgr
    )
  )

# Displaying model 7 table with interaction plot
model7 <- lm(auth ~ + cohort + polID + cohort*polID, data = df, weights = weight)
interact_plot(model7, pred = cohort, modx = polID, jnplot = TRUE)

The authoritarian scale scores for the three cohorts (PolID left, POlID right, and PolID moderates) show a declining tendency over time. In contrast, the political right has the sharpest slope and the largest decline in average authoritarian scales. The moderates come after this. The group who identified as politically left experienced the least amount of decline overall. All three groups have a similar overall trend of declining authoritarian values, though.

The lines converge, demonstrating that there is interaction between variables.

Task 5 Produce the model8 interaction plot from the tutorial, but again using the Netherlands data. Interpret.

# Displaying model8 with interaction plot
model8 <- lm(auth ~ religious + ID + religious*ID, data = df, weights = weight)
interact_plot(model8, pred = religious, modx = ID, jnplot = TRUE)

As religious identification rises, authoritarian attitudes rise within both the ID left and ID right groups to a relatively comparable degree, as seen by the generally favorable trend for both groups.

The greater intercept of the ID right cohort indicates that, in the absence of religion as an explanation variable, the ID right cohort would, on average, score higher on the authoritarian scale than the ID left cohort. There is no contact since the two lines are almost parallel rather than convergent.

END

Shah_Kashvi_Homework5

Kashvi Shah

2024-03-26