packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "MASS", "aod", "interactions", "kableExtra", "flextable", "scales", "ggeffects", "marginaleffects", "viridis", "knitr", "rmarkdown", "ggridges", "questionr")
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)
lapply(packages, library, character.only = TRUE)## Warning: package 'tidyverse' was built under R version 4.3.3## ── 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## Warning: package 'infer' was built under R version 4.3.3## Warning: package 'effects' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Warning: package 'survey' was built under R version 4.3.3## 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## Warning: package 'aod' was built under R version 4.3.3##
## Attaching package: 'aod'
##
## The following object is masked from 'package:survival':
##
## rats## Warning: package 'interactions' was built under R version 4.3.3##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows## Warning: package 'flextable' was built under R version 4.3.3##
## 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## Warning: package 'ggeffects' was built under R version 4.3.3##
## Attaching package: 'ggeffects'
##
## The following object is masked from 'package:interactions':
##
## johnson_neyman## Warning: package 'marginaleffects' was built under R version 4.3.3## Loading required package: viridisLite
##
## Attaching package: 'viridis'
##
## The following object is masked from 'package:scales':
##
## viridis_pal## Warning: package 'ggridges' was built under R version 4.3.3## Warning: package 'questionr' was built under R version 4.3.3## [[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"
##
## [[13]]
## [1] "ggeffects" "scales" "flextable" "kableExtra" "interactions"
## [6] "aod" "MASS" "survey" "survival" "Matrix"
## [11] "grid" "effects" "carData" "modelsummary" "fst"
## [16] "infer" "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [26] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [31] "datasets" "methods" "base"
##
## [[14]]
## [1] "marginaleffects" "ggeffects" "scales" "flextable"
## [5] "kableExtra" "interactions" "aod" "MASS"
## [9] "survey" "survival" "Matrix" "grid"
## [13] "effects" "carData" "modelsummary" "fst"
## [17] "infer" "lubridate" "forcats" "stringr"
## [21] "dplyr" "purrr" "readr" "tidyr"
## [25] "tibble" "ggplot2" "tidyverse" "stats"
## [29] "graphics" "grDevices" "utils" "datasets"
## [33] "methods" "base"
##
## [[15]]
## [1] "viridis" "viridisLite" "marginaleffects" "ggeffects"
## [5] "scales" "flextable" "kableExtra" "interactions"
## [9] "aod" "MASS" "survey" "survival"
## [13] "Matrix" "grid" "effects" "carData"
## [17] "modelsummary" "fst" "infer" "lubridate"
## [21] "forcats" "stringr" "dplyr" "purrr"
## [25] "readr" "tidyr" "tibble" "ggplot2"
## [29] "tidyverse" "stats" "graphics" "grDevices"
## [33] "utils" "datasets" "methods" "base"
##
## [[16]]
## [1] "knitr" "viridis" "viridisLite" "marginaleffects"
## [5] "ggeffects" "scales" "flextable" "kableExtra"
## [9] "interactions" "aod" "MASS" "survey"
## [13] "survival" "Matrix" "grid" "effects"
## [17] "carData" "modelsummary" "fst" "infer"
## [21] "lubridate" "forcats" "stringr" "dplyr"
## [25] "purrr" "readr" "tidyr" "tibble"
## [29] "ggplot2" "tidyverse" "stats" "graphics"
## [33] "grDevices" "utils" "datasets" "methods"
## [37] "base"
##
## [[17]]
## [1] "rmarkdown" "knitr" "viridis" "viridisLite"
## [5] "marginaleffects" "ggeffects" "scales" "flextable"
## [9] "kableExtra" "interactions" "aod" "MASS"
## [13] "survey" "survival" "Matrix" "grid"
## [17] "effects" "carData" "modelsummary" "fst"
## [21] "infer" "lubridate" "forcats" "stringr"
## [25] "dplyr" "purrr" "readr" "tidyr"
## [29] "tibble" "ggplot2" "tidyverse" "stats"
## [33] "graphics" "grDevices" "utils" "datasets"
## [37] "methods" "base"
##
## [[18]]
## [1] "ggridges" "rmarkdown" "knitr" "viridis"
## [5] "viridisLite" "marginaleffects" "ggeffects" "scales"
## [9] "flextable" "kableExtra" "interactions" "aod"
## [13] "MASS" "survey" "survival" "Matrix"
## [17] "grid" "effects" "carData" "modelsummary"
## [21] "fst" "infer" "lubridate" "forcats"
## [25] "stringr" "dplyr" "purrr" "readr"
## [29] "tidyr" "tibble" "ggplot2" "tidyverse"
## [33] "stats" "graphics" "grDevices" "utils"
## [37] "datasets" "methods" "base"
##
## [[19]]
## [1] "questionr" "ggridges" "rmarkdown" "knitr"
## [5] "viridis" "viridisLite" "marginaleffects" "ggeffects"
## [9] "scales" "flextable" "kableExtra" "interactions"
## [13] "aod" "MASS" "survey" "survival"
## [17] "Matrix" "grid" "effects" "carData"
## [21] "modelsummary" "fst" "infer" "lubridate"
## [25] "forcats" "stringr" "dplyr" "purrr"
## [29] "readr" "tidyr" "tibble" "ggplot2"
## [33] "tidyverse" "stats" "graphics" "grDevices"
## [37] "utils" "datasets" "methods" "base"france_data <- read.fst("france_data.fst")#context
CONTEXT 1. BY YEAR
francedf <- france_data %>%
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_,
TRUE ~ as.character(gndr)
)
)
table(francedf$gndr)##
## Female Male
## 10211 8827france_data <- francedf %>%
mutate(trust_category = case_when(
trstplc >= 0 & trstplc <= 3 ~ "Low trust",
trstplc >= 4 & trstplc <= 7 ~ "Medium trust",
trstplc >= 8 & trstplc <= 10 ~ "High trust",
TRUE ~ NA_character_ # Handle possible NA or out-of-range values
))
# Check the distribution of trust categories
table(france_data$trust_category)##
## High trust Low trust Medium trust
## 5380 2592 11019df <- france_data
df$year <- NA
replacements <- c(2001, 2003, 2005, 2007, 2009, 2011, 2013, 2015, 2017, 2019) #setting the years to odd numbers to determine the income.
for(i in 1:10){
df$year[df$essround == i] <- replacements[i]
}
trstplc_avg <- df %>%
group_by(year) %>%
summarize(trstplc_avg = mean(trstplc, na.rm = TRUE))
trstplc_plot_FR <- ggplot(trstplc_avg, aes(x = year, y = trstplc_avg)) +
geom_point(aes(size = trstplc_avg, color = trstplc_avg), alpha = 0.7) +
geom_line(aes(group = 1), color = "green", linetype = "dashed") +
labs(
x = "Survey Year",
y = "Average Trust in Police",
title = "Trend of Average Trust in Police Over Time in France") + # Added title here +
theme_minimal() +
theme(legend.position = "none")
# Y axis has income and x axis has the year. Years are arranged by odd numbers. With the highest year being 2019
print(trstplc_plot_FR)
- DESCRIPTIVE STATISTICS
france_data$weight <- france_data$dweight * france_data$pweight
survey_design <- svydesign(ids = ~1, data = france_data, weights = ~weight)france_descriptive <- france_data %>%
mutate(
trstplc = case_when (
trstplc %in% 0:2 ~ "No trust at all",
trstplc %in% 3:6 ~ "Medium Trust",
trstplc %in% 7:10 ~ "Complete trust",
trstplc %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(trstplc)
)
) %>%
filter(!is.na(trstplc))france_descriptive <- france_data %>%
mutate(
trstplc = ifelse(trstplc %in% c(7, 8, 9), NA, trstplc)
) %>%
filter(!is.na(trstplc))france_descriptive <- france_descriptive %>%
mutate(
gndr = ifelse(gndr %in% c(9), NA, gndr)
) %>%
filter(!is.na(gndr))france_descriptive <- france_descriptive %>%
mutate(
edulvlb = ifelse(edulvlb %in% c(77,88,99), NA, edulvlb)) %>%
filter(!is.na(edulvlb))france_descriptive <- france_descriptive %>%
mutate(
Environment = trstplc,
Gender = gndr,
Education = edulvlb
)france_table1 <- datasummary_skim(france_descriptive %>% dplyr::select(trstplc, gndr, edulvlb), title = "Table 1. Descriptive statistics of Trust in police, Gender, and Education", output = "flextable")## Warning: The histogram argument is only supported for (a) output types "default",
## "html", "kableExtra", or "gt"; (b) writing to file paths with extensions
## ".html", ".jpg", or ".png"; and (c) Rmarkdown, knitr or Quarto documents
## compiled to PDF (via kableExtra) or HTML (via kableExtra or gt). Use
## `histogram=FALSE` to silence this warning.france_table1
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
trstplc | 10 | 0 | 5.2 | 6.0 | 0.0 | 5.0 | 88.0 |
edulvlb | 18 | 0 | 376.6 | 376.9 | 0.0 | 321.0 | 8888.0 |
FINDINGS 1.NULL DISTRIBUTION
france_data <- france_data %>%
mutate(
# Clean 'trstplc' values
trstplc = ifelse(trstplc %in% c(77, 88, 99), NA, trstplc),
# Recode 'gndr'
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr == 9 ~ NA_character_, # Assuming '9' is a code for missing data
TRUE ~ as.character(gndr)
),
# Recode 'edulvla' to 'educ.ba'
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more", # Assuming '5' indicates bachelor's degree or higher for rounds < 5
TRUE ~ "No BA" # All other cases are considered 'No BA'
),
# Clean 'edulvla' values
edulvla = ifelse(edulvla %in% c(77, 88), NA_integer_, edulvla)
) %>%
mutate(
# Force 'educ.ba' to be treated as character to avoid factor conversion issues in ggplot
educ.ba = as.character(educ.ba)
) %>%
# Remove rows with missing data for the analysis
filter(!is.na(gndr), !is.na(educ.ba), !is.na(trstplc))# Specify the data
spec <- france_data %>%
specify(response = trstplc, explanatory = gndr)
# Perform the t-test
test_stat <- hypothesize(spec, null = "independence") %>%
calculate(stat = "t")## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "Female" - "Male", or divided in the order "Female" / "Male" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Female", "Male")` to the calculate() function.# Print the test statistic
print(test_stat)## Response: trstplc (numeric)
## Explanatory: gndr (factor)
## Null Hypothesis: independence
## # A tibble: 1 × 1
## stat
## <dbl>
## 1 4.83STEP 1: Calculating test statistic for our sample
test_stat <- france_data %>%
specify(explanatory = gndr, # change variable name for explanatory variable
response = trstplc) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "diff in means", order = c("Female", "Male"))## Message: The independence null hypothesis does not inform calculation of the
## observed statistic (a difference in means) and will be ignored.print(test_stat$stat)## [1] 0.1580604STEP 2: Simulating null distribution
null_distribution <- france_data %>%
specify(response = trstplc, explanatory = gndr) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "diff in means", order = c("Female", "Male"))STEP 3: Calculating p-value
p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "two-sided")## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation
## based on the number of `reps` chosen in the `generate()` step.
## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information.p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0Step 4: Visualization
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided")
null_distribution## Response: trstplc (numeric)
## Explanatory: gndr (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.0181
## 2 2 -0.0162
## 3 3 -0.0674
## 4 4 -0.0361
## 5 5 -0.0158
## 6 6 -0.0630
## 7 7 0.0558
## 8 8 0.0185
## 9 9 -0.0185
## 10 10 0.0842
## # ℹ 990 more rows##Interpreation: The values here some have a lower value than 0.05 and very few are larger than 0.05 therefore more suggest strong eveidence to reject the null hypothesis while there are some do not.
conf_int <- null_distribution %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int)
- REGRESSION MODEL
france_model1 <- france_data %>%
mutate(trstplc = na_if(trstplc, 7) %>% na_if(8) %>% na_if(9)) %>%
mutate(trstplc = 7 - trstplc) %>%
filter(!is.na(trstplc))france_data <- france_data %>%
# Recoding education with case_when for consistency
mutate(
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
# Ensure all NAs are handled uniformly
mutate(
edulvla = ifelse(edulvla %in% c(77, 88), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba) # Force educ.ba to be treated as character
)france_model1 <- france_model1 %>%
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr %in% c(9) ~ NA_character_,
TRUE ~ as.character(gndr) )
) %>%
mutate(educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
# Ensure all NAs are handled uniformly
mutate(
edulvla = ifelse(edulvla %in% c(77, 88), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba) # Force educ.ba to be treated as character
) %>%
filter(!is.na(gndr))table(france_model1$trstplc)##
## -3 1 2 3 4 5 6 7
## 730 2679 3271 1328 1064 690 299 539table(france_model1$educ.ba)##
## BA or more No BA
## 2124 8476table(france_model1$gndr)##
## Female Male
## 5635 4965france_clean <- france_data %>%
mutate(
trstplc= case_when(
trstplc%in% c(7,8,9) ~ NA_character_,
TRUE ~ as.character(trstplc))
) %>%
filter (!is.na(trstplc))france_clean <- france_data %>%
mutate (gndr = case_when (
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr %in% c(9) ~ NA_character_,
)) %>%
filter(!is.na(gndr))france_clean <- france_data %>%
mutate(
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"))
)# Regression model with 'gender' as the predictor
model_1 <- lm(trstplc ~ gndr, data = france_clean)
summary(model_1)##
## Call:
## lm(formula = trstplc ~ gndr, data = france_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0952 -1.0952 0.0629 1.9048 4.0629
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.09518 0.02220 274.554 < 2e-16 ***
## gndrMale -0.15806 0.03259 -4.849 1.25e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.24 on 18989 degrees of freedom
## Multiple R-squared: 0.001237, Adjusted R-squared: 0.001184
## F-statistic: 23.52 on 1 and 18989 DF, p-value: 1.249e-06# Model 2: Two predictors (gender & education) without interaction
model_2 <- lm(trstplc ~ gndr + educ.ba, data = france_clean)
summary(model_2)##
## Call:
## lm(formula = trstplc ~ gndr + educ.ba, data = france_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2877 -1.0414 0.1184 1.8721 4.1184
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.04136 0.02376 254.214 < 2e-16 ***
## gndrMale -0.15981 0.03256 -4.908 9.29e-07 ***
## educ.baBA or more 0.24634 0.03909 6.302 3.01e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.238 on 18988 degrees of freedom
## Multiple R-squared: 0.003321, Adjusted R-squared: 0.003216
## F-statistic: 31.64 on 2 and 18988 DF, p-value: 1.916e-14# Model 3: Interaction model involving both predictors
model_3 <- lm(trstplc ~ gndr + educ.ba + gndr * educ.ba, data = france_clean)
summary(model_3)##
## Call:
## lm(formula = trstplc ~ gndr + educ.ba + gndr * educ.ba, data = france_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2421 -1.0699 0.1517 1.8143 4.1517
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.06988 0.02508 242.033 < 2e-16 ***
## gndrMale -0.22157 0.03691 -6.003 1.97e-09 ***
## educ.baBA or more 0.11583 0.05366 2.159 0.030893 *
## gndrMale:educ.baBA or more 0.27794 0.07830 3.550 0.000387 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.237 on 18987 degrees of freedom
## Multiple R-squared: 0.003982, Adjusted R-squared: 0.003825
## F-statistic: 25.31 on 3 and 18987 DF, p-value: 2.483e-16# Create the models list
models_list <- list(model_1, model_2, model_3)
# Generate the table summary with a title
modelsummary(
models_list,
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
title = "Figure 1. Regression Model for Trust in Police")| (1) | (2) | (3) | |
|---|---|---|---|
| (Intercept) | 6.1 (0.0)*** | 6.0 (0.0)*** | 6.1 (0.0)*** |
| gndrMale | −0.2 (0.0)*** | −0.2 (0.0)*** | −0.2 (0.0)*** |
| educ.baBA or more | 0.2 (0.0)*** | 0.1 (0.1)* | |
| gndrMale × educ.baBA or more | 0.3 (0.1)*** | ||
| Num.Obs. | 18991 | 18991 | 18991 |
| R2 | 0.001 | 0.003 | 0.004 |
| R2 Adj. | 0.001 | 0.003 | 0.004 |
| AIC | 84530.2 | 84492.6 | 84482.0 |
| BIC | 84553.8 | 84524.0 | 84521.2 |
| Log.Lik. | −42262.117 | −42242.278 | −42235.979 |
| RMSE | 2.24 | 2.24 | 2.24 |
Coefficients Interpretation:
Intercept:
The intercept represents the average value of the response variable (authoritarian values scale) when all other predictor variables are zero. For Model 1, the intercept is 5.9 with a standard error of 0.1, indicating that the average authoritarian values score for the reference category of gender (omitted category) and education level is 5.9. The triple asterisks (***) denote that the coefficient is statistically significant at a very low p-value (under our alpha threshold of 0.05). gndr (Gender):
The coefficient for gender represents the change in the authoritarian values score associated with a one-unit change in gender (male vs. female), holding other variables constant. In Model 1, the coefficient for gender is 0.2 with a standard error of 0.1. This indicates that, on average, males have authoritarian values scores that are 0.2 points higher than females. The triple asterisks (***) denote that the coefficient is statistically significant at a very low p-value (under our alpha threshold of 0.05). educ.baBA or more (Education Level - Bachelor’s Degree or more):
The coefficient for the “BA or more” education level represents the change in authoritarian values score associated with having a Bachelor’s degree or higher compared to not having a Bachelor’s degree, holding other variables constant. In Model 1, the coefficient for “BA or more” education level is 0.3 with a standard error of 0.1. This indicates that, on average, individuals with a Bachelor’s degree or higher have authoritarian values scores that are 0.3 points higher than those without a Bachelor’s degree. The double asterisks (**) denote that the coefficient is statistically significant at a low p-value (0.01 < p < 0.05). gndr × educ.baBA or more (Interaction between Gender and Education Level):
The coefficient for the interaction between gender and education level represents the additional change in authoritarian values score associated with the combined effect of gender and having a Bachelor’s degree or higher, holding other variables constant. In Model 1, the coefficient for the interaction term is -0.3 with a standard error of 0.2. This indicates that the interaction between gender and having a Bachelor’s degree or higher is marginally significant, as the p-value is slightly above the alpha threshold (0.05 < p < 0.10). Adjusted R-Squared (Adj. R²):
The adjusted R-squared value measures the proportion of variability in the authoritarian values scale that is explained by the predictors in the model, adjusted for the number of predictors. For Model 1, the adjusted R-squared is 0.001, indicating that the predictors (gender and education level) explain only a very small proportion (0.1%) of the variability in authoritarian values scores. This suggests that other factors not included in the model may play a more significant role in shaping authoritarian values.
interaction_plot <- effect("gndr*educ.ba", model_3, na.rm=TRUE)
plot(interaction_plot,
main="Interaction Effect of Gender and Education on Trust in Police",
xlab="Education",
ylab="Trust in Police")
interaction_plot##
## gndr*educ.ba effect
## educ.ba
## gndr No BA BA or more
## Female 6.069876 6.185701
## Male 5.848307 6.242073- OUR VISUAL
GGPLOT
france_plot <- france_data %>%
mutate(
trstplc= case_when(
trstplc%in% c(7,8,9) ~ NA_character_,
TRUE ~ as.character(trstplc))
) %>%
filter (!is.na(trstplc))france_plot <- france_plot %>%
mutate (gndr = case_when (
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
gndr %in% c(9) ~ NA_character_,
)) %>%
filter(!is.na(gndr))france_plot <- france_plot %>%
mutate(
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(france_plot$trstplc, france_plot$gndr) %>%
cprop() %>%
as.data.frame() %>%
filter(Var1 != "Total",
Var2 != "All") %>%
ggplot(aes(x=Var1, y=Freq, fill=Var2)) +
geom_col(position = "dodge") +
labs(title="Trust in Police by Gender",
y = "Percentage of Trust",
x = "no trust - to - complete trust",
fill = "gender")
table(france_plot$impenv, france_plot$educ.ba) %>%
cprop() %>%
as.data.frame() %>%
filter(Var1 != "Total",
Var2 != "All") %>%
ggplot(aes(x=Var1, y=Freq, fill=Var2)) +
geom_col(position = "dodge") +
labs(title="Trust in Police by Education",
y = "Percentage of Trust",
x = "no trust - to - complete trust",
fill = "Education")