# 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"

Loading Data

sweden_data <- read_fst("sweden_data.fst")
ess <- read_fst("All-ESS-Data.fst")
sweden_data <- sweden_data %>%
  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(health = case_when(
    hlthhmp == 1 ~ 1,    # the number 1 is for better mental health 
    hlthhmp == 2 ~ 0,    # 0 setting to be a lower mental health 
    TRUE ~ NA_real_    # the rest of the data collected to NA 
  ))
sweden_data <- sweden_data %>% filter(!is.na(health))
sweden_data <- sweden_data %>%
  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(education = case_when(
    edlvdse == 1 ~ 1,    # the number 1 is for higer education 
    edlvdse == 2 ~ 0,    # 0 setting to be a lower education 
    TRUE ~ NA_real_    # the rest of the data collected to NA 
  ))
sweden_data <- sweden_data %>% filter(!is.na(education))

Recoding and Cleaning Data:

sweden_clean <- sweden_data %>% 
  mutate(
    edlvdse = ifelse(edlvdse == 2, 0, ifelse(edlvdse %in% c(7, 8, 9), NA, edlvdse)),
    hlthhmp = ifelse(hlthhmp %in% c(77, 88, 99), NA, hlthhmp),
  )
sweden_clean <- sweden_clean %>% filter(!is.na(edlvdse), !is.na(hlthhmp))
table(sweden_clean$edlvdse)
## 
##   0   1 
## 361  68
table(sweden_clean$hlthhmp)
## 
##   1   2 
## 111 318

Descriptive Statistics:

sweden_clean <- sweden_clean %>%
  mutate(
    health = health,
    hlthhmp = case_when(
      lrscale > 10 ~ NA_real_,     # Set values above 10 to NA
      TRUE ~ lrscale               # Keep other values as is
    ),
    education = case_when(
      edlvdse > 100 ~ NA_real_,       # Set values above 100 to NA
      edlvdse < 15 ~ NA_real_,        # Set values below 15 to NA
      TRUE ~ edlvdse                 # Keep other values as is
    )
  )
sweden_clean <- sweden_clean
sweden_clean <- sweden_clean %>% filter(!is.na(hlthhmp))
sweden_clean<- sweden_clean %>% filter(!is.na(edlvdse))

table1b <- datasummary_skim(sweden_clean  %>% dplyr::select(hlthhmp, edlvdse), 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
Table 1. Descriptive statistics for main variables

Unique (#)

Missing (%)

Mean

SD

Min

Median

Max

hlthhmp

11

0

4.7

2.8

0.0

5.0

10.0

edlvdse

2

0

0.1

0.4

0.0

0.0

1.0

Findings section Step 1:

test_stat <- sweden_clean %>%
  specify(explanatory = edlvdse, # change variable name for explanatory variable
          response = hlthhmp) %>% # change variable name for outcome of interest
  hypothesize(null = "independence") %>%
  calculate(stat = "correlation") 
## Message: The independence null hypothesis does not inform calculation of the
## observed statistic (a correlation) and will be ignored.
print(test_stat$stat)
## [1] -0.04621849

Step 2:

null_dist <- sweden_clean %>%
  specify(response = hlthhmp, explanatory = edlvdse) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "correlation")

null_dist
## Response: hlthhmp (numeric)
## Explanatory: edlvdse (numeric)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
##    replicate     stat
##        <int>    <dbl>
##  1         1 -0.0712 
##  2         2  0.0112 
##  3         3  0.00621
##  4         4  0.0387 
##  5         5 -0.0287 
##  6         6 -0.0312 
##  7         7 -0.0213 
##  8         8  0.0986 
##  9         9  0.0487 
## 10        10  0.0137 
## # ℹ 990 more rows

Step 3:

p_val <- null_dist %>% # 
  get_pvalue(obs_stat = test_stat, direction = "greater")
p_val
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1   0.829

Step 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)

Regression Model:

model1 <- lm(edlvdse ~ edlvdse, data = sweden_clean, weights = weight) 
## Warning in model.matrix.default(mt, mf, contrasts): the response appeared on
## the right-hand side and was dropped
## Warning in model.matrix.default(mt, mf, contrasts): problem with term 1 in
## model.matrix: no columns are assigned
model2 <- lm(edlvdse ~ edlvdse + hlthhmp, data = sweden_clean, weights = weight)
## Warning in model.matrix.default(mt, mf, contrasts): the response appeared on
## the right-hand side and was dropped

## Warning in model.matrix.default(mt, mf, contrasts): problem with term 1 in
## model.matrix: no columns are assigned
model3 <- lm(edlvdse ~ edlvdse + hlthhmp + edlvdse*hlthhmp, data = sweden_clean, weights = weight) 
## Warning in model.matrix.default(mt, mf, contrasts): the response appeared on
## the right-hand side and was dropped

## Warning in model.matrix.default(mt, mf, contrasts): problem with term 1 in
## model.matrix: no columns are assigned
models <- list(
   "Model 1" = lm(edlvdse ~ edlvdse, data = sweden_clean, weights = weight),
   "Model 2" = lm(edlvdse ~ edlvdse + hlthhmp, data = sweden_clean, weights = weight),
   "Model 3" = lm(edlvdse ~ edlvdse + hlthhmp + edlvdse*hlthhmp, data = sweden_clean, weights = weight))
## Warning in model.matrix.default(mt, mf, contrasts): the response appeared on
## the right-hand side and was dropped
## Warning in model.matrix.default(mt, mf, contrasts): problem with term 1 in
## model.matrix: no columns are assigned
## Warning in model.matrix.default(mt, mf, contrasts): the response appeared on
## the right-hand side and was dropped
## Warning in model.matrix.default(mt, mf, contrasts): problem with term 1 in
## model.matrix: no columns are assigned
## Warning in model.matrix.default(mt, mf, contrasts): the response appeared on
## the right-hand side and was dropped
## Warning in model.matrix.default(mt, mf, contrasts): problem with term 1 in
## model.matrix: no columns are assigned
modelsummary(models, fmt = 1,
  estimate  = c("{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL)
Model 1  Model 2  Model 3
(Intercept) 0.2 (0.0)*** 0.2 (0.0)*** 0.2 (0.0)***
hlthhmp 0.0 (0.0) 0.0 (0.0)***
edlvdse × hlthhmp 0.2 (0.0)***
Num.Obs. 410 410 410
R2 0.000 0.002 0.707
R2 Adj. 0.000 0.000 0.706
AIC 330.8 331.9 −168.8
BIC 338.9 343.9 −152.7
Log.Lik. −163.421 −162.928 88.382
RMSE 0.35 0.35 0.19