Assignment
Nicole and Kylie
2024-04-10
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "performance", "flextable", "broom", "scales", "ggeffects", "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.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## Warning: package 'infer' was built under R version 4.3.3## Warning: package 'effects' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Warning: package 'survey' was built under R version 4.3.3## 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## Warning: package 'ggeffects' was built under R version 4.3.3## Warning: package 'marginaleffects' was built under R version 4.3.3## [[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] "ggeffects" "scales" "broom" "flextable" "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] "marginaleffects" "ggeffects" "scales" "broom"
## [5] "flextable" "performance" "survey" "survival"
## [9] "Matrix" "grid" "effects" "carData"
## [13] "modelsummary" "fst" "infer" "lubridate"
## [17] "forcats" "stringr" "dplyr" "purrr"
## [21] "readr" "tidyr" "tibble" "ggplot2"
## [25] "tidyverse" "stats" "graphics" "grDevices"
## [29] "utils" "datasets" "methods" "base"#loading in the country of interest based on the report in which this case is Sweden.
sweden_data <-read.fst("sweden_data.fst")#main predictor in the research in our own interest will be income because it can determine someone’s attainment of success within education.
sweden_data <- sweden_data %>%
filter(cntry == "SE")
sweden_data$weight <- sweden_data$dweight * sweden_data$pweight
survey_design <- svydesign(ids = ~1, data = sweden_data, weights = ~weight)
sweden_data <- sweden_data %>%
mutate(income = case_when(
hincfel == 1 ~ 1, # the number one is in representation of higher income
hincfel == 2 ~ 0, # 0 setting to be a lower income
TRUE ~ NA_real_ # the rest of the data collected to NA
))
sweden_data <- sweden_data %>% filter(!is.na(income))The choice of predictor is education. Within the following code we insert what is the tabling sweden’s education. BA standing for Bachelor of arts, more than that would represent a higher education and No BA as lower education level.
swedendf <- sweden_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(swedendf$educ.ba)##
## No BA BA or more
## 11059 5406#second predictor being the focus area of study in Sweden. This can relate back to the main focus of income and why in such the areas of jobs determine certain successes within income in the country of Sweden. We are focusing on higher titles as they are more common to bring more wealth.
sweden_data <- swedendf %>%
mutate(field_of_study = recode(as.character(edufld),
'1' = "general no specific field", #some people may have jobs that required no degree
'2' = "law ", #law field jobs
'3' = "technical and engineering", #technical and engineering jobs or not even jobs but degrees obtained
'4' = "humanities",
'5' = "medical/healthservices/nursing/etc.",
'66' = NA_character_,
'99' = NA_character_,))
#the study fields do not just classify jobs but the specific categorization of future potential oppurtunites.
# check
table(sweden_data$field_of_study)##
## 10 11
## 377 51
## 12 13
## 170 54
## 14 6
## 114 325
## 7 77
## 278 3
## 8 88
## 609 16
## 9 general no specific field
## 706 897
## humanities law
## 1047 107
## medical/healthservices/nursing/etc. technical and engineering
## 125 159#let’s take a look at Sweden by itself, arranging the data to create a graph to illustrate our results.
df <- sweden_data
df$year <- NA
replacements <- c(2001, 2003, 2005, 2007, 2009, 2011, 2013, 2015, 2017, 2019) #setting the years to odd numbers to determine the income.
for(i in 1:10){
df$year[df$essround == i] <- replacements[i]
}
income_avg <- df %>%
group_by(year) %>%
summarize(income_avg = mean(income, na.rm = TRUE))
income_plot_SE <- ggplot(income_avg, aes(x = year, y = income_avg)) +
geom_point(aes(size = income_avg, color = income_avg), alpha = 0.7) +
geom_line(aes(group = 1), color = "green", linetype = "dashed") +
labs(
x = "Survey Year",
y = "In(Income)") +
theme_minimal() +
theme(legend.position = "none")
# Y axis has income and x axis has the year. Years are arranged by odd numbers. With the highest year being 2019
print(income_plot_SE)
#stepping into the categorical tables to code the same data related to
sweden.
swedendf <- swedendf %>%
mutate(income = case_when(
hincfel == 1 ~ "Yes", # 1 determines 'yes'
hincfel == 2 ~ "No", # 0 determines 'no'
TRUE ~ NA_character_ # Keep NA values without changing it
))
swedent1 <- swedendf %>%
# renanming the variables
mutate(
income = income,
education = educ.ba,
field = edufld
)
table1 <- datasummary_skim(swedent1 %>% dplyr::select(income, educ.ba, edufld), type = "categorical",
title = "Table 1. Descriptive statistics for main variables", output = "flextable")
table1
|
| N | % |
|---|---|---|---|
income | No | 5777 | 35.1 |
Yes | 10688 | 64.9 | |
educ.ba | No BA | 11059 | 67.2 |
BA or more | 5406 | 32.8 |
#stepping into the numerical tables. This is step 1 of the process.
test_stat <- sweden_data %>%
specify(explanatory = educ.ba, # changing the variable name for explanatory variable
response = income) %>% # changing the variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "t") # use of quotation marks to appropriate the test statistic## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "No BA" - "BA or more", or divided in the order "No BA" / "BA or more"
## for ratio-based statistics. To specify this order yourself, supply `order =
## c("No BA", "BA or more")` to the calculate() function.print(test_stat$stat)## t
## -15.33989the results to step one gave us the result of -15.3
stepping into step 2 of the process.
null_distribution <- sweden_data %>%
specify(explanatory = educ.ba,
response = income) %>%
hypothesize(null = "independence") %>%
generate(reps = 500, type = "permute") %>% # only line we are adding, use reps = 500 as standard but change if too computationally demanding. Added the 500 to not alterate or demand to much of the program as it crashes.
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 "No BA" - "BA or more", or divided in the order "No BA" / "BA or more"
## for ratio-based statistics. To specify this order yourself, supply `order =
## c("No BA", "BA or more")` to the calculate() function.we are at step 3
p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "two-sided") # would only replace test_stat if assigned another name in Step 1## 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 0the p value resulted in 0.
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided")
the results range around 0.
null_distribution## Response: income (numeric)
## Explanatory: educ.ba (factor)
## Null Hypothesis: independence
## # A tibble: 500 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.968
## 2 2 -0.305
## 3 3 0.563
## 4 4 0.251
## 5 5 0.0425
## 6 6 0.563
## 7 7 -0.305
## 8 8 -0.445
## 9 9 -0.201
## 10 10 -1.14
## # ℹ 490 more rowsconf_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)
the target predictor plus the 2 other predictors that help the main predictor.
model1 <- lm(income ~ educ.ba, data = df, weights = weight)
model2 <- lm(income ~ educ.ba + edufld, 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)| (1) | (2) | |
|---|---|---|
| (Intercept) | 0.6 (0.0)*** | 0.6 (0.0)*** |
| educ.baBA or more | 0.1 (0.0)*** | 0.1 (0.0)*** |
| edufld | 0.0 (0.0) | |
| Num.Obs. | 16465 | 5162 |
| R2 | 0.014 | 0.009 |
| R2 Adj. | 0.014 | 0.008 |
| AIC | 22183.0 | 7095.2 |
| BIC | 22206.2 | 7121.4 |
| Log.Lik. | −11088.524 | −3543.585 |
| F | 228.091 | |
| RMSE | 0.47 | 0.48 |
the visual of our own choice, with results in backup to the predictors chosen, which in this case are, education and education BA or none. We are also inserting gender as an extra for this graph, because gender can also be another predictor to our main idea of income and why and how genders are involved into the category.
sweden_data <- sweden_data %>%
mutate(
gndr = case_when(
gndr == 1 ~ "Male",
gndr == 2 ~ "Female",
TRUE ~ NA_character_
),
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
filter(!is.na(gndr) & !is.na(educ.ba))education_by_gender <- sweden_data %>%
group_by(gndr, educ.ba) %>%
summarise(count = n(), .groups = 'drop') %>%
mutate(percentage = count / sum(count) * 500)education_plot <- ggplot(education_by_gender, aes(x = educ.ba, y = count, fill = educ.ba)) +
geom_bar(stat = "identity", position = "dodge") +
facet_wrap(~gndr) +
labs(title = "Educational Attainment by Gender in sweden (2019)",
x = "Educational Level",
y = "Count of Respondents") +
scale_fill_manual(values = c("BA or more" ="Green", "No BA" = "Yellow")) +
theme_minimal()
education_plot
#Most accurate numbers and data collected can come from the website Link provided by the coding done throughout the semester the followling link allows us to see the data based on other european countries, link: https://www.europeansocialsurvey.org/methodology/ess-methodology/data-processing-and-archiving/weighting
#Hope the report was found to be interesting!