# 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.1
## ✔ 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"
ess <- read_fst("All-ESS-Data.fst")

#Outcome Variable: Happy

uk_data <- ess %>%
  filter(cntry == "GB")

uk_raw <- ess %>%
  filter(cntry == "GB")

uk_happy <- uk_data$happy
uk_happy_raw <- uk_happy


uk_happy[uk_happy %in% 77:99] <- NA
table(uk_happy)
## uk_happy
##    0    1    2    3    4    5    6    7    8    9   10 
##  113  112  240  502  670 1568 1541 3804 6099 3740 2557
mean_happy <- mean(uk_happy, na.rm = TRUE)
cat("Mean of Happy 'y' is", mean_happy, "\n")
## Mean of Happy 'y' is 7.472262
median_happy <- median(uk_happy, na.rm = TRUE)
cat("Median of Happy 'y' is", median_happy, "\n")
## Median of Happy 'y' is 8
sd_happy <- sd(uk_happy, na.rm = TRUE)
cat("Standard Deviation of Happy 'y':", sd_happy, "\n")
## Standard Deviation of Happy 'y': 1.915428
summary(uk_happy)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   7.000   8.000   7.472   9.000  10.000      33
ggplot(uk_data %>%
         mutate(happy = ifelse(happy %in% c(77, 88, 99), NA, happy)), aes(x = happy)) +
  geom_histogram(binwidth = 2) +
  labs(title = "Happiness Levels in the UK",
       x = "Happiness Levels (0-10)",
       y = "Count" )
## Warning: Removed 33 rows containing non-finite values (`stat_bin()`).

#Key Independent Variable: Imprich (Attitudes of Richness)

uk_imprich <- uk_data$imprich
table(uk_imprich)
## uk_imprich
##    1    2    3    4    5    6    7    8    9 
##  677 1631 2472 4040 8881 2792   47  136   64
uk_imprich[uk_imprich %in% 7:9] <- NA
table(uk_imprich)
## uk_imprich
##    1    2    3    4    5    6 
##  677 1631 2472 4040 8881 2792
mean_y_imprich <- mean(uk_imprich, na.rm = TRUE)
median_y_imprich <- mean(uk_imprich, na.rm = TRUE)
sd_imprich <- sd(uk_imprich, na.rm = TRUE)

cat("Mean of Imprich 'y':", mean_y_imprich, "\n")
## Mean of Imprich 'y': 4.326941
cat("Median of Imprich 'y':", median_y_imprich, "\n")
## Median of Imprich 'y': 4.326941
cat("Standard Deviation of Imprich 'y':")
## Standard Deviation of Imprich 'y':
#not effective#

ggplot(uk_data %>%
         mutate(imprich = ifelse(imprich %in% c(7, 8, 9), NA, imprich)), aes(x = imprich)) +
  geom_histogram(binwidth = 2) +
  labs(title = "Richness Attitudes in the UK",
       x = "Richness Levels (0-10)",
       y = "Count")
## Warning: Removed 486 rows containing non-finite values (`stat_bin()`).

#Key Independent Variable: Educational Attainment

table(uk_data$edulvla)
## 
##    1    2    3    5   55   77   88   99 
## 2506 1988 1046 3000  114    2   21   18
table(uk_data$edulvlb)
## 
##    0  113  212  213  222  229  313  321  323  413  423  510  520  620  720  800 
##   16 2595  583  772  142   77  527 1221 1049   42  373  102 1390 1556 1201  264 
## 5555 7777 8888 
##  279   28   67
uk_data <- uk_data %>%
  mutate(
    edulvla = na_if(edulvla, 55) %>%
      na_if(., 77) %>%
      na_if(., 88),
    
    edulvlb = na_if(edulvlb, 5555) %>%
      na_if(., 7777) %>%
      na_if(., 8888)
  ) %>%
  
  mutate(
    edulvla = case_when(
      essround < 5 & edulvla == 55 ~ NA_real_,
      TRUE ~ edulvla
    ),
    
    edulvlb = case_when(
      essround >= 5 & edulvlb %in% c(5555, 7777, 8888) ~ NA_real_,
      TRUE ~ edulvlb
    ),
    
    educ_attain = case_when(
      essround < 5 & edulvla == 5  ~ "BA",
      essround >= 5 & edulvlb > 600 ~ "BA",
      TRUE ~ "No BA"
    )
  )

table(uk_data$educ_attain)
## 
##    BA No BA 
##  6021 14958

#Table Summary

uk_clean <- uk_data %>%
  mutate(
    happy = case_when(
      happy %in% 0:5 ~ "Unhappy",
      happy %in% 6:10 ~ "Happiness",
      happy %in% c(77, 88, 99) ~ NA_character_
    ),
    
    imprich = case_when(
      imprich %in% 1:3 ~ "Important to Me",
      imprich %in% 4:6 ~ "Not that Important to Me"
    ),
  )%>%

mutate(
    edulvla = na_if(edulvla, 55) %>%
      na_if(., 77) %>%
      na_if(., 88),
    
    edulvlb = na_if(edulvlb, 5555) %>%
      na_if(., 7777) %>%
      na_if(., 8888)
  ) %>%
  
  mutate(
    edulvla = case_when(
      essround < 5 & edulvla == 55 ~ NA_real_,
      TRUE ~ edulvla
    ),
    
    edulvlb = case_when(
      essround >= 5 & edulvlb %in% c(5555, 7777, 8888) ~ NA_real_,
      TRUE ~ edulvlb
    ),
    
    EducationLevel = case_when(
      essround < 5 & edulvla == 5  ~ "BA",
      essround >= 5 & edulvlb > 600 ~ "BA",
      TRUE ~ "No BA"
    ),
    
    happy = ifelse(happy %in% c(77, 88, 99), NA, happy),
    imprich = ifelse(imprich %in% c(77, 88, 99), NA, imprich),
  )
table(uk_clean$EducationLevel)
## 
##    BA No BA 
##  6021 14958
table(uk_clean$happy)
## 
## Happiness   Unhappy 
##     17741      3205
table(uk_clean$imprich)
## 
##          Important to Me Not that Important to Me 
##                     4780                    15713

#FINAL Tables

#CLEANED
summary(uk_happy)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   7.000   8.000   7.472   9.000  10.000      33
#CLEANED
summary(uk_imprich)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.000   4.000   5.000   4.327   5.000   6.000     486
table1_cat <- datasummary_skim(uk_clean  %>% dplyr::select(happy, imprich, EducationLevel), type = "categorical", title = "Table 1. Categorical Variables Summary", output = "flextable")

table1_cat
Table 1. Categorical Variables Summary

N

%

happy

Happiness

17741

84.6

Unhappy

3205

15.3

NA

33

0.2

imprich

Important to Me

4780

22.8

Not that Important to Me

15713

74.9

NA

486

2.3

EducationLevel

BA

6021

28.7

No BA

14958

71.3

#Descriptive Statistics#

ggplot(uk_clean, aes(x = EducationLevel, y = uk_happy, fill = factor(EducationLevel))) +
  geom_boxplot() +
  labs(title = "Happiness Scores According to Educational Attainment",
       x = "Educational Attainment",
       y = "Happiness Scores") +
  scale_fill_viridis_d(name = "EducationalLevel", option = "A") +
  theme_minimal()
## Warning: Removed 33 rows containing non-finite values (`stat_boxplot()`).

ggplot(uk_clean, aes(x = imprich, y = uk_happy, fill = factor(imprich))) +
  geom_boxplot() +
  labs(title = "Attitudes Towards Richness and Happiness Levels",
       x = "Richness Attitudes",
       y = "Happiness Scores") +
  scale_fill_viridis_d(name = "EducationalLevel", option = "A") +
  theme_minimal()
## Warning: Removed 33 rows containing non-finite values (`stat_boxplot()`).

#CHI SQUARED#

tab_dat <- uk_clean
mytab <- table(tab_dat$happy, tab_dat$EducationLevel)
rsums <- rowSums(mytab)
csums <- colSums(mytab)
N <- sum(mytab)

ptab <- tcrossprod(rsums/N, csums/N)
cat("Table of Expected Proportions")
## Table of Expected Proportions
print(ptab)
##            [,1]      [,2]
## [1,] 0.24314598 0.6038415
## [2,] 0.04392553 0.1090870
ftab <- N * ptab
cat("Table of Expected Frequencies:\n")
## Table of Expected Frequencies:
print(ftab)
##           [,1]      [,2]
## [1,] 5092.9358 12648.064
## [2,]  920.0642  2284.936
alpha <- 0.05 
c_val <- qchisq(alpha, df = 1, lower.tail = FALSE)
cat("Critical Value", round(c_val, 3), "\n")
## Critical Value 3.841
test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's χ2:", round(test_stat, 4), "\n")
## Pearson's χ2: 143.2039
p_val <- pchisq(test_stat, df = 1, lower.tail = F)
cat("p-value:", round(p_val, 4), "\n")
## p-value: 0
print(chisq.test(mytab, correct = F))
## 
##  Pearson's Chi-squared test
## 
## data:  mytab
## X-squared = 143.2, df = 1, p-value < 2.2e-16
alpha <- 0.05 
c_val <- qchisq(alpha, df = 1, lower.tail = FALSE)
cat("Critical Value", round(c_val, 3), "\n")
## Critical Value 3.841

#CHI-SQUARED FOR IMPRICH

tab_dat_imprich <- uk_clean
mytab_imprich <- table(tab_dat_imprich$happy, tab_dat_imprich$imprich)
rsums_imprich <- rowSums(mytab_imprich)
csums_imprich <- colSums(mytab_imprich)
N_imprich <- sum(mytab_imprich)

ptab_imprich <- tcrossprod(rsums_imprich/N, csums_imprich/N)
cat("Table of Expected Proportions")
## Table of Expected Proportions
print(ptab_imprich)
##            [,1]      [,2]
## [1,] 0.18901487 0.6211167
## [2,] 0.03372817 0.1108332
ftab_imprich <- N * ptab_imprich
cat("Table of Expected Frequencies (imprich):\n")
## Table of Expected Frequencies (imprich):
print(ftab_imprich)
##           [,1]      [,2]
## [1,] 3959.1056 13009.911
## [2,]  706.4702  2321.513
test_stat_imprich <- sum((mytab_imprich - ftab_imprich)^2 / ftab_imprich)
cat("Pearson's χ2 (imprich):", round(test_stat_imprich, 4), "\n")
## Pearson's χ2 (imprich): 16.0043
p_val_imprich <- pchisq(test_stat_imprich, df = 1, lower.tail = F)
cat("p-value:", round(p_val_imprich, 4), "\n")
## p-value: 1e-04
print(chisq.test(mytab_imprich, correct = F))
## 
##  Pearson's Chi-squared test
## 
## data:  mytab_imprich
## X-squared = 4.8899, df = 1, p-value = 0.02701

#Hypothesis Testing#

# Check if uk_clean is loaded and contains the necessary variables
if (!"uk_clean" %in% ls() || !all(c("EducationLevel", "happy") %in% colnames(uk_clean))) {
  stop("The dataset 'uk_clean' is not loaded or does not contain necessary variables.")
}

# Load necessary packages
library(infer)
library(dplyr)

test_stat <- uk_clean %>%
  specify(explanatory = EducationLevel, response = happy, success = "Happiness") %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "Chisq")
## Warning: Removed 33 rows containing missing values.
#Simulate Null Distribution
null_distribution <- test_stat
#Calculate p-value of sample
p_val <- null_distribution %>%
  get_pvalue(obs_stat = test_stat, direction = "greater")
## Warning: The first row and first column value of the given `obs_stat` will be
## used.
p_val
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1    0.28
null_distribution %>%
  visualize() +
  shade_p_value(obs_stat = test_stat, direction = "greater")
## Warning: The first row and first column value of the given `obs_stat` will be
## used.

conf_int <- null_distribution %>%
  get_confidence_interval(level = 0.95, type = "percentile")


null_distribution %>%
  visualize() +
  shade_p_value(obs_stat = test_stat, direction = "greater") +
  shade_confidence_interval(endpoints = conf_int)
## Warning: The first row and first column value of the given `obs_stat` will be
## used.

#Testing Imprich#

#Testing Imprich

test_stat_rich <- uk_clean %>%
  specify(explanatory = imprich, response = happy, success = "Happiness") %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "Chisq")
## Warning: Removed 513 rows containing missing values.
#Simulate Null Distribution
null_distribution_rich <- test_stat_rich
#Calculate p-value of sample
p_val <- null_distribution_rich %>%
  get_pvalue(obs_stat = test_stat_rich, direction = "greater")
## Warning: The first row and first column value of the given `obs_stat` will be
## used.
p_val
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1   0.307
null_distribution_rich %>%
  visualize() +
  shade_p_value(obs_stat = test_stat_rich, direction = "greater")
## Warning: The first row and first column value of the given `obs_stat` will be
## used.

conf_int <- null_distribution_rich %>%
  get_confidence_interval(level = 0.95, type = "percentile")


null_distribution_rich %>%
  visualize() +
  shade_p_value(obs_stat = test_stat, direction = "greater") +
  shade_confidence_interval(endpoints = conf_int)
## Warning: The first row and first column value of the given `obs_stat` will be
## used.

#Evidence to retain the Null for both of them

#Creating Models

#Re-scaling uk_happy variable from nominal (happy/unhappy) to happiness score (0-10)

table(uk_data$happy)
## 
##    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
uk_data$happy <- ifelse(uk_data$happy %in% c(77, 88, 99), NA, uk_data$happy)
table(uk_data$happy)
## 
##    0    1    2    3    4    5    6    7    8    9   10 
##  113  112  240  502  670 1568 1541 3804 6099 3740 2557
uk_data$happyscale <- scales::rescale(uk_data$happy, na.rm = TRUE, to = c(0,50))

ggplot(uk_data, aes(x = happy)) + geom_histogram(bins = 30, fill = "#FFD580", color = "black") + 
  theme_minimal() +
  labs(title = "Distribution of Happiness Scale for Hungary",
       x = "Happiness Scale",
       y = "Count")
## Warning: Removed 33 rows containing non-finite values (`stat_bin()`).

#Re-coding education variable to introduce in uk_data_scaled sections to sort all information for LRM Models into one place

uk_education_scaled <- uk_data$educ_attain
table(uk_data$educ_attain)
## 
##    BA No BA 
##  6021 14958
#Computing the Model
uk_happyBA_model <- lm(happyscale ~ educ_attain, data = uk_data)

#Identifying computed coefficients
coefficients_uk_happyBA <- coef(uk_happyBA_model)
print(coefficients_uk_happyBA)
##      (Intercept) educ_attainNo BA 
##        38.084151        -1.013904

#NUMERICAL DISPLAY OF BA/HAPPY MODEL

#Not necessarily needed, so voided from the project entirely, supposed to be “displaying without scientific notation” code from Tutorial 7

tidy_model <- tidy(uk_happyBA_model) knitr::kable(tidy_model, digits = 3)

#Summary of model

summary(uk_happyBA_model)
## 
## Call:
## lm(formula = happyscale ~ educ_attain, data = uk_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.084  -3.084   2.930   6.916  12.930 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       38.0842     0.1234 308.704  < 2e-16 ***
## educ_attainNo BA  -1.0139     0.1461  -6.939 4.05e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.566 on 20944 degrees of freedom
##   (33 observations deleted due to missingness)
## Multiple R-squared:  0.002294,   Adjusted R-squared:  0.002246 
## F-statistic: 48.15 on 1 and 20944 DF,  p-value: 4.055e-12

#FORMULATING IMPRICH MODEL

table(uk_data$imprich)
## 
##    1    2    3    4    5    6    7    8    9 
##  677 1631 2472 4040 8881 2792   47  136   64
uk_data <- uk_data %>%
  mutate(
    imprich = case_when(
      imprich %in% 1:3 ~ "Important to Me",
      imprich %in% 4:6 ~ "Not Important to Me",
      imprich %in% c(7, 8, 9) ~ NA_character_,
      TRUE ~ as.character(imprich)
    )
  )
table(uk_data$imprich)
## 
##     Important to Me Not Important to Me 
##                4780               15713
uk_happyRich_model <- lm(happyscale ~ imprich, data = uk_data)

coefficients_uk_happyRich <- coef(uk_happyRich_model)
print(coefficients_uk_happyRich)
##                (Intercept) imprichNot Important to Me 
##                 36.8586387                  0.6976037

#MODEL SUMMARY OF IMPRICH/HAPPY + BA/HAPPY MODELS

modelsummary(
  list(uk_happyRich_model),
  fmt = 1,
  estimate = c ("{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_omit = "Intercept",
  title = "Model Summary for BA/Happy Model"
)
Model Summary for BA/Happy Model
 (1)
imprichNot Important to Me 0.7 (0.2)***
Num.Obs. 20466
R2 0.001
R2 Adj. 0.001
AIC 150407.4
BIC 150431.1
Log.Lik. −75200.680
RMSE 9.54 
models <- list(
  "Happy/Education Model" = lm(happyscale ~ educ_attain, data = uk_data),
  "Happy/Rich Model" = lm(happyscale ~ imprich, data = uk_data)
)

modelsummary(models, fmt = 1,
  estimate  = c("{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_rename = c("educ_attainNo BA" = "No BA", "imprichNot Important to Me" = "Richness not Important to Me"),
  title = "Table 2. Regression Models for happiness scales"
)
Table 2. Regression Models for happiness scales
Happy/Education Model Happy/Rich Model
(Intercept) 38.1 (0.1)*** 36.9 (0.1)***
No BA −1.0 (0.1)***
Richness not Important to Me 0.7 (0.2)***
Num.Obs. 20946 20466
R2 0.002 0.001
R2 Adj. 0.002 0.001
AIC 154049.0 150407.4
BIC 154072.8 150431.1
Log.Lik. −77021.490 −75200.680
RMSE 9.57 9.54

#SETTING UP MLR

#INTERACTION MODEL 

interactionModel <- lm(happyscale ~ imprich + educ_attain + imprich*educ_attain, data = uk_data)
models <- list(
  "Happy/Education Model" = lm(happyscale ~ educ_attain, data = uk_data),
  "Happy/Rich Model" = lm(happyscale ~ imprich, data = uk_data),
  "Interaction Model" = lm(happyscale ~ imprich + educ_attain + imprich*educ_attain, data = uk_data)
)

modelsummary(models, fmt = 1,
  estimate  = c("{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}",
                "{estimate} ({std.error}){stars}"),
  statistic = NULL,
  coef_rename = c("educ_attainNo BA" = "No BA", "imprichNot Important to Me" = "Richness not Important to Me", "imprichNot Important to Me:educ_attainNo BA" = "No BA x Richness not Important"),
  title = "Table 2. Regression Models for Happiness Scales"
)
Table 2. Regression Models for Happiness Scales
Happy/Education Model Happy/Rich Model  Interaction Model
(Intercept) 38.1 (0.1)*** 36.9 (0.1)*** 37.5 (0.3)***
No BA −1.0 (0.1)*** −0.9 (0.3)**
Richness not Important to Me 0.7 (0.2)*** 0.8 (0.3)**
No BA x Richness not Important −0.1 (0.3)
Num.Obs. 20946 20466 20466
R2 0.002 0.001 0.003
R2 Adj. 0.002 0.001 0.003
AIC 154049.0 150407.4 150361.0
BIC 154072.8 150431.1 150400.6
Log.Lik. −77021.490 −75200.680 −75175.488
RMSE 9.57 9.54 9.53 
interaction_plot <- effect("imprich*educ_attain", interactionModel, na.rm = TRUE)

plot(interaction_plot,
     main = "Interaction Plot",
     xlab = "Richness Attitudes",
     ylab = "Happiness Scale")