# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "performance", "flextable", "broom", "scales", "ggeffects", "marginaleffects") # 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
## 
## 
## 
## Attaching package: 'flextable'
## 
## 
## 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] "performance"  "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] "flextable"    "performance"  "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] "broom"        "flextable"    "performance"  "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] "scales"       "broom"        "flextable"    "performance"  "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] "ggeffects"    "scales"       "broom"        "flextable"    "performance" 
##  [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] "marginaleffects" "ggeffects"       "scales"          "broom"          
##  [5] "flextable"       "performance"     "survey"          "survival"       
##  [9] "Matrix"          "grid"            "effects"         "carData"        
## [13] "modelsummary"    "fst"             "infer"           "lubridate"      
## [17] "forcats"         "stringr"         "dplyr"           "purrr"          
## [21] "readr"           "tidyr"           "tibble"          "ggplot2"        
## [25] "tidyverse"       "stats"           "graphics"        "grDevices"      
## [29] "utils"           "datasets"        "methods"         "base"
ess <- read_fst("All-ESS-Data.fst")

##Introduction

france_data <- ess %>%
  filter(cntry == "FR")
france_data$weight <- france_data$dweight * france_data$pweight
survey_design <- svydesign(ids = ~1, data = france_data, weights = ~weight)
france_data <- france_data %>%
  mutate(protest = case_when(
    trstplc == 1 ~ 1,    
    trstplc == 2 ~ 0,    
    TRUE ~ NA_real_    
  ))
france_data <- france_data %>% filter(!is.na(trstplc))

Potential Prodictors

france_data <- france_data %>%
  
  mutate(feelincome = case_when(
    hincfel == 1 ~ "comfortable",  
    hincfel %in% 2:4 ~ "difficult", 
    TRUE ~ NA_character_            
  )) %>%
  
  mutate(
    age = ifelse(agea == 999, NA, agea),
    adulthood = ifelse(age >= 18 & age <= 30, "Early", ifelse(age >= 31, "Later", NA)),
    cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
    gen = case_when(
      yrbrn %in% 1900:1945 ~ "Interwar",
      yrbrn %in% 1946:1964 ~ "Baby Boomers",
      yrbrn %in% 1965:1979 ~ "Gen X",
      yrbrn %in% 1980:1996 ~ "Millennials",
      TRUE ~ as.character(cohort)
    ),
    gen = factor(gen, levels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
  ) %>%
 
  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_      
    ),
    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"))
  )

france_data$trstplc <- ifelse(france_data$trstplc > 10, NA, france_data$trstplc)


france_data$trust <- scales::rescale(france_data$trstplc, na.rm = TRUE, to = c(0, 100))


france_data$populist <- scales::rescale(france_data$trust, na.rm = TRUE, to=c(100,0))
table(france_data$populist)
## 
##    0   10   20   30   40   50   60   70   80   90  100 
##  730 1062 3588 3741 2679 3271 1328 1064  690  299  539

##Statisctics for Potential Predictors and shwos or Outcome Variable Trust

FRancedf <- france_data %>%
  dplyr::select(trust, feelincome, adulthood, gen, populist, polID, educ.ba, weight)


FRancedf <- na.omit(FRancedf)
head(FRancedf)
##      trust  feelincome adulthood          gen populist    polID    educ.ba
## 3310    30   difficult     Later        Gen X       70 Moderate      No BA
## 3312    60   difficult     Later        Gen X       40 Moderate      No BA
## 3314    60 comfortable     Later Baby Boomers       40     Left BA or more
## 3315    40 comfortable     Later Baby Boomers       60    Right      No BA
## 3316    50 comfortable     Later Baby Boomers       50    Right BA or more
## 3317    60 comfortable     Later        Gen X       40 Moderate BA or more
##        weight
## 3310 1.233645
## 3312 2.151479
## 3314 2.930847
## 3315 3.060282
## 3316 1.372110
## 3317 2.211179

The Outcome Variable and Predictors for Research Report

FRancedf <- france_data %>%
  mutate(trust = case_when(
    trstplc == 1 ~ "Yes",    
    trstplc == 0 ~ "No",     
    TRUE ~ NA_character_     
  ))

FRancedf <- FRancedf %>%
  
  mutate(
    TrustPolice = trust,      
    Generation = gen,      
    Ideology = polID
  )
tablee1 <- datasummary_skim(france_data  %>% dplyr::select(trust, gen , polID), type = "categorical",
                           title = "Table 1. Descriptive statistics for main variables", output = "flextable")

tablee1
Table 1. Descriptive statistics for main variables

N

%

gen

Interwar

3999

21.0

Baby Boomers

6351

33.4

Gen X

4821

25.3

Millennials

3308

17.4

polID

Left

4709

24.7

Moderate

8706

45.7

Right

4302

22.6

NA

1321

6.9

Hypothesis Testing

##STEP 1:

test_stat <- FRancedf %>%
  specify(explanatory = trust, 
          response = gen ) %>% 
  hypothesize(null = "independence") %>%
  calculate(stat = "Chisq") 
## Warning: Removed 18223 rows containing missing values.
print(test_stat$stat)
## X-squared 
##  5.843781

##STEP 2:

null_dist <- FRancedf %>%
  specify(response = gen, explanatory = trust) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "Chisq")
## Warning: Removed 18223 rows containing missing values.
null_dist
## Response: gen (factor)
## Explanatory: trust (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1 4.21 
##  2         2 1.02 
##  3         3 2.75 
##  4         4 2.31 
##  5         5 2.61 
##  6         6 2.82 
##  7         7 2.53 
##  8         8 0.435
##  9         9 3.64 
## 10        10 4.77 
## # ℹ 990 more rows

##STEP 3:

p_val <- null_dist %>% # 
  get_pvalue(obs_stat = test_stat, direction = "greater")
p_val
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1   0.118

##STEP 4:

conf_int <- null_dist%>%
  get_confidence_interval(level = 0.95, type = "percentile")


null_dist %>%
  visualize(data, bins = 10, method = "simulation", dens_color = "black") +
  shade_p_value(obs_stat = test_stat, direction = "greater") +
  shade_confidence_interval(endpoints = conf_int)

## Additional Test Stat Statistic

test_stat <- 1000

conf_int <- null_dist %>%
  get_confidence_interval(level = 0.95, type = "percentile")


null_dist %>%
  visualize(bins = 15, method = "simulation", dens_color = "black") +
  geom_vline(aes(xintercept = test_stat), color = "red", linetype = "dashed", linewidth = 1) +
  shade_p_value(obs_stat = test_stat, direction = "greater") +
  shade_confidence_interval(endpoints = conf_int) +
  coord_cartesian(xlim = c(0, max(null_dist$stat) * 1.1)) + # Set manual x-axis limits
  labs(title = "Simulation-Based Null Distribution",
       x = "Chi-Square",
       y = "Count")
## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf

##Regression Models

FRancedf <- FRancedf %>%
  mutate(
    trust = case_when(
      trust == "Yes" ~ 1,
      trust == "No" ~ 0,
      TRUE ~ NA_integer_ 
    )
  )
modell1 <- glm(trust ~ Generation, family = binomial, data = FRancedf, weights = round(weight))   
modell2 <- glm(trust ~ Generation + Ideology, family = binomial, data = FRancedf, weights = round(weight))  
modell3 <- glm(trust ~ Generation + Ideology + Generation*Ideology, family = binomial, data = FRancedf, weights = round(weight)) 
modelsummary(
  list(modell1, modell2, modell3),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL)
 (1)   (2)   (3)
(Intercept) −0.4 (0.1)** −0.3 (0.1)* −0.5 (0.2)*
GenerationBaby Boomers 0.0 (0.1) −0.1 (0.2) 0.1 (0.3)
GenerationGen X −0.2 (0.2) −0.3 (0.2)+ 0.0 (0.3)
GenerationMillennials −0.5 (0.2)*** −0.7 (0.2)*** −0.5 (0.3)+
IdeologyModerate 0.0 (0.1) 0.4 (0.3)
IdeologyRight 0.0 (0.1) 0.1 (0.3)
GenerationBaby Boomers × IdeologyModerate −0.3 (0.4)
GenerationGen X × IdeologyModerate −0.5 (0.4)
GenerationMillennials × IdeologyModerate −0.7 (0.4)+
GenerationBaby Boomers × IdeologyRight −0.2 (0.4)
GenerationGen X × IdeologyRight −0.6 (0.4)
GenerationMillennials × IdeologyRight 0.5 (0.4)
Num.Obs. 814 717 717
AIC 2973.2 2639.2 2634.9
BIC 2992.0 2666.7 2689.8
Log.Lik. −1482.613 −1313.611 −1305.428
RMSE 0.48 0.48 0.48

##Linear Regeression Model Output Table “SLR”, “MLR” and “Interation Between Both Predictors”

models <- list()
models[['SLR']]<- glm(trust ~ Generation, family = binomial, data = FRancedf, weights = round(weight))   
models[['MLR']]<- glm(trust ~ Generation + Ideology, family = binomial, data = FRancedf, weights = round(weight))  
models[['Interaction']]<- glm(trust ~ Generation + Ideology + Generation*Ideology, family = binomial, data = FRancedf, weights = round(weight)) 
modelsummary(models,
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  exponentiate = FALSE, ## NOTE: this is to exponentiate directly, otherwise leave out or set to FALSE
  coef_rename = c("GenerationBaby Boomers" = "Baby Boomers", "GenerationGen X" = "Gen X", "GenerationMillennials" = "Millennials", "IdeologyExtreme Left" = "Extreme Left", "IdeologyExtreme Right" = "Extreme Right", "IdeologyModerate Left" = "Moderate Left", "IdeologyModerate Right" = "Moderate Right"),
  title = 'Table 2. Regression models predicting trust in France')
Table 2. Regression models predicting trust in France
SLR MLR  Interaction
(Intercept) −0.4 (0.1)** −0.3 (0.1)* −0.5 (0.2)*
Baby Boomers 0.0 (0.1) −0.1 (0.2) 0.1 (0.3)
Gen X −0.2 (0.2) −0.3 (0.2)+ 0.0 (0.3)
Millennials −0.5 (0.2)*** −0.7 (0.2)*** −0.5 (0.3)+
IdeologyModerate 0.0 (0.1) 0.4 (0.3)
IdeologyRight 0.0 (0.1) 0.1 (0.3)
Baby Boomers:IdeologyModerate −0.3 (0.4)
Gen X:IdeologyModerate −0.5 (0.4)
Millennials:IdeologyModerate −0.7 (0.4)+
Baby Boomers:IdeologyRight −0.2 (0.4)
Gen X:IdeologyRight −0.6 (0.4)
Millennials:IdeologyRight 0.5 (0.4)
Num.Obs. 814 717 717
AIC 2973.2 2639.2 2634.9
BIC 2992.0 2666.7 2689.8
Log.Lik. −1482.613 −1313.611 −1305.428
RMSE 0.48 0.48 0.48

My Own Visual of Choice

model4 <- glm(trust ~ lrscale + Ideology + lrscale*Ideology, family = binomial, data = FRancedf, weights = round(weight)) # both predictors + interaction
mydf <- predict_response(model4, terms = c("lrscale", "Ideology"))
print(mydf, collapse_tables = TRUE)
## # Predicted probabilities of trust
## 
## lrscale | Ideology | Predicted |     95% CI
## -------------------------------------------
##       0 |     Left |      0.32 | 0.28, 0.37
##       2 |          |      0.38 | 0.34, 0.41
##       3 |          |      0.41 | 0.35, 0.46
##       5 |          |      0.46 | 0.36, 0.57
##       7 |          |      0.52 | 0.36, 0.68
##      10 |          |      0.61 | 0.37, 0.81
##       0 | Moderate |      0.40 | 0.16, 0.69
##       2 |          |      0.38 | 0.23, 0.56
##       3 |          |      0.37 | 0.26, 0.49
##       5 |          |      0.35 | 0.32, 0.38
##       7 |          |      0.33 | 0.22, 0.46
##      10 |          |      0.31 | 0.11, 0.61
##       0 |    Right |      0.87 | 0.55, 0.97
##       2 |          |      0.78 | 0.51, 0.93
##       3 |          |      0.73 | 0.48, 0.89
##       5 |          |      0.61 | 0.44, 0.75
##       7 |          |      0.46 | 0.38, 0.55
##      10 |          |      0.27 | 0.20, 0.35
## 
## Not all rows are shown in the output. Use `print(..., n = Inf)` to show
##   all rows.
ggplot(mydf, aes(x, predicted, colour = group)) + geom_line()

##Another Visual of Choice for Research Report

plot_comparisons(modell3, variables = "Generation", condition = "Ideology")

##Checking My Model Fit

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


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


FRancedf <- FRancedf %>% filter(!is.na(auth))
FRancedf <- FRancedf %>%
  mutate(

    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)
    ),
    gen = factor(gen, levels = c("1", "2", "3", "4"),
                 labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials")),
    LRscale = case_when(
      lrscale %in% c(77, 88, 99) ~ NA_real_,
      TRUE ~ lrscale
    )
  )
modell5 <- lm(auth ~ gen + lrscale + gen*lrscale, data = FRancedf, weights = weight)
check_model(modell5)