Research Report SOC222
2024-04-10
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.
Setting up
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "MASS", "aod", "interactions", "kableExtra", "performance", "flextable", "scales")
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: '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] "performance" "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] "flextable" "performance" "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"
##
## [[13]]
## [1] "scales" "flextable" "performance" "kableExtra" "interactions"
## [6] "aod" "MASS" "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"uk_data <- read.fst("uk_data.fst")##Cleaning Predictor Variables
UK_clean_test <- uk_data %>% # put all into new cleaned variable
mutate(
health = ifelse(health %in% c(7, 8, 9), NA, health), # For 'happy', set values 7, 8, and 9 to NA.
happy = ifelse(happy %in% c(77, 88, 99), NA, happy) # For 'happy', set values 77, 88, and 99 to NA.
)# Populist scale (Norris/Inglehart use low levels of trust in parliaments, politicians and parties as an indicator of populism)
# Setting values greater than 10 to NA
UK_clean_test$trstplt <- ifelse(UK_clean_test$trstplt > 10, NA, UK_clean_test$trstplt)
UK_clean_test$trstprl <- ifelse(UK_clean_test$trstprl > 10, NA, UK_clean_test$trstprl)
UK_clean_test$trstprt <- ifelse(UK_clean_test$trstprt > 10, NA, UK_clean_test$trstprt)
# Creating and rescaling the trust variable
UK_clean_test$trust <- scales::rescale(UK_clean_test$trstplt + UK_clean_test$trstprl + UK_clean_test$trstprt, na.rm = TRUE, to = c(0, 100))
# Fliping scale
UK_clean_test$populist <- scales::rescale(UK_clean_test$trust, na.rm = TRUE, to=c(100,0))UK_clean_test <- UK_clean_test %>% filter(!is.na(populist), !is.na(health), !is.na(happy))table(UK_clean_test$happy) #Checking for NA##
## 0 1 2 3 4 5 6 7 8 9 10
## 101 95 219 441 594 1333 1335 3325 5420 3262 2167table(UK_clean_test$health) #Checking for NA##
## 1 2 3 4 5
## 5281 7513 3992 1182 324table_summary <- datasummary_skim(UK_clean_test %>% dplyr::select(populist, health, happy))
table_summary #Table with distributions (not shown in report)| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max | ||
|---|---|---|---|---|---|---|---|---|
| populist | 31 | 0 | 62.5 | 21.1 | 0.0 | 60.0 | 100.0 | |
| health | 5 | 0 | 2.1 | 1.0 | 1.0 | 2.0 | 5.0 | |
| happy | 11 | 0 | 7.5 | 1.9 | 0.0 | 8.0 | 10.0 |
Notes on distribution in code instead:
For the populist variable: the distribution is close to a normal distribution. The majority of people are found in the middle of the populist scale, while there is an increase in people on the higher ends of the populist scale
For the health variable: the distribution is right skewed, as the tail is stretching out to the right. This means the mean is higher than the median, and the majority of people have good health, compared to the small percentage that has bad health
For the happy variable: the distribution is left skewed, as the tail is stretching to the left. This means the median is larger than the mean, and the majority of people are on the happier end (7,8,9,10) than the lower end.
##The table that will go in my report, without the distributions
table1 <- datasummary_skim(UK_clean_test %>% dplyr::select(populist, health, happy), title = "Table 1. Descriptive statistics for main variables", output = "flextable")## Warning: The histogram argument is only supported for (a) output types "default",
## "html", "kableExtra", or "gt"; (b) writing to file paths with extensions
## ".html", ".jpg", or ".png"; and (c) Rmarkdown, knitr or Quarto documents
## compiled to PDF (via kableExtra) or HTML (via kableExtra or gt). Use
## `histogram=FALSE` to silence this warning.table1
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
populist | 31 | 0 | 62.5 | 21.1 | 0.0 | 60.0 | 100.0 |
health | 5 | 0 | 2.1 | 1.0 | 1.0 | 2.0 | 5.0 |
happy | 11 | 0 | 7.5 | 1.9 | 0.0 | 8.0 | 10.0 |
flextable::save_as_docx(table1, path = "Table1.docx",
width = 7.0, height = 7.0)
set_flextable_defaults(fonts_ignore=TRUE)
print(table1, preview = "docx") # Saving Table## a flextable object.
## col_keys: ` `, `Unique (#)`, `Missing (%)`, `Mean`, `SD`, `Min`, `Median`, `Max`
## header has 1 row(s)
## body has 3 row(s)
## original dataset sample:
## Unique (#) Missing (%) Mean SD Min Median Max
## 1 populist 31 0 62.5 21.1 0.0 60.0 100.0
## 2 health 5 0 2.1 1.0 1.0 2.0 5.0
## 3 happy 11 0 7.5 1.9 0.0 8.0 10.0//I am aware that it would be better to rescale it so the higher numbers mean better health like the happiness one, here is my work trying to complete this. However, Professor said it was not nesscessary.
UK_clean_test\(health1 <- scales::rescale(UK_clean_test\)health, na.rm = TRUE)
UK_clean_test\(healthrescale <- scales::rescale(UK_clean_test\)health1, na.rm = TRUE)
Calculate Test Statistic
test_stat <- UK_clean_test %>%
specify(explanatory = health, # explanatory variable
response = populist) %>% # outcome/response variable
hypothesize(null = "independence") %>% #optional for test statistic / can ignore
calculate(stat = "slope") # picked because both health and populist variables are numeric, therefore a linear regression needs to be run (there is only one numeric variable for response variable and one for the explanatory variable). In the visual when a 'slope' is visualized in a linear regression, I will be able to see how populist attitudes change as health increases towards 5.## Message: The independence null hypothesis does not inform calculation of the
## observed statistic (a slope) and will be ignored.print(test_stat$stat)## health
## 3.285947Simulate Null distrubution
Figure1_null_distribution <- UK_clean_test %>%
specify(explanatory = health,
response = populist) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "slope")P value of sample
p_val <- Figure1_null_distribution %>% # calculating p value
get_pvalue(obs_stat = test_stat, direction = "two-sided") # two sided is usually the direction for lienar regressions## 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##Visualization of test statisitc of sample with null distrubution
Figure1_null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided")
Figure1_null_distribution # visualization of null distribution without confidence interval ## Response: populist (numeric)
## Explanatory: health (numeric)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.0348
## 2 2 0.232
## 3 3 -0.257
## 4 4 -0.182
## 5 5 -0.161
## 6 6 0.0577
## 7 7 -0.133
## 8 8 0.211
## 9 9 -0.102
## 10 10 -0.0532
## # ℹ 990 more rowsconf_int <- Figure1_null_distribution %>%
get_confidence_interval(level = 0.95, type = "percentile") #adding confidence interval of 95%
Figure1_null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int) +
labs(title = "Figure 1. Simulation-Based Null Distribution",
x = "slope",
y = "count")
Figure1_null_distribution # null distribution with confidence interval## Response: populist (numeric)
## Explanatory: health (numeric)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.0348
## 2 2 0.232
## 3 3 -0.257
## 4 4 -0.182
## 5 5 -0.161
## 6 6 0.0577
## 7 7 -0.133
## 8 8 0.211
## 9 9 -0.102
## 10 10 -0.0532
## # ℹ 990 more rows//Showing work
ggsave(filename = “Figure2.pdf”, plot = Figure1_null_distribution, device = pdf, width = 10, height = 10) //My attempt to save Figure 1.
Figure2 <- datasummary_skim(Figure1_null_distribution, %>% dplyr::select(populist, health, happy), title = “Figure 1. Null Distribution”, output = “flextable”)
//I have to resort to a screenshot. I tried different ways of saving this figure. However, the only ‘save’ code that I was able to find was how to save the previous “table 1”
Modelling Output
UK_clean_test$weight <- UK_clean_test$dweight * UK_clean_test$pweight
survey_design <- svydesign(ids = ~1, data = UK_clean_test, weights = ~weight) # Weighing the data#putting SLR, MLR and MLR + interaction in model variables
model1 <- lm(populist ~ health, data = UK_clean_test, weights = weight)
model2 <- lm(populist ~ health + happy, data = UK_clean_test, weights = weight)
model3 <- lm(populist ~ health + happy + health*happy, data = UK_clean_test, weights = weight) #Putting the 3 models in a list under one variable name
models <- list(
"Model 1" = lm(populist ~ health, data = UK_clean_test, weights = weight),
"Model 2" = lm(populist ~ health + happy, data = UK_clean_test, weights = weight),
"Model 3" = lm(populist ~ health + happy + health*happy, data = UK_clean_test, weights = weight))# Regression models table
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL,
title = 'Figure 2. Regression models for populist attitudes') # Title for the table (variables do not need to be renamed here)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 54.6 (0.4)*** | 72.0 (0.8)*** | 71.4 (1.6)*** |
| health | 3.5 (0.2)*** | 2.4 (0.2)*** | 2.7 (0.6)*** |
| happy | −2.0 (0.1)*** | −1.9 (0.2)*** | |
| health × happy | 0.0 (0.1) | ||
| Num.Obs. | 18292 | 18292 | 18292 |
| R2 | 0.024 | 0.054 | 0.054 |
| R2 Adj. | 0.024 | 0.053 | 0.053 |
| AIC | 165278.2 | 164722.4 | 164724.2 |
| BIC | 165301.6 | 164753.7 | 164763.3 |
| Log.Lik. | −82636.093 | −82357.223 | −82357.110 |
| RMSE | 20.85 | 20.52 | 20.51 |
**before weighting the data, the interaction effect acts according to my hypothesis and I get different coefficients. To my knowledge, weighting the data makes it more accurate. If I was not supposed to weigh the data, I am aware it would change the coefficients.
check_model(model3) //Although not required, I wanted to check the model fit. Gives error “package ‘see’ required for model diagnostic plots”. I tried adding the performance package at the start.
#Visual of Choice
interaction_plot <- effect("health*happy", model3, na.rm=TRUE) # Choosing the interaction plot because I want to showcase in a visual example the effects happiness also has on populist attitudes. I find it interesting that the y intercepts gets lower as happiness goes up
plot(interaction_plot,
main="Figure 3. Interaction effect",
xlab="Health",
ylab="Populist Attitudes Scale")
interaction_plot##
## health*happy effect
## happy
## health 0 2 5 8 10
## 1 74.10048 70.17833 64.29510 58.41188 54.48973
## 2 76.80590 72.81242 66.82220 60.83198 56.83850
## 3 79.51132 75.44651 69.34930 63.25208 59.18728
## 4 82.21674 78.08060 71.87639 65.67219 61.53605
## 5 84.92216 80.71469 74.40349 68.09229 63.88482/I tried saving saving interaction plot. I do not see code in the tutorial that enables me to do this
ggsave(filename = “Figure3.pdf”, plot = interaction_plot, device = pdf, width = 10, height = 10)