Final Project
2024-04-11
Lauren Petsche and Jack Zhang
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "performance", "broom", "MASS", "aod", "interactions", "kableExtra", "flextable", "scales", "dplyr")
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.5
## ✔ 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] "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] "broom" "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] "MASS" "broom" "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] "aod" "MASS" "broom" "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] "interactions" "aod" "MASS" "broom" "performance"
## [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] "kableExtra" "interactions" "aod" "MASS" "broom"
## [6] "performance" "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"
##
## [[13]]
## [1] "flextable" "kableExtra" "interactions" "aod" "MASS"
## [6] "broom" "performance" "survey" "survival" "Matrix"
## [11] "grid" "effects" "carData" "modelsummary" "fst"
## [16] "infer" "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [26] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [31] "datasets" "methods" "base"
##
## [[14]]
## [1] "scales" "flextable" "kableExtra" "interactions" "aod"
## [6] "MASS" "broom" "performance" "survey" "survival"
## [11] "Matrix" "grid" "effects" "carData" "modelsummary"
## [16] "fst" "infer" "lubridate" "forcats" "stringr"
## [21] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [26] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [31] "utils" "datasets" "methods" "base"
##
## [[15]]
## [1] "scales" "flextable" "kableExtra" "interactions" "aod"
## [6] "MASS" "broom" "performance" "survey" "survival"
## [11] "Matrix" "grid" "effects" "carData" "modelsummary"
## [16] "fst" "infer" "lubridate" "forcats" "stringr"
## [21] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [26] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [31] "utils" "datasets" "methods" "base"loading data
norway_data <- read_fst("norway_data.fst")df <- norway_datadf$weight <- df$dweight * df$pweight
survey_design <- svydesign(ids = ~1, data = df, weights = ~weight)df$happy <- ifelse(df$happy > 10, NA, df$yrbrn)df <- df %>%
mutate(
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")),
age = ifelse(agea == 999, NA, agea),
cohort = ifelse(yrbrn < 1930 | yrbrn > 2000, NA, yrbrn),
gen = case_when(
yrbrn %in% 1900:1945 ~ "1",
yrbrn %in% 1946:1964 ~ "2",
yrbrn %in% 1965:1979 ~ "3",
yrbrn %in% 1980:1996 ~ "4",
TRUE ~ as.character(cohort)
),
gen = factor(gen,
levels = c("1", "2", "3", "4"),
labels = c("Interwar", "Baby Boomers", "Gen X", "Millennials"))
)
table(df$gen)##
## Interwar Baby Boomers Gen X Millennials
## 2752 5177 4324 3240df <- df %>%
mutate(
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"))
) regression models
model1 <- lm(happy ~ educ.ba, data = df, weights = weight)
model2 <- lm(happy ~ educ.ba + gen, data = df, weights = weight)
model3 <- lm(happy ~ educ.ba + gen + educ.ba*gen, data = df, weights = weight)
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) | 1980.9 (3.6)*** | 1936.2 (0.1)*** | 1935.6 (0.1)*** |
| educ.baBA or more | 7.1 (6.0) | 0.0 (0.1) | 2.6 (0.2)*** |
| genBaby Boomers | 19.2 (0.1)*** | 19.7 (0.2)*** | |
| genGen X | 35.4 (0.1)*** | 35.7 (0.2)*** | |
| genMillennials | 51.1 (0.1)*** | 52.6 (0.2)*** | |
| educ.baBA or more × genBaby Boomers | −2.4 (0.3)*** | ||
| educ.baBA or more × genGen X | −2.1 (0.3)*** | ||
| educ.baBA or more × genMillennials | −5.3 (0.3)*** | ||
| Num.Obs. | 16043 | 15472 | 15472 |
| R2 | 0.000 | 0.907 | 0.909 |
| R2 Adj. | 0.000 | 0.907 | 0.909 |
| AIC | 235337.6 | 96795.5 | 96502.7 |
| BIC | 235360.7 | 96841.4 | 96571.5 |
| Log.Lik. | −117665.820 | −48391.770 | −48242.342 |
| RMSE | 329.24 | 5.53 | 5.49 |
pdf1 <- modelsummary(list(model1, model2, model3), fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
coef_rename = c("genBaby Boomers" = "Baby Boomers", "genGen X" = "Gen X", "genMillennials" = "Millennials", "educ.baBA or more" = "BA or More"),
title = 'Table 1. Regression Models for Happiness')
pdf1| (1) | (2) | (3) | |
|---|---|---|---|
| (Intercept) | 1980.9 (3.6)*** | 1936.2 (0.1)*** | 1935.6 (0.1)*** |
| BA or More | 7.1 (6.0) | 0.0 (0.1) | 2.6 (0.2)*** |
| Baby Boomers | 19.2 (0.1)*** | 19.7 (0.2)*** | |
| Gen X | 35.4 (0.1)*** | 35.7 (0.2)*** | |
| Millennials | 51.1 (0.1)*** | 52.6 (0.2)*** | |
| BA or More:Baby Boomers | −2.4 (0.3)*** | ||
| BA or More:Gen X | −2.1 (0.3)*** | ||
| BA or More:Millennials | −5.3 (0.3)*** | ||
| Num.Obs. | 16043 | 15472 | 15472 |
| R2 | 0.000 | 0.907 | 0.909 |
| R2 Adj. | 0.000 | 0.907 | 0.909 |
| AIC | 235337.6 | 96795.5 | 96502.7 |
| BIC | 235360.7 | 96841.4 | 96571.5 |
| Log.Lik. | −117665.820 | −48391.770 | −48242.342 |
| RMSE | 329.24 | 5.53 | 5.49 |
flextable::save_as_docx(pdf1, path = "table1final.docx",
width = 7.0, height = 7.0)
set_flextable_defaults(fonts_ignore=TRUE)
print(pdf1, preview = "pdf")## <table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
## <caption>Table 1. Regression Models for Happiness</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> </th>
## <th style="text-align:center;"> (1) </th>
## <th style="text-align:center;"> (2) </th>
## <th style="text-align:center;"> (3) </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> (Intercept) </td>
## <td style="text-align:center;"> 1980.9 (3.6)*** </td>
## <td style="text-align:center;"> 1936.2 (0.1)*** </td>
## <td style="text-align:center;"> 1935.6 (0.1)*** </td>
## </tr>
## <tr>
## <td style="text-align:left;"> BA or More </td>
## <td style="text-align:center;"> 7.1 (6.0) </td>
## <td style="text-align:center;"> 0.0 (0.1) </td>
## <td style="text-align:center;"> 2.6 (0.2)*** </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Baby Boomers </td>
## <td style="text-align:center;"> </td>
## <td style="text-align:center;"> 19.2 (0.1)*** </td>
## <td style="text-align:center;"> 19.7 (0.2)*** </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Gen X </td>
## <td style="text-align:center;"> </td>
## <td style="text-align:center;"> 35.4 (0.1)*** </td>
## <td style="text-align:center;"> 35.7 (0.2)*** </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Millennials </td>
## <td style="text-align:center;"> </td>
## <td style="text-align:center;"> 51.1 (0.1)*** </td>
## <td style="text-align:center;"> 52.6 (0.2)*** </td>
## </tr>
## <tr>
## <td style="text-align:left;"> BA or More:Baby Boomers </td>
## <td style="text-align:center;"> </td>
## <td style="text-align:center;"> </td>
## <td style="text-align:center;"> −2.4 (0.3)*** </td>
## </tr>
## <tr>
## <td style="text-align:left;"> BA or More:Gen X </td>
## <td style="text-align:center;"> </td>
## <td style="text-align:center;"> </td>
## <td style="text-align:center;"> −2.1 (0.3)*** </td>
## </tr>
## <tr>
## <td style="text-align:left;box-shadow: 0px 1.5px"> BA or More:Millennials </td>
## <td style="text-align:center;box-shadow: 0px 1.5px"> </td>
## <td style="text-align:center;box-shadow: 0px 1.5px"> </td>
## <td style="text-align:center;box-shadow: 0px 1.5px"> −5.3 (0.3)*** </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Num.Obs. </td>
## <td style="text-align:center;"> 16043 </td>
## <td style="text-align:center;"> 15472 </td>
## <td style="text-align:center;"> 15472 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> R2 </td>
## <td style="text-align:center;"> 0.000 </td>
## <td style="text-align:center;"> 0.907 </td>
## <td style="text-align:center;"> 0.909 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> R2 Adj. </td>
## <td style="text-align:center;"> 0.000 </td>
## <td style="text-align:center;"> 0.907 </td>
## <td style="text-align:center;"> 0.909 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> AIC </td>
## <td style="text-align:center;"> 235337.6 </td>
## <td style="text-align:center;"> 96795.5 </td>
## <td style="text-align:center;"> 96502.7 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> BIC </td>
## <td style="text-align:center;"> 235360.7 </td>
## <td style="text-align:center;"> 96841.4 </td>
## <td style="text-align:center;"> 96571.5 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Log.Lik. </td>
## <td style="text-align:center;"> −117665.820 </td>
## <td style="text-align:center;"> −48391.770 </td>
## <td style="text-align:center;"> −48242.342 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> RMSE </td>
## <td style="text-align:center;"> 329.24 </td>
## <td style="text-align:center;"> 5.53 </td>
## <td style="text-align:center;"> 5.49 </td>
## </tr>
## </tbody>
## </table>interaction plot
model3 <- lm(happy ~ + educ.ba + cohort + educ.ba*cohort, data = df, weights = weight)
plot3 <- interact_plot(model3, pred = cohort, modx = educ.ba, jnplot = TRUE)
plot3
ggsave(filename = "plot1.pdf", plot = plot3, device = "pdf", width = 10, height = 8)null distribution
norway_task1 <- norway_data %>%
mutate(
ID = case_when(
happy %in% 0:4 ~ "unhappy",
happy %in% 6:10 ~ "happy",
happy %in% c(77, 88) ~ NA_character_,
TRUE ~ 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"))
) %>%
filter(!is.na(ID), !is.na(educ.ba)) mytab <- table(norway_task1$ID, norway_task1$educ.ba)
mytab##
## No BA BA or more
## happy 9425 5384
## unhappy 359 145rsums <- rowSums(mytab)
csums <- colSums(mytab)
N <- sum(mytab)
expected <- outer(rsums, csums) / N
expected## No BA BA or more
## happy 9461.9771 5347.0229
## unhappy 322.0229 181.9771alpha <- 0.05
df5 <- (nrow(mytab) - 1) * (ncol(mytab) - 1)
critical_value <- qchisq(1 - alpha, df5)
observed <- mytab
chi_squared <- sum((observed - expected)^2 / expected)
p_value <- pchisq(chi_squared, df5, lower.tail = FALSE)
cat("Critical Value:", critical_value, "\n")## Critical Value: 3.841459cat("Chi-Squared Statistic:", chi_squared, "\n")## Chi-Squared Statistic: 12.15985cat("P-Value:", p_value, "\n")## P-Value: 0.0004882927edu_stat <- norway_task1 %>%
specify(explanatory = ID,
response = educ.ba,
success = "BA or more") %>%
hypothesize(null = "independence") %>%
calculate(stat = "chisq") ## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "happy" - "unhappy", or divided in the order "happy" / "unhappy" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("happy", "unhappy")` to the calculate() function.print(edu_stat$stat)## X-squared
## 11.83323null_edu <- norway_task1 %>%
specify(explanatory = ID,
response = educ.ba,
success = "BA or more") %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "chisq")## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "happy" - "unhappy", or divided in the order "happy" / "unhappy" for
## ratio-based statistics. To specify this order yourself, supply `order =
## c("happy", "unhappy")` to the calculate() function.null_edu## Response: educ.ba (factor)
## Explanatory: ID (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 1.78e- 1
## 2 2 1.38e+ 0
## 3 3 3.95e-28
## 4 4 2.41e+ 0
## 5 5 1.86e+ 0
## 6 6 8.06e- 1
## 7 7 2.06e- 2
## 8 8 1.08e- 1
## 9 9 6.46e- 1
## 10 10 5.46e- 2
## # ℹ 990 more rowsp_val <- null_edu %>%
get_pvalue(obs_stat = edu_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 0conf_int <- null_edu%>%
get_confidence_interval(level = 0.95, type = "percentile")
null_edu %>%
visualize(data, bins = 15, method = "simulation", dens_color = "black") +
shade_p_value(obs_stat = edu_stat, direction = "greater") +
shade_confidence_interval(endpoints = conf_int)## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
edu_stat <- 1000
conf_int <- null_edu %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_edu %>%
visualize(bins = 10, method = "simulation", dens_color = "black") +
geom_vline(aes(xintercept = edu_stat), color = "red", linetype = "dashed", size = 1) +
shade_p_value(obs_stat = edu_stat, direction = "greater") +
shade_confidence_interval(endpoints = conf_int) +
coord_cartesian(xlim = c(0, max(null_edu$stat) * 1.1)) +
labs(title = "Null Distribution for Norway",
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
In the visualization of the null distibution the red line indicates evidence under the hypothesis testing to reject the null.
Critical Value (3.841459):
Chi-Squared Statistic (12.15985):
The Chi squared statistic value is greater than the critical value which indicates a deviation from what would be expected is there was no association between the variables. This also indicates that there is some evidence to reject the null hypothesis.
P-Value (0.0004882927)
This P value would provide evidence to reject the null hypothesis. This value is smaller than the significance level of (α = 0.05) which provides evidence to reject the null and support the alternative hypothesis