FINAL ASSIGNMENT
Rylee Boccongelle
2024-04-09
# 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"Part C Descriptive Statistics of Variables
uk_data <- read_fst("uk_data.fst")# For our outcome variable we are doinf the variable happy
uk_happy <- uk_data %>%
filter(cntry == "GB") %>%
dplyr:: select(happy)
uk_happy$y <- uk_happy$happy
table(uk_happy$y)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88
## 113 112 240 502 670 1568 1541 3804 6099 3740 2557 7 26#removing the no answer observations
uk_happy$y[uk_happy$y %in% 88:77] <- NA
table(uk_happy$y)##
## 0 1 2 3 4 5 6 7 8 9 10
## 113 112 240 502 670 1568 1541 3804 6099 3740 2557# This is a descriptive statistics table for the outcome variable "happiness" using the cleaned data we did above
tablehappy <- datasummary_skim ((uk_happy$y), title = "Table 1. Descriptive Statistics of 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.5 | 1.9 | 0.0 | 8.0 | 10.0 |
# For our secondary explanatory variable we are doing lrscale variable which is how respondents fall on the political scale, as with the happiness variable we need to clean the data to ensure we remove the NA responses
table(uk_data$lrscale)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88
## 487 305 802 1728 1935 7778 1959 1837 1149 271 366 75 2287#removing the no answer observations
uk_data <- uk_data %>%
mutate(
lrscale = na_if(lrscale, 77) %>%
na_if(., 88)
)
table(uk_data$lrscale)##
## 0 1 2 3 4 5 6 7 8 9 10
## 487 305 802 1728 1935 7778 1959 1837 1149 271 366#Now that we have cleaned the data for the trust in politicians variable (trustplt) above we can now create a descriptive statistics table
tablepolitical <- datasummary_skim ((uk_data$lrscale), title = "Table 2. Descriptive Statistics of Political Scale", 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 | 11 | 5.0 | 1.9 | 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 or no BA. First we need to clean it, like we did for the happiness variable.
table(uk_data$edulvla)##
## 1 2 3 5 55 77 88 99
## 2506 1988 1046 3000 114 2 21 18# Cleaning data for variable "edulvla"; have to remove no response observations to not skew data
uk_data <- uk_data %>%
mutate(
edulvla = na_if(edulvla, 55) %>%
na_if(., 77) %>%
na_if(., 88) %>%
na_if(., 99),
)
table(uk_data$edulvla)##
## 1 2 3 5
## 2506 1988 1046 3000uk_data <- uk_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(55, 77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba)
)
table(uk_data$educ.ba)##
## BA or more No BA
## 6395 14584tableeducation <- datasummary_skim ((uk_data$edulvla), title = "Table 3. Descriptive Statistics of 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 | 5 | 59 | 2.9 | 1.7 | 1.0 | 2.0 | 5.0 |
#This is a table of the descriptive statistics for to compare all three variables of interesr. This table will include the 2 explanatory variables (happiness and trust in politicians) and it will include our outcome variable educational attain of BA or no BA.
uk_clean <- uk_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(55, 77, 88, 99), NA, edulvla),
)
uk_clean <- uk_clean %>% filter(!is.na(happy), !is.na(lrscale), !is.na(edulvla))
table1descriptive <- datasummary_skim((uk_clean %>% dplyr::select(happy, lrscale, edulvla)),title = "Table 1: Descriptive Statistics of the Variables Happiness, Trust in Politicians, and Education Level from the ESS dataset in the UK", 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.5 | 1.9 | 0.0 | 8.0 | 10.0 |
lrscale | 11 | 0 | 5.1 | 1.8 | 0.0 | 5.0 | 10.0 |
edulvla | 4 | 0 | 3.0 | 1.7 | 1.0 | 2.0 | 5.0 |
``
Null Distribution and test statistic
#recoding education outcome variable
uk_data <- uk_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(55, 77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba)
)
table(uk_data$educ.ba)##
## BA or more No BA
## 6021 14958#Clean happiness variable
uk_data <- uk_data %>%
mutate(
happiness = ifelse(happy %in% c(88, 77), NA, happy)
)
table(uk_data$happiness)##
## 0 1 2 3 4 5 6 7 8 9 10
## 113 112 240 502 670 1568 1541 3804 6099 3740 2557#recoding lrsce
uk_data <- uk_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, 99) ~ NA_character_, # Values 77, 88, 99 in 'lrscale' are set as NA
TRUE ~ NA_character_ # Any other values not explicitly matched above
),
)
table(uk_data$ID)##
## Left Right
## 3322 3623# the POLID variable is cleaned (see below) the NA represents moderate values (THIS IS FOR THE DESCRIPTIVE STATISTICS PORTION OF THE ASSIGNMENT IT WILL NOT LET ME KNIT UNLESS I PUT IT HERE IT IS NOT PART OF THE NULL DISTRIBUTION STEPS)
table1b <- datasummary_skim(uk_data %>% dplyr::select(ID, educ.ba), type = "categorical", title = "Table 2. Descriptive statistics for explanatory variables", output = "flextable")
table1b
|
| N | % |
|---|---|---|---|
ID | Left | 3322 | 15.8 |
Right | 3623 | 17.3 | |
NA | 14034 | 66.9 | |
educ.ba | BA or more | 6021 | 28.7 |
No BA | 14958 | 71.3 |
# First we are calculating the test statistic using happiness as our outcome variable and education of BA or No BA as our oexplanatory variable. Based on the chart shown in lecture, happiness is a numerical variable since it was measured on a scale of 1-10 and education is a categorical varaible based on BA or No BA, therefore a t test statistic would be appropriate.
test_stat <- uk_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
## 3.338891print(test_stat$stat)## t
## 3.338891#now we are simulating the null distribution
null_distribution <- uk_data %>%
specify(explanatory = educ.ba,
response = happy) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
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.# Now we are calculating the 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 0#Through this code, we are now going to visualize the null distribution and using the two-sided direction as this is appropriate for a t statistic.
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: happy (numeric)
## Explanatory: educ.ba (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.0701
## 2 2 -3.36
## 3 3 1.96
## 4 4 1.42
## 5 5 -0.437
## 6 6 -0.260
## 7 7 0.264
## 8 8 2.05
## 9 9 0.877
## 10 10 0.222
## # ℹ 990 more rows# Based on the visualization of the null distribution, it is evident that under frequentist assumptions that there is strong evidence to reject the null hypothesis. To come to this conclusion, it is firstly apparent that the red line which represents the test statistic is outside of the null distribution and outside of the confidence intreval.Additionally, the p value which is 0.002 which is very small. The common threshold for the P-value is 0.05, since our p-value is under this there is evidence to reject the null hypothesis and there's a less than 5 percet chance that the result occured if the null was true.Regression Modelling
# first we are going to ensure that the data is weighted
uk_data$weight <- uk_data$dweight * uk_data$pweight
survey_design <- svydesign(ids = ~1, data = uk_data, weights = ~weight)# First we will produce a SLR Model using education of BA or no BA as our explanatory variable
model1 <- lm(happy ~ educ.ba, data = uk_data, weights = weight)
modelsummary(
list(model1),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "Table 3: SLR Model - Education Level and Happiness",
statistic = NULL)| Â (1) | |
|---|---|
| (Intercept) | 7.8 (0.0)*** |
| educ.baNo BA | −0.2 (0.1)** |
| Num.Obs. | 20979 |
| R2 | 0.000 |
| R2 Adj. | 0.000 |
| AIC | 113295.6 |
| BIC | 113319.4 |
| Log.Lik. | −56644.785 |
| RMSE | 3.65 |
# From this chart outputted it is evidence that the coefficient from the variable educ.ba for no BA of -0.2 means that on average those from the UK who do not have a BA score 0.2 points (on a scale from 0-10) lower on the happiness scale than those from the UK who have a BA. Additionally, there are 3 asterisks stars (*) which tells us that there is statistical significance between the variables. However, all though the variable is statistically, signifigant, the adjusted R squared value of 0.2 percent indicates that the educ.ba variable only explains 0.2 percent of the variance for happiness levels in the UK after accounting for other predictors. Therefore this tells us that there are other factors outside of education that point to shaping happiness levels in the UK.# Now we will make a MLR model with both explanatory variables - education and trust in politicians
model2 <- lm(happy ~ educ.ba + ID, data = uk_data, weights = weight)
modelsummary(
list(model2),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "MLR Model - Education Level, Political Identification, and Happiness",
statistic = NULL)| Â (1) | |
|---|---|
| (Intercept) | 7.5 (0.1)*** |
| educ.baNo BA | −0.2 (0.1)** |
| IDRight | 0.5 (0.1)*** |
| Num.Obs. | 6945 |
| R2 | 0.009 |
| R2 Adj. | 0.008 |
| AIC | 33373.2 |
| BIC | 33400.6 |
| Log.Lik. | −16682.591 |
| RMSE | 2.71 |
# From this chart outputted, it is evident that those from the UK who fall on the right of th political spectrum, score on average, 0.5 points higher on the happiness scale than those who fall on the ideological left. Additionally of important note, there are 3 asterisks indicating that this relationship between happiness and political identification is statistically significant. # Now we will make an interaction model involving both explanatory variablesmodel3 <- lm(happy ~ educ.ba + ID + educ.ba*ID, data = uk_data, weights = weight)
modelsummary(
list(model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "MLR Interaction Model: Education Level, Political Identification, and Happiness",
statistic = NULL)| Â (1) | |
|---|---|
| (Intercept) | 7.5 (0.1)*** |
| educ.baNo BA | −0.2 (0.1)+ |
| IDRight | 0.5 (0.1)*** |
| educ.baNo BA × IDRight | 0.0 (0.1) |
| Num.Obs. | 6945 |
| R2 | 0.009 |
| R2 Adj. | 0.008 |
| AIC | 33375.1 |
| BIC | 33409.3 |
| Log.Lik. | −16682.545 |
| RMSE | 2.71 |
# From this table outputted, it is evident that those without a BA and fall on the right of the political ideological spectrum fall on average, 0.1 points higher on the happiness scale.#Finally, we will create a model summary table of all three models we have generarted above (SLR, MLR, AND MLR Interaction effect).
modelsummary(
list(model1, model2, model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "Table 6: Model Summary Table for Education Level, Political Identification & Happiness",
statistic = NULL)| Â (1) | Â Â (2) | Â Â (3) | |
|---|---|---|---|
| (Intercept) | 7.8 (0.0)*** | 7.5 (0.1)*** | 7.5 (0.1)*** |
| educ.baNo BA | −0.2 (0.1)** | −0.2 (0.1)** | −0.2 (0.1)+ |
| IDRight | 0.5 (0.1)*** | 0.5 (0.1)*** | |
| educ.baNo BA × IDRight | 0.0 (0.1) | ||
| Num.Obs. | 20979 | 6945 | 6945 |
| R2 | 0.000 | 0.009 | 0.009 |
| R2 Adj. | 0.000 | 0.008 | 0.008 |
| AIC | 113295.6 | 33373.2 | 33375.1 |
| BIC | 113319.4 | 33400.6 | 33409.3 |
| Log.Lik. | −56644.785 | −16682.591 | −16682.545 |
| RMSE | 3.65 | 2.71 | 2.71 |
Visual Output of Choice
interaction_plot <- effect("educ.ba*ID", model3, na.rm=TRUE)
plot(interaction_plot,
main="Interaction Effect",
xlab="ID",
ylab="Happiness Scale")
model4 <- lm(happy ~ + educ.ba + lrscale + educ.ba*lrscale, data = uk_data, weights = weight)
interact_plot(model4, pred = lrscale, modx = educ.ba, jnplot = TRUE)
cat_plot(model3, pred = ID, modx = educ.ba, jnplot = TRUE)