Research Report All Coding FINAL
Gurjyot Kaur
2024-04-12
# List of packages
packages <- c("tidyverse", "modelsummary", "flextable", "infer",
"fst", "viridis", "knitr", "kableExtra", "rmarkdown", "ggridges", "questionr")
# 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
##
## Attaching package: 'flextable'
##
##
## The following object is masked from 'package:purrr':
##
## compose## Warning: package 'infer' was built under R version 4.3.3## Loading required package: viridisLite
##
## Attaching package: 'kableExtra'
##
## The following objects are masked from 'package:flextable':
##
## as_image, footnote
##
## The following object is masked from 'package:dplyr':
##
## group_rows## [[1]]
## [1] "lubridate" "forcats" "stringr" "dplyr" "purrr" "readr"
## [7] "tidyr" "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "modelsummary" "lubridate" "forcats" "stringr" "dplyr"
## [6] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [11] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [16] "datasets" "methods" "base"
##
## [[3]]
## [1] "flextable" "modelsummary" "lubridate" "forcats" "stringr"
## [6] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [11] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [16] "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "infer" "flextable" "modelsummary" "lubridate" "forcats"
## [6] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [11] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "fst" "infer" "flextable" "modelsummary" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[6]]
## [1] "viridis" "viridisLite" "fst" "infer" "flextable"
## [6] "modelsummary" "lubridate" "forcats" "stringr" "dplyr"
## [11] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [16] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [21] "datasets" "methods" "base"
##
## [[7]]
## [1] "knitr" "viridis" "viridisLite" "fst" "infer"
## [6] "flextable" "modelsummary" "lubridate" "forcats" "stringr"
## [11] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [16] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [21] "utils" "datasets" "methods" "base"
##
## [[8]]
## [1] "kableExtra" "knitr" "viridis" "viridisLite" "fst"
## [6] "infer" "flextable" "modelsummary" "lubridate" "forcats"
## [11] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [16] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [21] "grDevices" "utils" "datasets" "methods" "base"
##
## [[9]]
## [1] "rmarkdown" "kableExtra" "knitr" "viridis" "viridisLite"
## [6] "fst" "infer" "flextable" "modelsummary" "lubridate"
## [11] "forcats" "stringr" "dplyr" "purrr" "readr"
## [16] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [21] "graphics" "grDevices" "utils" "datasets" "methods"
## [26] "base"
##
## [[10]]
## [1] "ggridges" "rmarkdown" "kableExtra" "knitr" "viridis"
## [6] "viridisLite" "fst" "infer" "flextable" "modelsummary"
## [11] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [16] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"
##
## [[11]]
## [1] "questionr" "ggridges" "rmarkdown" "kableExtra" "knitr"
## [6] "viridis" "viridisLite" "fst" "infer" "flextable"
## [11] "modelsummary" "lubridate" "forcats" "stringr" "dplyr"
## [16] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [21] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [26] "datasets" "methods" "base"# Load Dataset:
sweden_data <- read.fst("sweden_data.fst")DESCRIPTIVE STATISTICS TABLE
Looking at the variables before cleaning:
table(sweden_data$stflife) # we do not need 5 (moderate satisfaction), 77 (refusal), 88 (don't know), or 99 (no answer)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88 99
## 67 50 155 282 375 1047 1034 2939 5630 3924 2679 8 22 4table(sweden_data$hincfel) # we can see 7 (refusal), 8 (don't know), and 9 (no answer) need to be omitted##
## 1 2 3 4 7 8 9
## 10688 5777 1264 330 36 91 30table(sweden_data$edulvla) # 55 (other)- it is on the ess variable description, not in this sample but just in case, 77 (refusal), 88 (don't know), and 99 (no answer) need to be omitted##
## 1 2 3 5 77 88 99
## 1132 1336 3356 1855 4 9 12Cleaning the variables:
sweden_clean <- sweden_data %>%
mutate(
stflife = ifelse(stflife %in% c(5, 77, 88, 99), NA, stflife), # For 'stflife,' set values 5, 77, 88, or 99 to NA.
hincfel = ifelse(hincfel %in% c(7, 8, 9), NA, hincfel), # For 'hincfel,' set values 7, 8, and 9 to NA.
edulvla = ifelse(edulvla %in% c(55, 77, 88, 99), NA, edulvla), # For 'edulvla,' set values 55, 77, 88, and 99 to NA.
)Removing NAs completely:
sweden_clean <- sweden_clean %>% filter(!is.na(stflife), !is.na(hincfel), !is.na(edulvla))Making sure the variables are cleaned:
table(sweden_clean$stflife) ##
## 0 1 2 3 4 6 7 8 9 10
## 28 21 59 107 155 438 1220 2362 1698 1116table(sweden_clean$hincfel)##
## 1 2 3 4
## 4156 2466 481 101table(sweden_clean$edulvla)##
## 1 2 3 5
## 1045 1225 3165 1769unique(sweden_clean$edulvla)## [1] 1 5 3 2unique(sweden_clean$hincfel)## [1] 2 1 3 4unique(sweden_clean$stflife)## [1] 8 10 9 6 7 4 3 2 1 0# Double-checking using unique to see any NAsTable summary:
table1b <- datasummary_skim(sweden_clean %>% dplyr::select(stflife, hincfel, edulvla), 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.table1b
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
stflife | 10 | 0 | 8.0 | 1.6 | 0.0 | 8.0 | 10.0 |
hincfel | 4 | 0 | 1.5 | 0.7 | 1.0 | 1.0 | 4.0 |
edulvla | 4 | 0 | 3.0 | 1.3 | 1.0 | 3.0 | 5.0 |
Not cleaned table:
table_not_cleaned <- datasummary_skim(sweden_data %>% dplyr::select(stflife, hincfel, edulvla))
table_not_cleaned| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max | ||
|---|---|---|---|---|---|---|---|---|
| stflife | 14 | 0 | 7.9 | 3.8 | 0.0 | 8.0 | 99.0 | |
| hincfel | 7 | 0 | 1.6 | 0.9 | 1.0 | 1.0 | 9.0 | |
| edulvla | 8 | 58 | 3.3 | 5.2 | 1.0 | 3.0 | 99.0 |
Clearer labels + Final numeric table:
sweden_clean <- sweden_clean %>%
# Relabeling variables for clarity
mutate(
"Satisfaction with life as a whole (stflife)" = stflife, # Renaming 'stflife' to 'Satisfication with life as a whole'
"Feeling about household's income nowadays (hincfel)" = hincfel, # Renaming 'hincfel' to 'Feeling about household's income nowadays'
"Highest level of education (edulvla)" = edulvla # Renaming 'edulvla' to 'Highest level of education'
)table1b <- datasummary_skim(sweden_clean %>% dplyr::select("Satisfaction with life as a whole (stflife)", "Feeling about household's income nowadays (hincfel)", "Highest level of education (edulvla)"), title = "Table 1. Descriptive statistics of stflife, hincfel, and edulvla in Sweden", 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.table1b
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
Satisfaction with life as a whole (stflife) | 10 | 0 | 8.0 | 1.6 | 0.0 | 8.0 | 10.0 |
Feeling about household's income nowadays (hincfel) | 4 | 0 | 1.5 | 0.7 | 1.0 | 1.0 | 4.0 |
Highest level of education (edulvla) | 4 | 0 | 3.0 | 1.3 | 1.0 | 3.0 | 5.0 |
Saving the numeric table:
set_flextable_defaults(fonts_ignore=TRUE)
flextable::save_as_docx(table1b, path = "Table3.docx",
width = 7.0, height = 7.0) The data can also be inputted into a categorical summary:
sweden_clean2 <- sweden_clean %>%
mutate(
stflife = case_when(
stflife %in% 0:4 ~ "Poor satisfaction", # Values 0 to 3 in 'stflife' are categorized as "Poor Satisfaction"
stflife %in% 6:10 ~ "High satisfaction", # Values 7 to 10 in 'stflife' are categorized as "Right"
stflife %in% c(5, 77, 88, 99) ~ NA_character_, # Values 5, 77, 88, 99 in 'stfife' are set as NA
TRUE ~ NA_character_
),
hincfel = case_when(
hincfel == 1 ~ "Comfortable income", # 1 is recoded to Comfortable income
hincfel == 2 ~ "Coping with income", # 2 is recoded to Coping with income
hincfel %in% 3:4 ~ "Difficult with income", # 3 and 4 are recoded to Difficult with income
hincfel %in% c(7, 8, 9) ~ NA_character_, # Values 7, 8, 9 in 'hincfel' are set as NA
TRUE ~ NA_character_
),
# Recording edulvla
edulvla = case_when(
edulvla == 1 ~ "Below secondary",
edulvla %in% 2:3 ~ "Secondary education",
edulvla %in% 4:5 ~ "Post-secondary education",
edulvla %in% c(55, 77, 88, 99) ~ NA_character_, # Values 55, 77, 88, 99 in 'edulvla' are set as NA
TRUE ~ NA_character_
))
sweden_clean2 <- sweden_clean2 %>% filter(!is.na(stflife), !is.na(hincfel), !is.na(edulvla)) # ensuring there are no NAstable1_cat <- datasummary_skim(sweden_clean2 %>% dplyr::select(stflife, hincfel, edulvla), type = "categorical", title = "Table 1. Categorical Variables Summary", output = "flextable")
table1_cat
|
| N | % |
|---|---|---|---|
stflife | High satisfaction | 6834 | 94.9 |
Poor satisfaction | 370 | 5.1 | |
hincfel | Comfortable income | 4156 | 57.7 |
Coping with income | 2466 | 34.2 | |
Difficult with income | 582 | 8.1 | |
edulvla | Below secondary | 1045 | 14.5 |
Post-secondary education | 1769 | 24.6 | |
Secondary education | 4390 | 60.9 |
Saving the categorical table:
set_flextable_defaults(fonts_ignore=TRUE)
flextable::save_as_docx(table1_cat, path = "Categorical table3.docx",
width = 20.0, height = 10.0) INFER PACKAGE SHOWING NULL DISTRIBUTION
The outcome/response variable (life satisfication) is a categorical variables with 2 categories (Poor or High). The main predictor (income) is a categorical variable with 3 categories (Below secondary, secondary, or post-secondary). Thus, we can use a Chi-square test of independence.
Step 1 -
ess <- read_fst("All-ESS-Data.fst")sweden_data <- ess %>%
# Filter for France
filter(cntry == "SE") %>%
# Recoding stflife and cleaning trstplt (will use later)
mutate(
stflife = case_when(
stflife %in% 0:4 ~ "Poor satisfication", # Values 0 to 3 in 'stflife' are categorized as "Left"
stflife %in% 6:10 ~ "High satisfication", # Values 7 to 10 in 'stflife' are categorized as "Right"
stflife %in% c(5, 77, 88, 99) ~ NA_character_ # Values 77, 88, 99 in 'stfife' are set as NA
),
# Recoding 'yrbrn' (year born) into a new variable 'Gen'
hincfel = case_when(
hincfel == 1 ~ "Living comfortably on present income",
hincfel == 2 ~ "Coping on present income",
hincfel %in% 3:4 ~ "Difficult with income",
hincfel %in% c(7, 8, 9) ~ NA_character_ # Values 7, 8, 9 in 'hincfel' are set as NA
),
# Recording edulvla
edulvla = case_when(
edulvla == 1 ~ "Lower secondary",
edulvla %in% 2:3 ~ "Secondary education",
edulvla %in% 4:5 ~ "Tertiary education",
edulvla %in% c(55, 77, 88, 99) ~ NA_character_, # Values 7, 8, 9 in 'hincfel' are set as NA
TRUE ~ NA_character_
))
sweden_data <- sweden_data %>% filter(!is.na(stflife), !is.na(hincfel), !is.na(edulvla))First variable -
sweden_data <- sweden_data %>% filter(!is.na(stflife))unique(sweden_data$stflife)## [1] "High satisfication" "Poor satisfication"Second variable -
sweden_data <- sweden_data %>% filter(!is.na(hincfel))unique(sweden_data$hincfel)## [1] "Coping on present income"
## [2] "Living comfortably on present income"
## [3] "Difficult with income"Third variable -
sweden_data <- sweden_data %>% filter(!is.na(edulvla))unique(sweden_data$edulvla)## [1] "Lower secondary" "Tertiary education" "Secondary education"Finding critical value for stflife and edulvla:
test_stat <- sweden_data %>%
specify(explanatory = hincfel,
response = stflife) %>%
hypothesize(null = "independence") %>%
calculate(stat = "Chisq")
print(test_stat$stat)## X-squared
## 349.0627Step 2 -
null_dist <- sweden_data %>%
specify(response = stflife, explanatory = hincfel) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")
null_dist## Response: stflife (factor)
## Explanatory: hincfel (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 5.08
## 2 2 1.33
## 3 3 0.0434
## 4 4 0.639
## 5 5 14.3
## 6 6 1.47
## 7 7 0.842
## 8 8 3.99
## 9 9 0.559
## 10 10 2.47
## # ℹ 990 more rowsStep 3 -
p_val <- null_dist %>% #
get_pvalue(obs_stat = test_stat, direction = "greater")## 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 -
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)## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
The test statistic is too extreme so the default visualization does not display well. Here is a closer look:
test_stat <- 960
conf_int <- null_dist %>%
get_confidence_interval(level = 0.95, type = "percentile")
# Visualize the null distribution with manual x-axis limits
null_dist %>%
visualize(bins = 15, method = "simulation", dens_color = "black") +
geom_vline(aes(xintercept = test_stat), color = "red", linetype = "dashed", size = 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)) +
labs(title = "Figure 1. Simulation-Based Null Distribution",
x = "Chi-Square",
y = "Count")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
The plot above indicates that test statistic falls extremely out the simulated null distribution, showing evidence to reject the null hypothesis.
REGRESSION MODELS
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey","MASS", "aod", "interactions", "kableExtra", "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)## 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## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools##
## Attaching package: 'scales'## The following object is masked from 'package:viridis':
##
## viridis_pal## 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## [[1]]
## [1] "fstcore" "questionr" "ggridges" "rmarkdown" "kableExtra"
## [6] "knitr" "viridis" "viridisLite" "fst" "infer"
## [11] "flextable" "modelsummary" "lubridate" "forcats" "stringr"
## [16] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [21] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [26] "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "fstcore" "questionr" "ggridges" "rmarkdown" "kableExtra"
## [6] "knitr" "viridis" "viridisLite" "fst" "infer"
## [11] "flextable" "modelsummary" "lubridate" "forcats" "stringr"
## [16] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [21] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [26] "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "fstcore" "questionr" "ggridges" "rmarkdown" "kableExtra"
## [6] "knitr" "viridis" "viridisLite" "fst" "infer"
## [11] "flextable" "modelsummary" "lubridate" "forcats" "stringr"
## [16] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [21] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [26] "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "fstcore" "questionr" "ggridges" "rmarkdown" "kableExtra"
## [6] "knitr" "viridis" "viridisLite" "fst" "infer"
## [11] "flextable" "modelsummary" "lubridate" "forcats" "stringr"
## [16] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [21] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [26] "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "effects" "carData" "fstcore" "questionr" "ggridges"
## [6] "rmarkdown" "kableExtra" "knitr" "viridis" "viridisLite"
## [11] "fst" "infer" "flextable" "modelsummary" "lubridate"
## [16] "forcats" "stringr" "dplyr" "purrr" "readr"
## [21] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [26] "graphics" "grDevices" "utils" "datasets" "methods"
## [31] "base"
##
## [[6]]
## [1] "survey" "survival" "Matrix" "grid" "effects"
## [6] "carData" "fstcore" "questionr" "ggridges" "rmarkdown"
## [11] "kableExtra" "knitr" "viridis" "viridisLite" "fst"
## [16] "infer" "flextable" "modelsummary" "lubridate" "forcats"
## [21] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [26] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [31] "grDevices" "utils" "datasets" "methods" "base"
##
## [[7]]
## [1] "MASS" "survey" "survival" "Matrix" "grid"
## [6] "effects" "carData" "fstcore" "questionr" "ggridges"
## [11] "rmarkdown" "kableExtra" "knitr" "viridis" "viridisLite"
## [16] "fst" "infer" "flextable" "modelsummary" "lubridate"
## [21] "forcats" "stringr" "dplyr" "purrr" "readr"
## [26] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [31] "graphics" "grDevices" "utils" "datasets" "methods"
## [36] "base"
##
## [[8]]
## [1] "aod" "MASS" "survey" "survival" "Matrix"
## [6] "grid" "effects" "carData" "fstcore" "questionr"
## [11] "ggridges" "rmarkdown" "kableExtra" "knitr" "viridis"
## [16] "viridisLite" "fst" "infer" "flextable" "modelsummary"
## [21] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [26] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [31] "stats" "graphics" "grDevices" "utils" "datasets"
## [36] "methods" "base"
##
## [[9]]
## [1] "interactions" "aod" "MASS" "survey" "survival"
## [6] "Matrix" "grid" "effects" "carData" "fstcore"
## [11] "questionr" "ggridges" "rmarkdown" "kableExtra" "knitr"
## [16] "viridis" "viridisLite" "fst" "infer" "flextable"
## [21] "modelsummary" "lubridate" "forcats" "stringr" "dplyr"
## [26] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [31] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [36] "datasets" "methods" "base"
##
## [[10]]
## [1] "interactions" "aod" "MASS" "survey" "survival"
## [6] "Matrix" "grid" "effects" "carData" "fstcore"
## [11] "questionr" "ggridges" "rmarkdown" "kableExtra" "knitr"
## [16] "viridis" "viridisLite" "fst" "infer" "flextable"
## [21] "modelsummary" "lubridate" "forcats" "stringr" "dplyr"
## [26] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [31] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [36] "datasets" "methods" "base"
##
## [[11]]
## [1] "performance" "interactions" "aod" "MASS" "survey"
## [6] "survival" "Matrix" "grid" "effects" "carData"
## [11] "fstcore" "questionr" "ggridges" "rmarkdown" "kableExtra"
## [16] "knitr" "viridis" "viridisLite" "fst" "infer"
## [21] "flextable" "modelsummary" "lubridate" "forcats" "stringr"
## [26] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [31] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [36] "utils" "datasets" "methods" "base"
##
## [[12]]
## [1] "performance" "interactions" "aod" "MASS" "survey"
## [6] "survival" "Matrix" "grid" "effects" "carData"
## [11] "fstcore" "questionr" "ggridges" "rmarkdown" "kableExtra"
## [16] "knitr" "viridis" "viridisLite" "fst" "infer"
## [21] "flextable" "modelsummary" "lubridate" "forcats" "stringr"
## [26] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [31] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [36] "utils" "datasets" "methods" "base"
##
## [[13]]
## [1] "broom" "performance" "interactions" "aod" "MASS"
## [6] "survey" "survival" "Matrix" "grid" "effects"
## [11] "carData" "fstcore" "questionr" "ggridges" "rmarkdown"
## [16] "kableExtra" "knitr" "viridis" "viridisLite" "fst"
## [21] "infer" "flextable" "modelsummary" "lubridate" "forcats"
## [26] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [31] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [36] "grDevices" "utils" "datasets" "methods" "base"
##
## [[14]]
## [1] "scales" "broom" "performance" "interactions" "aod"
## [6] "MASS" "survey" "survival" "Matrix" "grid"
## [11] "effects" "carData" "fstcore" "questionr" "ggridges"
## [16] "rmarkdown" "kableExtra" "knitr" "viridis" "viridisLite"
## [21] "fst" "infer" "flextable" "modelsummary" "lubridate"
## [26] "forcats" "stringr" "dplyr" "purrr" "readr"
## [31] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [36] "graphics" "grDevices" "utils" "datasets" "methods"
## [41] "base"
##
## [[15]]
## [1] "ggeffects" "scales" "broom" "performance" "interactions"
## [6] "aod" "MASS" "survey" "survival" "Matrix"
## [11] "grid" "effects" "carData" "fstcore" "questionr"
## [16] "ggridges" "rmarkdown" "kableExtra" "knitr" "viridis"
## [21] "viridisLite" "fst" "infer" "flextable" "modelsummary"
## [26] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [31] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [36] "stats" "graphics" "grDevices" "utils" "datasets"
## [41] "methods" "base"
##
## [[16]]
## [1] "marginaleffects" "ggeffects" "scales" "broom"
## [5] "performance" "interactions" "aod" "MASS"
## [9] "survey" "survival" "Matrix" "grid"
## [13] "effects" "carData" "fstcore" "questionr"
## [17] "ggridges" "rmarkdown" "kableExtra" "knitr"
## [21] "viridis" "viridisLite" "fst" "infer"
## [25] "flextable" "modelsummary" "lubridate" "forcats"
## [29] "stringr" "dplyr" "purrr" "readr"
## [33] "tidyr" "tibble" "ggplot2" "tidyverse"
## [37] "stats" "graphics" "grDevices" "utils"
## [41] "datasets" "methods" "base"sweden_data <- ess %>%
# Filter for Sweden
filter(cntry == "SE")
sweden_data$weight <- sweden_data$dweight * sweden_data$pweight
survey_design <- svydesign(ids = ~1, data = sweden_data, weights = ~weight)
sweden_data <- sweden_data %>%
mutate(
stflife = case_when(
stflife %in% 0:4 ~ "Poor satisfication", # Values 0 to 3 in 'stflife' are categorized as "Left"
stflife %in% 6:10 ~ "High satisfication", # Values 7 to 10 in 'stflife' are categorized as "Right"
stflife %in% c(5, 77, 88, 99) ~ NA_character_, # Values 5, 77, 88, 99 in 'stfife' are set as NA
TRUE ~ as.character(stflife)
),
stflife = factor(stflife, levels = c("Poor satisfication", "High satisfication"),
labels = c("Poor satisfaction", "High satisfaction")
),
hincfel = case_when(
hincfel == 1 ~ "Comfortable",
hincfel == 2 ~ "Coping",
hincfel %in% 3:4 ~ "Difficult",
hincfel %in% c(7, 8, 9) ~ NA_character_, # Values 7, 8, 9 in 'hincfel' are set as NA
TRUE ~ NA_character_
),
hincfel = factor(hincfel, levels = c("Difficult", "Coping", "Comfortable"),
labels = c("Difficult", "Coping", "Comfortable")
),
# Recording edulvla
edulvla = case_when(
edulvla == 1 ~ "Below secondary",
edulvla %in% 2:3 ~ "Secondary education",
edulvla %in% 4:5 ~ "Tertiary education",
edulvla %in% c(55, 77, 88, 99) ~ NA_character_, # Values 55, 77, 88, 99 in 'hincfel' are set as NA
TRUE ~ NA_character_
),
edulvla = factor(edulvla, levels = c("Below secondary", "Secondary education", "Tertiary education"),
labels = c("Below secondary", "Secondary education", "Tertiary education")
))
sweden_data <- sweden_data %>% filter(!is.na(stflife), !is.na(hincfel), !is.na(edulvla))table(sweden_data$stflife) # checking to see if they are ordinal##
## Poor satisfaction High satisfaction
## 370 6834table(sweden_data$hincfel)##
## Difficult Coping Comfortable
## 582 2466 4156table(sweden_data$edulvla)##
## Below secondary Secondary education Tertiary education
## 1045 4390 1769sweden_data <- sweden_data %>%
mutate(
stflife = case_when(
stflife == "Poor satisfaction" ~ 1,
stflife == "High satisfaction" ~ 0,
TRUE ~ NA_integer_
)
)table(sweden_data$stflife) # checking if it is 0 and 1##
## 0 1
## 6834 370model1 <- glm(stflife ~ hincfel, family = binomial, data = sweden_data) # single predictor
model2 <- glm(stflife ~ hincfel + edulvla, data = sweden_data) # two predictors
model3 <- glm(stflife ~ hincfel + edulvla + hincfel*edulvla, family = binomial, data = sweden_data) # both predictors + interactionmodelsummary(
list(model1, model2, model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
exponentiate = TRUE)| (1) | (2) | (3) | |
|---|---|---|---|
| (Intercept) | 0.3 (0.0)*** | 1.2 (0.0)*** | 0.2 (0.1)*** |
| hincfelCoping | 0.3 (0.0)*** | 0.9 (0.0)*** | 0.3 (0.1)** |
| hincfelComfortable | 0.1 (0.0)*** | 0.8 (0.0)*** | 0.1 (0.0)*** |
| edulvlaSecondary education | 1.0 (0.0) | 1.3 (0.4) | |
| edulvlaTertiary education | 1.0 (0.0) | 0.9 (0.3) | |
| hincfelCoping × edulvlaSecondary education | 0.7 (0.3) | ||
| hincfelComfortable × edulvlaSecondary education | 0.9 (0.4) | ||
| hincfelCoping × edulvlaTertiary education | 1.0 (0.5) | ||
| hincfelComfortable × edulvlaTertiary education | 1.0 (0.5) | ||
| Num.Obs. | 7204 | 7204 | 7204 |
| R2 | 0.049 | ||
| AIC | 2669.3 | −1672.3 | 2676.3 |
| BIC | 2689.9 | −1631.1 | 2738.2 |
| Log.Lik. | −1331.641 | 842.175 | −1329.152 |
| RMSE | 0.22 | 0.22 | 0.22 |
Improving display:
models <- list()
models[['SLR']]<- glm(stflife ~ hincfel, family = binomial, data = sweden_data) # single predictor
models[['MLR']]<- glm(stflife ~ hincfel + edulvla, family = binomial, data = sweden_data) # two predictors
models[['Interaction']]<- glm(stflife ~ hincfel + edulvla + hincfel*edulvla, family = binomial, data = sweden_data) # both predictors + interaction
modelsummary(models,
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
exponentiate = TRUE,
coef_rename = c("hincfelCoping" = "Coping with income", "hincfelComfortable" = "Comfortable", "edulvlaSecondary education" = "Secondary", "edulvlaTertiary education" = "Tertiary"),
title = 'Table 2. Regression models predicting life satisfaction')| SLR | MLR | Interaction | |
|---|---|---|---|
| (Intercept) | 0.3 (0.0)*** | 0.2 (0.0)*** | 0.2 (0.1)*** |
| Coping with income | 0.3 (0.0)*** | 0.3 (0.0)*** | 0.3 (0.1)** |
| Comfortable | 0.1 (0.0)*** | 0.1 (0.0)*** | 0.1 (0.0)*** |
| Secondary | 1.1 (0.2) | 1.3 (0.4) | |
| Tertiary | 0.9 (0.2) | 0.9 (0.3) | |
| Coping with income:Secondary | 0.7 (0.3) | ||
| Comfortable:Secondary | 0.9 (0.4) | ||
| Coping with income:Tertiary | 1.0 (0.5) | ||
| Comfortable:Tertiary | 1.0 (0.5) | ||
| Num.Obs. | 7204 | 7204 | 7204 |
| AIC | 2669.3 | 2670.0 | 2676.3 |
| BIC | 2689.9 | 2704.4 | 2738.2 |
| Log.Lik. | −1331.641 | −1329.994 | −1329.152 |
| RMSE | 0.22 | 0.22 | 0.22 |
Interaction effect visual
interaction_plot <- effect("hincfel*edulvla", model3, na.rm=TRUE)
plot(interaction_plot,
main="Figure 2. Interaction effect",
xlab="edulvla",
ylab="stflife") # Scale goes from high to poor satisfaction because high was numbered 0 and poor was numbered 1