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

Instructions: Start a new R markdown for the homework and call it “Yourlastname_Firstname_Homework_5”.

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.

packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "MASS", "aod", "interactions", "kableExtra", "flextable", "scales")

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
## 
## 
## 
## Attaching package: 'MASS'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 
## 
## Attaching package: 'aod'
## 
## 
## The following object is masked from 'package:survival':
## 
##     rats
## 
## 
## 
## Attaching package: 'kableExtra'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
## 
## 
## 
## Attaching package: 'flextable'
## 
## 
## The following objects are masked from 'package:kableExtra':
## 
##     as_image, footnote
## 
## 
## The following object is masked from 'package:purrr':
## 
##     compose
## 
## 
## 
## Attaching package: 'scales'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
## [[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] "MASS"         "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] "aod"          "MASS"         "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"        
## 
## [[9]]
##  [1] "interactions" "aod"          "MASS"         "survey"       "survival"    
##  [6] "Matrix"       "grid"         "effects"      "carData"      "modelsummary"
## [11] "fst"          "infer"        "lubridate"    "forcats"      "stringr"     
## [16] "dplyr"        "purrr"        "readr"        "tidyr"        "tibble"      
## [21] "ggplot2"      "tidyverse"    "stats"        "graphics"     "grDevices"   
## [26] "utils"        "datasets"     "methods"      "base"        
## 
## [[10]]
##  [1] "kableExtra"   "interactions" "aod"          "MASS"         "survey"      
##  [6] "survival"     "Matrix"       "grid"         "effects"      "carData"     
## [11] "modelsummary" "fst"          "infer"        "lubridate"    "forcats"     
## [16] "stringr"      "dplyr"        "purrr"        "readr"        "tidyr"       
## [21] "tibble"       "ggplot2"      "tidyverse"    "stats"        "graphics"    
## [26] "grDevices"    "utils"        "datasets"     "methods"      "base"        
## 
## [[11]]
##  [1] "flextable"    "kableExtra"   "interactions" "aod"          "MASS"        
##  [6] "survey"       "survival"     "Matrix"       "grid"         "effects"     
## [11] "carData"      "modelsummary" "fst"          "infer"        "lubridate"   
## [16] "forcats"      "stringr"      "dplyr"        "purrr"        "readr"       
## [21] "tidyr"        "tibble"       "ggplot2"      "tidyverse"    "stats"       
## [26] "graphics"     "grDevices"    "utils"        "datasets"     "methods"     
## [31] "base"        
## 
## [[12]]
##  [1] "scales"       "flextable"    "kableExtra"   "interactions" "aod"         
##  [6] "MASS"         "survey"       "survival"     "Matrix"       "grid"        
## [11] "effects"      "carData"      "modelsummary" "fst"          "infer"       
## [16] "lubridate"    "forcats"      "stringr"      "dplyr"        "purrr"       
## [21] "readr"        "tidyr"        "tibble"       "ggplot2"      "tidyverse"   
## [26] "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
## [31] "methods"      "base"

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.

germany_data <- read_fst("germany_data.fst")
gd <- germany_data
gd <- gd %>%
  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
gd$auth <- scales::rescale(gd$behave + 
                      gd$secure + 
                      gd$safety + 
                      gd$tradition + 
                      gd$rules, to=c(0,100), na.rm=TRUE)


gd <- gd %>% filter(!is.na(auth))
gd <- gd %>%
  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
    )
  )
model3 <- lm(auth ~ polID + religious, data = gd, weights = weight)
model4 <- lm(auth ~ polID + religious + polID*religious, data = gd, weights = weight)
modelsummary(
  list(model3, model4),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL)
 (1)   (2)
(Intercept) 57.5 (0.7)*** 59.6 (0.9)***
polIDModerate 5.3 (0.7)*** 2.6 (1.1)*
polIDRight 9.3 (1.1)*** 3.4 (1.9)+
religious 0.3 (0.1)** −0.3 (0.2)+
polIDModerate × religious 0.8 (0.2)**
polIDRight × religious 1.5 (0.4)***
Num.Obs. 2855 2855
R2 0.035 0.041
R2 Adj. 0.034 0.039
AIC 24680.0 24665.5
BIC 24709.7 24707.2
Log.Lik. −12334.978 −12325.764
RMSE 17.13 17.13

polIDModerate(0.7): With a coefficient of 0.7, the graph displays that Germans identifying as moderate politically is 0.7 points higher than the unnamed identifier (left).

polIDRight(1.1): With the table data, Germans aligning right politically is 1.1 points higher than identifying as left wing.

religious(0.1): people who are religious, and identify with either of the two variables (moderate and right) are 0.1 points higher than those who do not alight politically or religiously with those variables (non-religious and left wing)

Two of the three coefficients have a *** marking next to their coefficients meaning that the p value result being due to chance is quite low. the last coefficient (religious) has ** markings beside the coefficient, showing that there is a higher probability of the results being due to chance and incorrectly rejecting 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.

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     59.61010 58.93398 57.91980 56.90562 56.22950
##   Moderate 62.21299 63.14317 64.53843 65.93370 66.86388
##   Right    63.00456 65.26702 68.66072 72.05442 74.31689

Takeaways:

What I notice about these plots is that the political affiliation for being considered left wing is an example of interaction as that line is non parallel to the rest.

The rest of the plot lines in the moderate category and right category are more parallel to each other, but if the lines were to continue outwards they would eventually intersect.

These three plot lines show that there is interactions between the variables due to their non parallel lines

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.

netherlands_data <- read_fst("netherlands_data.fst")
nd <- netherlands_data
nd <- nd %>%
  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 ))


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


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

table(nd$religion)
## 
##    No   Yes 
## 10913  6750
nd <- nd %>%
  mutate(ID = case_when(
    lrscale >= 0 & lrscale <= 4 ~ "Left",
    lrscale >= 6 & lrscale <= 10 ~ "Right",
    lrscale > 10 ~ NA_character_, 
    TRUE ~ NA_character_  
  ))
table(nd$ID)
## 
##  Left Right 
##  5449  7010
model5 <- lm(auth ~ religion + ID + behave, data = nd, weights = weight)
model6 <- lm(auth ~ religion + ID + behave + religion*behave, data = nd, weights = weight)
modelsummary(
  list(model5, model6),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_omit = "Intercept")
 (1)   (2)
religionYes 5.9 (0.7)*** 8.5 (2.6)**
IDRight 2.6 (0.6)*** 2.6 (0.6)***
behave 9.6 (0.3)*** 9.8 (0.4)***
religionYes × behave −0.6 (0.6)
Num.Obs. 1176 1176
R2 0.533 0.534
R2 Adj. 0.532 0.532
AIC 8983.2 8984.1
BIC 9008.5 9014.6
Log.Lik. −4486.598 −4486.068
RMSE 10.91 10.91

religionYes (0.7): A coefficient of 0.7 displays that there is a 0.7 point increase for individuals who are religious compared to individuals who are not religious.

IDRight (2.6): Those who identify as religious are most likely to identify as right winged by 2.6 points on the authoritarian scale than those who identify as the un-referenced/omitted category.(Left, moderate, and non-religious).

behave(9.6): for the authoritarian scale, a coefficient of 9.6 indicates that individuals are most likely to be associated with authoritarian behaviours and discipline methods.

All three coefficients have a *** marking meaning that the p-values are extremely low and the probability of an error is unlikely. Meaning that there is a less than 0.1% chance that an error has occurred due to chance, and the rejection of the null due to those circumstances.

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

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

interaction_plot
## 
##  religion*behave effect
##         behave
## religion        1        2        4        5        6
##      No  29.92494 39.75404 59.41224 69.24134 79.07045
##      Yes 37.83663 47.04872 65.47289 74.68497 83.89706
interaction_plot
## 
##  religion*behave effect
##         behave
## religion        1        2        4        5        6
##      No  29.92494 39.75404 59.41224 69.24134 79.07045
##      Yes 37.83663 47.04872 65.47289 74.68497 83.89706
cat_plot(model6, pred = behave, modx = religion, jnplot = TRUE)

Looking at the way the graph is displayed, the lines show two relatively parallel lines indicating no interaction between the two. if you let the lines continue upwards towards a higher rate for authoritarian values, the lines will inevitably intersect along the way.

Task 4

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

the codes i was running were working fine before i started to knit the document, and then it was giving me error messages. i’m not quite sure what i’m doing wrong or how to fix them to get it to knit, so here’s the work and observations that i had when it was working

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


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


nd <- nd %>% filter(!is.na(auth))
interaction_plot
## 
##  religion*behave effect
##         behave
## religion        1        2        4        5        6
##      No  29.92494 39.75404 59.41224 69.24134 79.07045
##      Yes 37.83663 47.04872 65.47289 74.68497 83.89706

interact_plot(model7, pred = cohort, modx = polID, jnplot = TRUE) interaction_effect <- effects::effect(“cohort*polID”, model7, na.rm = TRUE)

plot(interaction_plot, main=“Interaction effect”, xlab=“cohort”, ylab=“Authoritarian attitudes scale”)

Task 5

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

netherlands_data <- read_fst(“netherlands_data.fst”) nd <- netherlands_data

interaction_plot <- effect(“religious*ID”, model8, na.rm=TRUE)

plot(interaction_plot, main=“Interaction effect”, xlab=““, ylab=”Authoritarian attitudes scale”)

interaction_plot

End