Research Project
Jordan MacLean
2024-04-11
Outcome Variable (Dependent Variable): Trust in Parliament ‘trstprl’
Predictor/Explanatory Variable (Independent Variable): Educational Attainment ‘edu.ba’ ‘edulvla’ ‘edulvlb’
Categorical Variable (Independent Variable): Feeling of Household Income ‘hincfel’
RQ: How does trust in France’s parliament depend on educational attainment
# Solves any issues concerning 'vector memory' errors and clears the environment
rm(list=ls()); gc() ## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 527816 28.2 1173344 62.7 NA 669402 35.8
## Vcells 972181 7.5 8388608 64.0 16384 1851829 14.2# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "performance", "flextable", "broom", "scales", "ggeffects", "marginaleffects", "haven", "lme4", "interactions") # 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
##
##
##
## Attaching package: 'interactions'
##
##
## The following object is masked from 'package:ggeffects':
##
## johnson_neyman## [[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"
##
## [[13]]
## [1] "haven" "marginaleffects" "ggeffects" "scales"
## [5] "broom" "flextable" "performance" "survey"
## [9] "survival" "Matrix" "grid" "effects"
## [13] "carData" "modelsummary" "fst" "infer"
## [17] "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble"
## [25] "ggplot2" "tidyverse" "stats" "graphics"
## [29] "grDevices" "utils" "datasets" "methods"
## [33] "base"
##
## [[14]]
## [1] "lme4" "haven" "marginaleffects" "ggeffects"
## [5] "scales" "broom" "flextable" "performance"
## [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] "interactions" "lme4" "haven" "marginaleffects"
## [5] "ggeffects" "scales" "broom" "flextable"
## [9] "performance" "survey" "survival" "Matrix"
## [13] "grid" "effects" "carData" "modelsummary"
## [17] "fst" "infer" "lubridate" "forcats"
## [21] "stringr" "dplyr" "purrr" "readr"
## [25] "tidyr" "tibble" "ggplot2" "tidyverse"
## [29] "stats" "graphics" "grDevices" "utils"
## [33] "datasets" "methods" "base"# Reading in the France data set
france_data <- read_fst("france_data.fst")# Here we code specifically for our research questions outcome variable, Trust in Parliament
france_data <- france_data %>%
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 %>% filter(!is.na(trstprl))# Here, we assign years based on ESS round, specifically to the varaible france_data. Our research thus contains data from the years 2002 to 2020
france_data$year <- NA
replacements <- c(2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020)
for (i in 1:10) {
france_data$year[france_data$essround == i] <- replacements [i]
}# Coding a ggplot solely on France Trust in Parliament from 2002-2020
trstprl_avg <- france_data %>%
group_by(year) %>%
summarize(trstprl_avg = mean(trstprl, na.rm = TRUE))
trstprl_plot_FR <- ggplot(trstprl_avg, aes(x = year, y = trstprl_avg)) +
geom_point(aes(size = trstprl_avg, color = trstprl_avg), alpha = 0.7) +
geom_line(aes(group = 1), color = "blue", linetype = "dashed") +
labs(
x = "Survey Year",
y = "Trust In Parliament") +
theme_minimal() +
theme(legend.position = "none")
# Displaying the ggplot
print(trstprl_plot_FR)
# Coding for the variables
france_data <- france_data %>%
mutate(
feelincome = case_when(
hincfel == 1 ~ "comfortable",
hincfel %in% 2:4 ~ "difficult",
TRUE ~ NA_character_
),
edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
),
# Ensure educ.ba is treated as a factor with specified levels
educ.ba = factor(educ.ba, levels = c("No BA", "BA or more"))
)
# Now correct the assignment to handle values greater than 10 for 'trstprl'
france_data$trstprl <- ifelse(france_data$trstprl > 10, NA, france_data$trstprl)
# Assuming the operation to create and rescale 'trust' is correct and intentional:
france_data$trust <- scales::rescale(france_data$trstprl, na.rm = TRUE, to = c(0, 100))# Checking to see if there is data within the 'trust' variable
table(france_data$trust)##
## 0 10 20 30 40 50 60 70 80 90 100
## 1826 767 1772 2363 2308 4468 1983 1716 1050 200 221francedf <- france_data %>%
dplyr::select(trust, feelincome, educ.ba, weight)
# Checking the first few rows to verify the selection
francedf <- na.omit(francedf)
head(francedf)## trust feelincome educ.ba weight
## 3310 20 difficult No BA 1.233645
## 3312 30 difficult No BA 2.151479
## 3313 30 difficult BA or more 1.245936
## 3314 70 comfortable BA or more 2.930847
## 3315 20 comfortable No BA 3.060282
## 3316 50 comfortable BA or more 1.372110# Creating models for our regression model table
model1 <- lm(trust ~ educ.ba, data = francedf, weights = weight)
model2 <- lm(trust ~ feelincome, data = francedf, weights = weight)
model3 <- lm(trust ~ educ.ba*feelincome, data = francedf, weight = weight)
# Using modelsummary to display the models
modelsummary(
list(model1, model2, model3),
title = "Regression Model Analysis of Trust in Parliament: The Role of Educational Attainment and Feelings of Household Income",
fmt = 1,
estimate = "{estimate} ({std.error}){stars}",
statistic = NULL
)| Â (1) | Â Â (2) | Â Â (3) | |
|---|---|---|---|
| (Intercept) | 39.8 (0.2)*** | 48.2 (0.3)*** | 45.1 (0.4)*** |
| educ.baBA or more | 11.7 (0.5)*** | 8.8 (0.7)*** | |
| feelincomedifficult | −9.0 (0.4)*** | −7.4 (0.5)*** | |
| educ.baBA or more × feelincomedifficult | 1.9 (0.9)* | ||
| Num.Obs. | 15375 | 15375 | 15375 |
| R2 | 0.041 | 0.033 | 0.060 |
| R2 Adj. | 0.041 | 0.033 | 0.059 |
| AIC | 141627.5 | 141757.8 | 141330.0 |
| BIC | 141650.4 | 141780.7 | 141368.2 |
| Log.Lik. | −70810.744 | −70875.912 | −70660.015 |
| RMSE | 23.01 | 23.12 | 22.81 |
# Creating a numeric table
# Correcting the code to convert 'feelincome' from categorical to numeric so that we can include it in the summary table (did not work as a categorical)
# We solved this issue by creating a numeric variable 'feelincome_numeric' with the 'francedf'
francedf <- francedf %>%
mutate(
feelincome_numeric = case_when(
feelincome == "comfortable" ~ 1,
feelincome == "difficult" ~ 2,
TRUE ~ NA_real_ # This line of code puts NA values as missing data
)
)
# Instead of making feelings = feelincome, from our new numeric varaible it is Feelings = feelincome_numeric
francet1 <- francedf %>%
mutate(
Trust = as.numeric(trust), # Ensure Trust is numeric
Feelings = feelincome_numeric, # Use the newly created numeric 'Feelings'
Education = as.numeric(educ.ba) # Ensure Education is numeric, if this is not already the case
)
# Then we generate a summary table of the numeric variables
table1b <- datasummary_skim(francet1 %>% select(Trust, Feelings, Education),
type = "numeric",
title = "Table 1. Descriptive Statistics for Main Variables",
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.# Now we output the table
table1b
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
Trust | 11 | 0 | 42.2 | 23.5 | 0.0 | 50.0 | 100.0 |
Feelings | 2 | 0 | 1.7 | 0.5 | 1.0 | 2.0 | 2.0 |
Education | 2 | 0 | 1.2 | 0.4 | 1.0 | 1.0 | 2.0 |
# Code for printing the table as a PDF
set_flextable_defaults(fonts_ignore=TRUE)
flextable::save_as_docx(table1b, path = "Table1.docx", # change name to whatever you would like
width = 7.0, height = 7.0) # play around with dimensions
# For PDFHYPOTHESIS TESTING WITH INFER
STEP 1 CALCULATE THE TEST STATISTIC OF YOUR SAMPLE
# Calculating the test statistic
test_stat <- france_data %>%
specify(explanatory = educ.ba,
response = trust) %>%
hypothesize(null = "independence") %>%
calculate(stat = "t")## Warning: Removed 364 rows containing missing values.## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "No BA" - "BA or more", or divided in the order "No BA" / "BA or more"
## for ratio-based statistics. To specify this order yourself, supply `order =
## c("No BA", "BA or more")` to the calculate() function.# Printing the test statistic
print(test_stat$stat)## t
## -28.14463# the t-statistic measures the difference btween the observed sample mean and the population mean (or another sample mean, in the case of two-sample tests) in units of standard error
# t statistic is -26.21425STEP 2 SIMULATE THE NULL DISTRIBUTION
# Simulating the null distribution
null_distribution <- france_data %>%
specify(response = trust, explanatory = educ.ba) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "t")## Warning: Removed 364 rows containing missing values.## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "No BA" - "BA or more", or divided in the order "No BA" / "BA or more"
## for ratio-based statistics. To specify this order yourself, supply `order =
## c("No BA", "BA or more")` to the calculate() function.STEP 3 CALCULATE THE P-VALUE OF YOUR SAMPLE
# Calculating the p-value of our sample
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.# Printing the p value
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0# Really, it's essentially a p-value of 0, not 0 - but, still, under the assumptions of frequentist hypothesis testing, this means there is strong evidence to reject nullSTEP 4 VISUALIZE
# Visualizing the null distribution along the observed test-statistic (red line), indicating clear evidence to reject the null under assumptions of frequentist hypothesis testing
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two_sided")
null_distribution## Response: trust (numeric)
## Explanatory: educ.ba (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.286
## 2 2 0.682
## 3 3 -0.778
## 4 4 1.03
## 5 5 0.939
## 6 6 0.568
## 7 7 0.247
## 8 8 0.189
## 9 9 0.0303
## 10 10 -1.02
## # ℹ 990 more rows# Visualizing with confidence intervals
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)
# Gaining a summary of model3
summary(model3)##
## Call:
## lm(formula = trust ~ educ.ba * feelincome, data = francedf, weights = weight)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -160.382 -25.177 3.797 23.075 162.419
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.1228 0.3914 115.272 <2e-16 ***
## educ.baBA or more 8.7856 0.6656 13.199 <2e-16 ***
## feelincomedifficult -7.3950 0.4601 -16.074 <2e-16 ***
## educ.baBA or more:feelincomedifficult 1.8758 0.9350 2.006 0.0449 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37.31 on 15371 degrees of freedom
## Multiple R-squared: 0.05958, Adjusted R-squared: 0.0594
## F-statistic: 324.6 on 3 and 15371 DF, p-value: < 2.2e-16# Interaction Plot
model3 <- lm(trust ~ educ.ba*feelincome, data = francedf, weights = francedf$weight)
interaction_plot <- effect("educ.ba*feelincome", model3)
# Ploting the interaction effect
plot(interaction_plot,
main="Interaction Effect",
xlab="Educational Attainment",
ylab="Trust in Parliament") # Making sure the labels match the variables
#The lines are not exactly parallel as one line is more steep than the other though appearing very minisculeinteraction_plot##
## educ.ba*feelincome effect
## feelincome
## educ.ba comfortable difficult
## No BA 45.1228 37.72778
## BA or more 53.9084 48.38916# Gaining a summary of model3 once again
summary(model3)##
## Call:
## lm(formula = trust ~ educ.ba * feelincome, data = francedf, weights = francedf$weight)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -160.382 -25.177 3.797 23.075 162.419
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.1228 0.3914 115.272 <2e-16 ***
## educ.baBA or more 8.7856 0.6656 13.199 <2e-16 ***
## feelincomedifficult -7.3950 0.4601 -16.074 <2e-16 ***
## educ.baBA or more:feelincomedifficult 1.8758 0.9350 2.006 0.0449 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37.31 on 15371 degrees of freedom
## Multiple R-squared: 0.05958, Adjusted R-squared: 0.0594
## F-statistic: 324.6 on 3 and 15371 DF, p-value: < 2.2e-16# Reintroducing the model3 so that it is recognized
model3 <- glm(trust ~ educ.ba*feelincome, data = francedf, weight = weight) # both predictors + interaction# Graphing the ggplot
ggplot(francedf, aes(x = educ.ba, y = trust, fill = factor(feelincome))) +
geom_boxplot() +
labs(title = "Trust in Parliament by Educational Attainment and Feelings of Household Income",
x = "Educational Attainment",
y = "Trust in Parliament") +
scale_fill_viridis_d(name = "Household Income Feeling", option = "A") +
theme_minimal()