## ── 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.5.0 ✔ 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
## Version 2.0.0 of `modelsummary`, to be released soon, will introduce a
## breaking change: The default table-drawing package will be `tinytable`
## instead of `kableExtra`. All currently supported table-drawing packages
## will continue to be supported for the foreseeable future, including
## `kableExtra`, `gt`, `huxtable`, `flextable, and `DT`.
##
## You can always call the `config_modelsummary()` function to change the
## default table-drawing package in persistent fashion. To try `tinytable`
## now:
##
## config_modelsummary(factory_default = 'tinytable')
##
## To set the default back to `kableExtra`:
##
## config_modelsummary(factory_default = 'kableExtra')
##
## 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
#Descriptive statistics with datasummary_skim
table_summary <- datasummary_skim(ireland_data[c("gender_numeric", "happy_numeric", "educ.ba_numeric")], title="Table 1: Descriptive Statistics for Gender, Educ.ba, and Happy")
table_summary
Table 1: Descriptive Statistics for Gender, Educ.ba, and Happy
Unique
Missing Pct.
Mean
SD
Min
Median
Max
gender_numeric
2
0
1.5
0.5
1.0
2.0
2.0
happy_numeric
11
0
7.5
1.9
0.0
8.0
10.0
educ.ba_numeric
2
0
0.3
0.4
0.0
0.0
1.0
# For educ.ba
#Checking variable types to determine test statistic
is.numeric(ireland_data$educ.ba)
## [1] FALSE
is.numeric(ireland_data$happy)
## [1] TRUE
# Step 1: Calculate the test statistic of your sample
test_stat <- ireland_data %>%
specify(explanatory = educ.ba,
response = happy) %>%
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 "BA or more" - "No BA", or divided in the order "BA or more" / "No BA"
## for ratio-based statistics. To specify this order yourself, supply `order =
## c("BA or more", "No BA")` to the calculate() function.
## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "BA or more" - "No BA", or divided in the order "BA or more" / "No BA"
## for ratio-based statistics. To specify this order yourself, supply `order =
## c("BA or more", "No BA")` 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.
## Warning in (function (mapping = NULL, data = NULL, stat = "identity", position = "identity", : All aesthetics have length 1, but the data has 1000 rows.
## ℹ Did you mean to use `annotate()`?
## Warning in (function (mapping = NULL, data = NULL, stat = "identity", position = "identity", : All aesthetics have length 1, but the data has 1000 rows.
## ℹ Did you mean to use `annotate()`?
#Regression Modelling
#Model 1: educ.ba main single predictor
#Model 2: multivariate linear regression educ.ba + gender
#Model 3: interaction model
#Making 'No BA' the reference category
ireland_data<-ireland_data %>%
mutate(
educ.ba = factor(educ.ba, levels = c("No BA", "BA or more")))
# Running Regression Models
model1 <-lm(happy ~ educ.ba, data = ireland_data)
model2 <-lm(happy ~ educ.ba + gender, data = ireland_data, weights=weight)
model3 <-lm(happy ~ educ.ba + gender + educ.ba*gender, data = ireland_data, 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,
coef_rename = c("educ.baBA or more" = "BA or more", "genderMale" = "Male"),
title = 'Regression model for level of education and gender on happiness for Ireland')
Regression model for level of education and gender on happiness for
Ireland
 (1)
  (2)
  (3)
(Intercept)
7.4 (0.0)***
7.0 (0.1)***
7.0 (0.1)***
BA or more
0.3 (0.0)***
0.3 (0.1)**
0.3 (0.1)*
Male
−0.3 (0.1)***
−0.3 (0.1)**
BA or more:Male
0.0 (0.2)
Num.Obs.
22110
2386
2386
R2
0.006
0.009
0.009
R2 Adj.
0.006
0.008
0.008
AIC
90438.2
11229.8
11231.7
BIC
90462.2
11252.9
11260.6
Log.Lik.
−45216.111
−5610.878
−5610.873
RMSE
1.87
1.95
1.95
# 3rd visual: Interaction plot for model 3
interaction_plot <- effect("educ.ba*gender", model3, na.rm=TRUE)
plot(interaction_plot,
main="Interaction Effect: Education Level & Gender on Happiness in Ireland",
xlab="Level of Education",
ylab="Happiness Levels")
interaction_plot
##
## educ.ba*gender effect
## gender
## educ.ba Female Male
## No BA 6.965602 6.662492
## BA or more 7.249267 6.926717
