Research Report 1
Luna Marin, Shamiso Choga, and Favour Ayoola-Ola
2024-04-11
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "MASS", "aod", "interactions", "kableExtra", "flextable", "scales") # 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: '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"Section 1: Calculating Descriptive Statistics for the Dataset - Analyzing Variables
france_data <- read_fst("france_data.fst")#For our outcome variable, we are focusing on happiness.
fr_happy <- france_data %>%
filter(cntry == "FR") %>%
dplyr:: select(happy)
fr_happy$y <- fr_happy$happy
table(fr_happy$y)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88
## 90 57 185 346 514 2349 1869 4003 5591 2386 1622 14 12#Excluding observations with no answers.
fr_happy$y[fr_happy$y %in% 88:77] <- NA
table(fr_happy$y)##
## 0 1 2 3 4 5 6 7 8 9 10
## 90 57 185 346 514 2349 1869 4003 5591 2386 1622#This table provides descriptive statistics for the outcome variable "happiness" using the cleaned data from earlier.
tablehappy <- datasummary_skim ((fr_happy$y), title = "Table 1. Descriptive Statistics of the Happiness Variable", 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.tablehappy
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
data | 12 | 0 | 7.2 | 1.8 | 0.0 | 8.0 | 10.0 |
#Our secondary explanatory variable is the political scale (lrscale), which measures how respondents position themselves on the political spectrum. Similar to the happiness variable, we need to clean this data by removing NA responses to ensure accuracy in our analysis.
table(france_data$lrscale)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88
## 945 499 1269 1996 1838 5274 1594 1749 1458 409 686 381 940#Excluding responses with no answer observations.
fr_data <- france_data %>%
mutate(
lrscale = na_if(lrscale, 77) %>%
na_if(., 88)
)
table(fr_data$lrscale)##
## 0 1 2 3 4 5 6 7 8 9 10
## 945 499 1269 1996 1838 5274 1594 1749 1458 409 686#With the data for the trust in politicians variable now cleaned, we can proceed to create a descriptive statistics table.
tablepolitical <- datasummary_skim ((fr_data$lrscale), title = "Table 2. Descriptive Statistics for the Political Scale Variable", 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.tablepolitical
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
data | 12 | 7 | 4.9 | 2.4 | 0.0 | 5.0 | 10.0 |
#For the main explanatory variable, we are using the education variable "edulvla," which measures whether participants have a BA degree or no BA. As with the happiness variable, we need to clean this data first by removing missing or invalid responses to ensure a clear analysis.
table(fr_data$edulvla)##
## 1 2 3 4 5 77 88
## 1449 1018 2879 13 2003 2 4#Clean the 'edulvla' variable by replacing values 77 and 88 with NA
france_data <- france_data %>%
mutate(
edulvla = na_if(edulvla, 77) %>%
na_if(., 88)
)
#Display the frequency table of the cleaned 'edulvla' variable
table(france_data$edulvla)##
## 1 2 3 4 5
## 1449 1018 2879 13 2003#With the data now cleaned, we can proceed to create a descriptive statistics table.
tableeducation <- datasummary_skim ((fr_data$edulvla), title = "Table 3. Descriptive Statistics for the Highest Level of Education Variable", 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.tableeducation
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
data | 8 | 61 | 3.1 | 2.7 | 1.0 | 3.0 | 88.0 |
#This table presents the descriptive statistics for comparing all three variables of interest: the two explanatory variables (happiness and trust in politicians) and the outcome variable (educational attainment of BA or no BA).
fr_clean <- fr_data %>%
mutate(
happy = ifelse(happy %in% c(77, 88), NA, happy),
lrsclae = ifelse(trstplt %in% c(77, 88), NA, lrscale),
edulvla = ifelse(edulvla %in% c( 77, 88), NA, edulvla),
)
fr_clean <- fr_clean %>% filter(!is.na(happy), !is.na(lrscale), !is.na(edulvla))
table1descriptive <- datasummary_skim((fr_clean %>% dplyr::select(happy, lrscale, edulvla)),title = "Table 4: Descriptive Statistics of the Variables Happiness, Trust in Politicians, and Education Level from the ESS dataset in France", 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.table1descriptive
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
happy | 11 | 0 | 7.2 | 1.8 | 0.0 | 8.0 | 10.0 |
lrscale | 11 | 0 | 4.8 | 2.3 | 0.0 | 5.0 | 10.0 |
edulvla | 5 | 0 | 3.0 | 1.4 | 1.0 | 3.0 | 5.0 |
Section 2: Using the infer package, we can generate the null distribution (along with a confidence interval) and calculate the test statistic.
#Recoding the "Highest Level of Education" variable (edulvla) into the categories "BA" or "no BA".
fr_data <- fr_data %>%
mutate (
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
mutate(
edulvla = ifelse(edulvla %in% c(77, 88), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(7777, 8888), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba)
)
table(fr_data$educ.ba)##
## BA or more No BA
## 4252 14786#Confirming that the "happy" variable is cleaned. No recoding is necessary, as the numerical scale from 0 to 10 effectively describes the happiness level.
fr_data <- fr_data %>%
mutate(
happiness = ifelse(happy %in% c(88, 77), NA, happy)
)
table(fr_data$happiness)##
## 0 1 2 3 4 5 6 7 8 9 10
## 90 57 185 346 514 2349 1869 4003 5591 2386 1622#Recoding the "lrscale" variable to categorize respondents' political leanings on the left-right scale.
fr_data <- fr_data %>%
mutate(
# Recoding political orientation
ID = case_when(
lrscale %in% 0:3 ~ "Left", # Values 0 to 3 in 'lrscale' are categorized as "Left"
lrscale %in% 7:10 ~ "Right", # Values 7 to 10 in 'lrscale' are categorized as "Right"
lrscale %in% c(77, 88) ~ NA_character_, # Values 77, 88 in 'lrscale' are set as NA
TRUE ~ NA_character_ # Any other values not explicitly matched above
),
)
table(fr_data$ID)##
## Left Right
## 4709 4302# We begin by calculating the test statistic using happiness (measured on a scale from 1 to 10) as the outcome variable and education (categorized as BA or No BA) as the explanatory variable. A two-sample t-test is appropriate given the categorical and numeric nature of the variables. This test is chosen because education level is our main predictor and the response variable is happiness.
test_stat <- fr_data %>%
specify(explanatory = educ.ba,
response = happy) %>%
hypothesize(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 "BA or more" - "No BA", or divided in the order "BA or more" / "No BA"
## for ratio-based statistics. To specify this order yourself, supply `order =
## c("BA or more", "No BA")` to the calculate() function.print(test_stat$stat)## t
## 7.749007#Now, we simulate the null distribution to evaluate the relationship between happiness and education level.
null_distribution <- fr_data %>%
specify(explanatory = educ.ba,
response = happiness) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "t")## Warning: Removed 26 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 "BA or more" - "No BA", or divided in the order "BA or more" / "No BA"
## for ratio-based statistics. To specify this order yourself, supply `order =
## c("BA or more", "No BA")` to the calculate() function.#Now, we calculate the p-value to assess the statistical significance of the observed relationship between happiness and education level.
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 0#This code visualizes the null distribution using a two-sided direction, which is appropriate for a t-test.
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)
null_distribution## Response: happiness (numeric)
## Explanatory: educ.ba (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.147
## 2 2 -0.371
## 3 3 -0.211
## 4 4 -0.595
## 5 5 -0.953
## 6 6 0.391
## 7 7 0.127
## 8 8 0.236
## 9 9 -1.16
## 10 10 2.31
## # ℹ 990 more rowsInterpretation: The null distribution visualization shows compelling evidence against the null hypothesis under frequentist assumptions. This is clear because the red line, which denotes the test statistic, falls outside both the null distribution and its confidence interval. This separation indicates a strong reason to reject the null hypothesis.
Section 3: Analyzing the data using regression models, including single, multivariate, and interaction effects.
# To begin, we will apply appropriate weights to the data to ensure the analysis accurately represents the sample. This weighting will help account for any biases in the dataset and make the results more representative of the population.
fr_data$weight <- fr_data$dweight * fr_data$pweight
survey_design <- svydesign(ids = ~1, data = fr_data, weights = ~weight)# First, we will create a simple linear regression (SLR) model with the education level (BA or no BA) as the explanatory variable. This model will help us understand the relationship between the highest level of education and our outcome variable, such as happiness.
model1 <- lm(happiness ~ educ.ba, data = fr_data, weights = weight)
modelsummary(
list(model1),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "Table 5: Simple Linear Regression Model (SLR) - Relationship between Education Level and Happiness",
statistic = NULL)| (1) | |
|---|---|
| (Intercept) | 7.6 (0.0)*** |
| educ.baNo BA | −0.4 (0.0)*** |
| Num.Obs. | 19012 |
| R2 | 0.007 |
| R2 Adj. | 0.007 |
| AIC | 77124.5 |
| BIC | 77148.1 |
| Log.Lik. | −38559.258 |
| RMSE | 1.79 |
Interpretation: The coefficient for the “No BA” education level in France is -0.4 with a standard error of 0.0, denoted by *** to indicate a high level of statistical significance. This negative coefficient suggests that individuals in France without a BA degree experience a decrease of 0.4 points in happiness compared to those with a BA or higher degree. However, the model’s R-squared (R2) value of 0.007 indicates that the model explains only 0.7% of the variance in happiness. This is a very low proportion, which suggests that the model has limited ability to explain the variation in happiness among individuals in France. In summary, the model points to a slight decrease in happiness for individuals without a BA degree compared to those with a BA or higher degree in France. Nonetheless, the model’s low explanatory power implies that there may be other, unaccounted-for factors that influence happiness more significantly in this context.
# In this section, we will create a Multiple Linear Regression (MLR) model to examine the relationship between the outcome variable happiness and two explanatory variables: education level (BA or no BA) and trust in politicians.
model2 <- lm(happiness ~ educ.ba + ID, data = fr_data, weights = weight)
modelsummary(
list(model2),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "Table 6: Multiple Linear Regression (MLR) Model - Relationship between Education Level, Political Identification, and Happiness",
statistic = NULL)| (1) | |
|---|---|
| (Intercept) | 7.6 (0.0)*** |
| educ.baNo BA | −0.4 (0.0)*** |
| IDRight | 0.2 (0.0)*** |
| Num.Obs. | 9000 |
| R2 | 0.011 |
| R2 Adj. | 0.011 |
| AIC | 36209.9 |
| BIC | 36238.3 |
| Log.Lik. | −18100.942 |
| RMSE | 1.76 |
Interpretation: Overall, the model suggests that in France, individuals who identify as right-leaning politically tend to be slightly happier compared to centrist individuals, scoring on average 0.2 points higher on the happiness scale. Additionally, individuals without a BA degree tend to be slightly less happy compared to those with a BA degree, scoring on average 0.4 points lower on the happiness scale. Both of these conclusions are based on the statistically significant coefficients: the coefficient for “No BA” education level is -0.4, and the coefficient for right-leaning political identification is 0.2.
#We will now develop an interaction model that includes both explanatory variables: education level (BA or no BA) and political identification (on a left-right scale). This model will assess how the combination of these two factors affects the dependent variable, happiness. By considering the interaction between education and political identification, we can understand how the relationship between education and happiness varies depending on political leanings.
model3 <- lm(happiness ~ educ.ba + ID + educ.ba*ID, data = fr_data, weights = weight)
modelsummary(
list(model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "Table 7: MLR Interaction Model: Relationship between Education Level, Political Identification, and Happiness",
statistic = NULL)| (1) | |
|---|---|
| (Intercept) | 7.6 (0.0)*** |
| educ.baNo BA | −0.4 (0.1)*** |
| IDRight | 0.1 (0.1) |
| educ.baNo BA × IDRight | 0.1 (0.1) |
| Num.Obs. | 9000 |
| R2 | 0.012 |
| R2 Adj. | 0.011 |
| AIC | 36209.9 |
| BIC | 36245.5 |
| Log.Lik. | −18099.973 |
| RMSE | 1.76 |
Interpretation: From the model results, it is evident that individuals in France who do not have a BA degree and align themselves with the right of the political spectrum tend to have a 0.1-point increase in happiness on average. This interaction effect between education level and political identification in France, represented by the coefficient of 0.1, suggests that the combined effect of having no BA and identifying as right-leaning politically leads to a slight increase in happiness.
# Finally, we will compile a summary table that presents an overview of the three models we have created: the simple linear regression (SLR), multiple linear regression (MLR), and MLR with interaction effects.
modelsummary(
list(model1, model2, model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "Table 8: Linear Regression Models: SLR, MLR, & MLR Interaction for Education Level, Political Identification, & Happiness",
statistic = NULL)| (1) | (2) | (3) | |
|---|---|---|---|
| (Intercept) | 7.6 (0.0)*** | 7.6 (0.0)*** | 7.6 (0.0)*** |
| educ.baNo BA | −0.4 (0.0)*** | −0.4 (0.0)*** | −0.4 (0.1)*** |
| IDRight | 0.2 (0.0)*** | 0.1 (0.1) | |
| educ.baNo BA × IDRight | 0.1 (0.1) | ||
| Num.Obs. | 19012 | 9000 | 9000 |
| R2 | 0.007 | 0.011 | 0.012 |
| R2 Adj. | 0.007 | 0.011 | 0.011 |
| AIC | 77124.5 | 36209.9 | 36209.9 |
| BIC | 77148.1 | 36238.3 | 36245.5 |
| Log.Lik. | −38559.258 | −18100.942 | −18099.973 |
| RMSE | 1.79 | 1.76 | 1.76 |
Section 4: Presenting a Third Visualization.
# In this section, we have included an interaction plot to examine the effect of the interaction between educational level (BA or no BA) and political identification (ID) on happiness.
interaction_plot <- effect("educ.ba*ID", model3, na.rm=TRUE)
plot(interaction_plot,
main="Interaction Effect",
xlab="ID",
ylab="Happiness Scale")
interaction_plot##
## educ.ba*ID effect
## ID
## educ.ba Left Right
## BA or more 7.605636 7.70281
## No BA 7.168581 7.38515