getwd()
## [1] "/Users/gabrieladrzymala/Desktop/SOC222 R STUFF "

#Loading packages

# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "performance", "flextable", "broom", "scales", "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)
## ── 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: 'flextable'
## 
## 
## 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] "performance"  "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] "flextable"    "performance"  "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] "broom"        "flextable"    "performance"  "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] "scales"       "broom"        "flextable"    "performance"  "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] "marginaleffects" "scales"          "broom"           "flextable"      
##  [5] "performance"     "survey"          "survival"        "Matrix"         
##  [9] "grid"            "effects"         "carData"         "modelsummary"   
## [13] "fst"             "infer"           "lubridate"       "forcats"        
## [17] "stringr"         "dplyr"           "purrr"           "readr"          
## [21] "tidyr"           "tibble"          "ggplot2"         "tidyverse"      
## [25] "stats"           "graphics"        "grDevices"       "utils"          
## [29] "datasets"        "methods"         "base"

#Loading data

estonia_data <- read.fst("estonia_data.fst")
library(fst)
library(modelsummary)
estonia_data <- estonia_data %>%
  filter(cntry == "EE")
estonia_data$weight <- estonia_data$dweight * estonia_data$pweight
survey_design <- svydesign(ids = ~1, data = estonia_data, weights = ~weight)
estonia_education <- estonia_data %>%
  filter(cntry == "EE") %>%
  dplyr::select(edulvla)
estonia_income <- estonia_data %>%
  filter(cntry == "EE") %>%
  dplyr::select(hincfel)
estonia_gender <- estonia_data %>%
  filter(cntry == "EE") %>%
  dplyr::select(gndr)

#RECODING VARIABLES

estoniadf <- estonia_data %>%
  filter(cntry == "EE")
estonia_data$weight <- estonia_data$dweight * estonia_data$pweight
survey_design <- svydesign(ids = ~1, data = estonia_data, weights = ~weight)
estonia_data <- estonia_data %>%
  mutate(income = case_when(
    hincfel == 1 ~ 1,    
    hincfel == 2 ~ 2,
    hincfel == 3 ~ 3,
    hincfel == 4 ~ 4,
    TRUE ~ NA_real_    # the rest of the data collected to NA 
  ))
estonia_data <- estonia_data %>% filter(!is.na(income))
estonia_data <- estonia_data %>%
  mutate(
    gndr = case_when (
      gndr == 1 ~ "Male",    
    gndr == 2 ~ "Female", 
      gndr %in% c(9) ~ NA_character_,
      TRUE ~ as.character(gndr)
  )) %>%
  filter (!is.na(gndr))
estoniadf <- estonia_data %>%
  mutate(
    educ.ba = case_when(
      essround < 5 & edulvla == 5 ~ "BA or more", #setting the lowest number at 5 
      essround >= 5 & edulvlb > 500 ~ "BA or more", #setting the highest number to 500
      TRUE ~ "No BA"
    ),
    edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla), #setting intergers as 77,88,99
    edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb), #setting the rest of the integers to 5555,7777,8888
    educ.ba = factor(educ.ba, levels = c("No BA", "BA or more"))
  )
table(estoniadf$educ.ba)
## 
##      No BA BA or more 
##      11644       5145

#Categorical Descriptive Table

df <- estoniadf %>%
  mutate(income = case_when(
    hincfel == 1 ~ "Living comfortably on present income",    
    hincfel == 2 ~ "Coping on present income", 
    hincfel == 3 ~ "Difficult on present income",
    hincfel == 4 ~ "Very difficult on present income",
    TRUE ~ NA_character_     
  ))

df <- estoniadf %>%
  # renanming the variables 
  mutate(
    hincfel = income,      
    education = educ.ba,      
    gender = gndr
  )
table1 <- datasummary_skim(df  %>% dplyr::select(income, educ.ba, gndr), type = "categorical",
                           title = "Table 1. Main Variables", output = "flextable")

table1
Table 1. Main Variables

N

%

educ.ba

No BA

11644

69.4

BA or more

5145

30.6

gndr

Female

9629

57.4

Male

7160

42.6

set_flextable_defaults(fonts_ignore=TRUE)
flextable::save_as_docx(table1, path = "Table1.docx", 
                       width = 7.0, height = 7.0)

#Step 1: Calculate the test statistic of your sample

test_stat <- estonia_data %>%
  specify(explanatory = gndr, 
          response = income) %>% 
  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 "Female" - "Male", or divided in the order "Female" / "Male" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Female", "Male")` to the calculate() function.
print(test_stat$stat)
##       t 
## 8.66783

#Step 2: Simulate the null distribution

null_distribution <- estonia_data %>%
  specify(explanatory = gndr,
          response = income) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 500, 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 "Female" - "Male", or divided in the order "Female" / "Male" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("Female", "Male")` to the calculate() function.

#Step 3: Calculate the p-value of your sample

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

#Step 4: Visualize without confidence intervals

null_distribution %>%
  visualize() +
  shade_p_value(obs_stat = test_stat, direction = "two-sided")

null_distribution
## Response: income (numeric)
## Explanatory: gndr (factor)
## Null Hypothesis: independence
## # A tibble: 500 × 2
##    replicate    stat
##        <int>   <dbl>
##  1         1  1.01  
##  2         2 -2.11  
##  3         3  0.0369
##  4         4  0.950 
##  5         5 -0.592 
##  6         6 -0.712 
##  7         7 -0.248 
##  8         8 -0.652 
##  9         9  0.402 
## 10        10 -0.389 
## # ℹ 490 more rows

#Step 4: Visualize with confidence intervals

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)

#Linear regression models

estoniadf <- estoniadf %>%
  mutate(income = case_when(
    hincfel == 1 ~ "Living comfortably on present income",    
    hincfel == 2 ~ "Coping on present income", 
    hincfel == 3 ~ "Difficult on present income",
    hincfel == 4 ~ "Very difficult on present income",
    TRUE ~ NA_character_     # Keep NA values without changing it
  ))
model1 <- lm(hincfel ~ gndr, data = df, weights = weight)
model2 <- lm(hincfel ~ gndr + edulvla, data = df, weights = weight)
modelsummary(
  list(model1, model2), 
  fmt = 1,
  estimate  = c( "{estimate} ({std.error}){stars}", # what is inputed here is what is supposed to be shown on the chart. 
                  "{estimate} ({std.error}){stars}"),
  statistic = NULL,

   title = "Figure 1: Regression model for Income, Gender, and Education in Estonia")
Figure 1: Regression model for Income, Gender, and Education in Estonia
 (1)   (2)
(Intercept) 2.3 (0.0)*** 2.7 (0.0)***
gndrMale −0.1 (0.0)*** −0.1 (0.0)***
edulvla −0.1 (0.0)***
Num.Obs. 16789 5106
R2 0.005 0.031
R2 Adj. 0.004 0.031
AIC 38809.5 11492.0
BIC 38832.7 11518.2
Log.Lik. −19401.763 −5742.016
RMSE 0.77 0.75

#Did not get third model in the table so did some troubleshooting

model3 <- lm(hincfel ~ gndr + edulvla + gndr*edulvla, data = df, weights = weight)
modelsummary(
  list(model1, model2, model3),
  fmt = 3,
  estimate = c("{estimate} ({std.error}){stars}",
               "{estimate} ({std.error}){stars}",
               "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  title = "Figure 2: Regression model for Income, Gender, and Education in Estonia")
Figure 2: Regression model for Income, Gender, and Education in Estonia
 (1)   (2)   (3)
(Intercept) 2.287 (0.008)*** 2.732 (0.031)*** 2.765 (0.039)***
gndrMale −0.104 (0.012)*** −0.122 (0.021)*** −0.206 (0.060)***
edulvla −0.099 (0.008)*** −0.109 (0.011)***
gndrMale × edulvla 0.026 (0.017)
Num.Obs. 16789 5106 5106
R2 0.005 0.031 0.032
R2 Adj. 0.004 0.031 0.031
AIC 38809.5 11492.0 11491.8
BIC 38832.7 11518.2 11524.5
Log.Lik. −19401.763 −5742.016 −5740.903
RMSE 0.77 0.75 0.75

#Visual of choice

estonia_plot <- estonia_data %>%
  mutate(
    hincfel = case_when(
      hincfel == 1 ~ "Living comfortably on present income",    
    hincfel == 2 ~ "Coping on present income", 
    hincfel == 3 ~ "Difficult on present income",
    hincfel == 4 ~ "Very difficult on present income",
    TRUE ~ NA_character_)
    ) %>%
  filter(!is.na(hincfel))
estonia_plot <- estonia_plot %>%
  mutate (gndr = case_when (
    gndr == 1 ~ "Male",
    gndr == 2 ~ "Female",
    gndr %in% c(9) ~ NA_character_,
  )) %>%
  filter(!is.na(gndr))
estonia_plot <- estonia_plot %>%
  mutate(
        educ.ba = case_when(
      essround < 5 & edulvla == 5 ~ "BA or more", #setting the lowest number at 5 
      essround >= 5 & edulvlb > 500 ~ "BA or more", #setting the highest number to 500
      TRUE ~ "No BA"
    ),
    edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla), #setting intergers as 77,88,99
    edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb), #setting the rest of the integers to 5555,7777,8888
    educ.ba = factor(educ.ba, levels = c("No BA", "BA or more"))
    ) %>%
  filter(!is.na(edulvla))
library(ggplot2)
library(jtools)
## 
## Attaching package: 'jtools'
## The following object is masked from 'package:flextable':
## 
##     theme_apa
# Generate the table and calculate conditional proportions
table(estonia_data$hincfel, estonia_data$gndr) %>%
 
  as.data.frame() %>%
  filter(Var1 != "Total",
         Var2 != "All") %>%
  ggplot(aes(x=Var1, y=Freq, fill=Var2)) +
  geom_col(position = "dodge") +
  # Manually specify colors for each gender category
  scale_fill_manual(values = c("Male" = "blue", "Female" = "pink")) + 
  labs(title="Income in Estonia by Gender",
       y = "Number of respondents",
       x = "Living comfortably on present income to very difficult on present income ",
       fill = "Gender")