Loading Data

packages <- c("tidyverse", "infer", "fst", "modelsummary", "survey", "interactions", "scales", "broom", "dplyr", "glmnet") 

new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_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: 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: 'scales'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
## 
## 
## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools
## 
## Loaded glmnet 4.1-8
## [[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] "survey"       "survival"     "Matrix"       "grid"         "modelsummary"
##  [6] "fst"          "infer"        "lubridate"    "forcats"      "stringr"     
## [11] "dplyr"        "purrr"        "readr"        "tidyr"        "tibble"      
## [16] "ggplot2"      "tidyverse"    "stats"        "graphics"     "grDevices"   
## [21] "utils"        "datasets"     "methods"      "base"        
## 
## [[6]]
##  [1] "interactions" "survey"       "survival"     "Matrix"       "grid"        
##  [6] "modelsummary" "fst"          "infer"        "lubridate"    "forcats"     
## [11] "stringr"      "dplyr"        "purrr"        "readr"        "tidyr"       
## [16] "tibble"       "ggplot2"      "tidyverse"    "stats"        "graphics"    
## [21] "grDevices"    "utils"        "datasets"     "methods"      "base"        
## 
## [[7]]
##  [1] "scales"       "interactions" "survey"       "survival"     "Matrix"      
##  [6] "grid"         "modelsummary" "fst"          "infer"        "lubridate"   
## [11] "forcats"      "stringr"      "dplyr"        "purrr"        "readr"       
## [16] "tidyr"        "tibble"       "ggplot2"      "tidyverse"    "stats"       
## [21] "graphics"     "grDevices"    "utils"        "datasets"     "methods"     
## [26] "base"        
## 
## [[8]]
##  [1] "broom"        "scales"       "interactions" "survey"       "survival"    
##  [6] "Matrix"       "grid"         "modelsummary" "fst"          "infer"       
## [11] "lubridate"    "forcats"      "stringr"      "dplyr"        "purrr"       
## [16] "readr"        "tidyr"        "tibble"       "ggplot2"      "tidyverse"   
## [21] "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
## [26] "methods"      "base"        
## 
## [[9]]
##  [1] "broom"        "scales"       "interactions" "survey"       "survival"    
##  [6] "Matrix"       "grid"         "modelsummary" "fst"          "infer"       
## [11] "lubridate"    "forcats"      "stringr"      "dplyr"        "purrr"       
## [16] "readr"        "tidyr"        "tibble"       "ggplot2"      "tidyverse"   
## [21] "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
## [26] "methods"      "base"        
## 
## [[10]]
##  [1] "glmnet"       "broom"        "scales"       "interactions" "survey"      
##  [6] "survival"     "Matrix"       "grid"         "modelsummary" "fst"         
## [11] "infer"        "lubridate"    "forcats"      "stringr"      "dplyr"       
## [16] "purrr"        "readr"        "tidyr"        "tibble"       "ggplot2"     
## [21] "tidyverse"    "stats"        "graphics"     "grDevices"    "utils"       
## [26] "datasets"     "methods"      "base"
hungary_data <- read_fst("hungary_data.fst")
hungary_data <- hungary_data %>%
  filter(cntry == "HU")
hungary_data$weight <- hungary_data$dweight * hungary_data$pweight
survey_design <- svydesign(ids = ~1, data = hungary_data, weights = ~weight)
hungary_data <- hungary_data %>%
  mutate(healthcare = case_when(
    medtrnp == 0 ~ 0,    # Not marked to 0
    medtrnp == 1 ~ 1,    # Marked to 1
    TRUE ~ NA_real_    # Everything else to NA
  ))
hungary_data <- hungary_data %>% filter(!is.na(healthcare))
hungary_data <- hungary_data %>%
  # Recoding income and education
  mutate(
    income = case_when(
      hinctnta %in% 1:2 ~ "Low-income",                    
      hinctnta %in% 3:5 ~ "Medium-income",                     
      hinctnta %in% 6:10 ~ "High-income",                  
      hinctnta %in% c(77, 88, 99) ~ NA_character_      
    ),
    educ.ba = case_when(
      essround == 5 & edulvla == 5 ~ "BA or more",
      essround >= 5 & edulvlb > 600 ~ "BA or more",
      TRUE ~ "No BA"
    ),
    edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla),
    edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
    educ.ba = factor(educ.ba, levels = c("No BA", "BA or more"))
  )
hungarydf <- hungary_data %>%
  dplyr::select(healthcare, income, educ.ba)

hungarydf <- na.omit(hungarydf)
head(hungarydf)
##   healthcare        income    educ.ba
## 1          0 Medium-income BA or more
## 2          0   High-income      No BA
## 3          0   High-income      No BA
## 4          0 Medium-income      No BA
## 5          0 Medium-income      No BA
## 6          0 Medium-income      No BA

Descriptive table

hungary_data <- hungary_data %>%
  mutate(
    Healthcare = healthcare,
    Income = case_when(
      hinctnta %in% 1:2 ~ 1, # Low-income                   
      hinctnta %in% 3:5 ~ 2, # Medium-income                     
      hinctnta %in% 6:10 ~ 3, # High-income                  
      hinctnta %in% c(77, 88, 99) ~ NA_integer_,     
      TRUE ~ NA_real_
    ),
    Educ.ba = case_when(
      essround == 5 & edulvla == 5 ~ 1, # BA or more
      essround >= 5 & edulvlb > 600 ~ 1, # BA or more"
      TRUE ~ 0 # No BA
    ),
  )
df <- hungary_data
df <- df %>% filter(!is.na(Income))
df <- df %>% filter(!is.na(educ.ba))

table1 <- datasummary_skim(df  %>% dplyr::select(Healthcare, Income, Educ.ba), title = "Table 1. Descriptive Statistics for Main Variables (Healthcare, Income, and Education)", 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.
table1
Table 1. Descriptive Statistics for Main Variables (Healthcare, Income, and Education)

Unique (#)

Missing (%)

Mean

SD

Min

Median

Max

Healthcare

2

0

0.0

0.1

0.0

0.0

1.0

Income

3

0

2.5

0.7

1.0

3.0

3.0

Educ.ba

2

0

0.2

0.4

0.0

0.0

1.0

Null Distribution Model

ess <- hungary_data
ess <- ess %>%
  mutate(healthcare = case_when(
    medtrnp == 1 ~ 1,    # Marked
    medtrnp == 0 ~ 0,    # Not Marked
    TRUE ~ NA_real_    
  ))
ess <- ess %>% filter(!is.na(healthcare))
ess <- ess %>%
  mutate(
    income = case_when(
    hinctnta %in% 1:2 ~ "Low-income",
    hinctnta %in% 3:5 ~ "Middle Income",
    hinctnta %in% 6:10 ~ "High-income",
    hinctnta %in% c(77, 88, 99) ~ NA_character_
    ),
    
    educ.ba = case_when(
      essround == 7 & edulvla == 7 ~ "BA or more",
      essround == 7 & edulvlb > 600 ~ "BA or more",
      TRUE ~ "No BA"
    ),
    
    edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla),
    edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
    educ.ba = factor(educ.ba, levels = c("No BA", "BA or more")),
  )
test_stat <- ess %>%
  specify(explanatory = income, 
          response = healthcare) %>%     
  calculate(stat = "F") 
## Warning: Removed 491 rows containing missing values.
print(test_stat$stat)
## [1] 1.584746
test_stat <- ess %>%
  specify(explanatory = edulvlb, 
          response = healthcare) %>%     
  calculate(stat = "slope") 
## Warning: Removed 8 rows containing missing values.
print(test_stat$stat)
##      edulvlb 
## 7.129716e-06
null_dist <- ess %>%
  specify(explanatory = hinctnta, 
          response = medtrnp) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "slope")

null_dist
## Response: medtrnp (numeric)
## Explanatory: hinctnta (numeric)
## # A tibble: 1,000 × 2
##    replicate       stat
##        <int>      <dbl>
##  1         1 -0.0000971
##  2         2 -0.0000621
##  3         3 -0.0000969
##  4         4 -0.000119 
##  5         5 -0.0000897
##  6         6 -0.000144 
##  7         7 -0.000151 
##  8         8 -0.000115 
##  9         9 -0.000130 
## 10        10 -0.000104 
## # ℹ 990 more rows
p_val <- null_dist %>% # 
  get_pvalue(obs_stat = test_stat, direction = "two-sided")
p_val
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1    0.08
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 = "two-sided") +
  shade_confidence_interval(endpoints = conf_int)

Regression Models

ess$weight <- ess$dweight * ess$pweight
ess <- ess %>%
  dplyr::select(healthcare, income, educ.ba, weight)
table(ess$healthcare)
## 
##    0    1 
## 1683   15
model1 <- lm(healthcare ~ income, data = ess, weights = weight)  # single predictor 

model2 <- lm(healthcare ~ income + educ.ba, data = ess, weights = weight) # two predictors 

model3 <- lm(healthcare ~ income + educ.ba + income*educ.ba, data = ess, weights = weight) # both predictors + interaction
modelsummary(
  list(model1, model2, model3),
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL)
 (1)   (2)   (3)
(Intercept) 0.0 (0.0)* 0.0 (0.0)+ 0.0 (0.0)*
incomeLow-income 0.0 (0.0)+ 0.0 (0.0)+ 0.0 (0.0)
incomeMiddle Income 0.0 (0.0) 0.0 (0.0) 0.0 (0.0)
educ.baBA or more 0.0 (0.0) 0.0 (0.0)
incomeLow-income × educ.baBA or more 0.2 (0.0)***
incomeMiddle Income × educ.baBA or more 0.0 (0.0)+
Num.Obs. 1207 1207 1207
R2 0.003 0.003 0.017
R2 Adj. 0.001 0.000 0.013
AIC −2051.8 −2049.8 −2063.7
BIC −2031.4 −2024.4 −2028.0
Log.Lik. 1029.888 1029.924 1038.856
RMSE 0.10 0.10 0.10

Visual Of Choice

hist(ess$healthcare)