project
2024-03-31
Data packet loading
# List of packages
packages <- c("tidyverse", "fst", "modelsummary") # 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.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## [[1]]
## [1] "lubridate" "forcats" "stringr" "dplyr" "purrr" "readr"
## [7] "tidyr" "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "fst" "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [7] "readr" "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [13] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "modelsummary" "fst" "lubridate" "forcats" "stringr"
## [6] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [11] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [16] "utils" "datasets" "methods" "base"Data processing
belgium_data <- read.fst("belgium_data.fst")num_variables <- length(names(belgium_data))
num_variables## [1] 2665table(belgium_data$polintr)##
## 1 2 3 4 8 9
## 1645 6356 5931 3512 3 4p1,17451
belgium_data$polintr[belgium_data$polintr %in% c(8, 9)] <- NAtable(belgium_data$polintr)##
## 1 2 3 4
## 1645 6356 5931 3512sum(!is.na(belgium_data$polintr))## [1] 17444p2,17451
table(belgium_data$clsprty)##
## 1 2 7 8 9
## 8755 8549 23 115 9belgium_data$clsprty[belgium_data$clsprty %in% c(7, 8, 9)] <- NA
table(belgium_data$clsprty)##
## 1 2
## 8755 8549sum(!is.na(belgium_data$clsprty))## [1] 17304p3
In data processing, I initially selected three predictors, among which “ginfdif” stands for “Government should reduce differences in income levels,” as an alternative. If the first two predictors show independence with the outcome variable, then I would use the “ginfdif” variable for manipulation.
table(belgium_data$gincdif)##
## 1 2 3 4 5 7 8 9
## 4433 7740 2451 2191 523 6 101 6belgium_data$gincdif[belgium_data$gincdif %in% c(7, 8, 9)] <- NA
table(belgium_data$gincdif)##
## 1 2 3 4 5
## 4433 7740 2451 2191 523length(belgium_data$gincdif)## [1] 17451sum(!is.na(belgium_data$gincdif))## [1] 17338ov 17451
table(belgium_data$stfgov)##
## 0 1 2 3 4 5 6 7 8 9 10 77 88 99
## 987 593 1333 1960 2150 3596 3076 2381 829 146 92 5 296 7belgium_data$stfgov[belgium_data$stfgov %in% c(77, 88, 99)] <- NA
table(belgium_data$stfgov)##
## 0 1 2 3 4 5 6 7 8 9 10
## 987 593 1333 1960 2150 3596 3076 2381 829 146 92sum(!is.na(belgium_data$stfgov))## [1] 17143Data processing: Ensure that the predictor variables have the same number of valid observations
num_obs_polintr <- sum(!is.na(belgium_data$polintr))
num_obs_clsprty <- sum(!is.na(belgium_data$clsprty))
num_obs_gincdif <- sum(!is.na(belgium_data$gincdif))
num_obs_stfgov <- sum(!is.na(belgium_data$stfgov))num_obs_polintr## [1] 17444num_obs_clsprty## [1] 17304num_obs_gincdif## [1] 17338num_obs_stfgov## [1] 17143# Find the minimum number of observations among all variables
min_num_obs <- min(num_obs_polintr, num_obs_clsprty, num_obs_gincdif, num_obs_stfgov)
min_num_obs## [1] 17143# Retain the data with the minimum number of observations from each variable
belgium_data <- belgium_data[complete.cases(belgium_data$polintr, belgium_data$clsprty, belgium_data$gincdif, belgium_data$stfgov), ]num_obs_polintr <- sum(!is.na(belgium_data$polintr))
num_obs_clsprty <- sum(!is.na(belgium_data$clsprty))
num_obs_gincdif <- sum(!is.na(belgium_data$gincdif))
num_obs_stfgov <- sum(!is.na(belgium_data$stfgov))num_obs_polintr## [1] 16935num_obs_clsprty## [1] 16935num_obs_gincdif## [1] 16935num_obs_stfgov## [1] 16935Generate descriptive statistical charts.
library(dplyr)
# Mapping numerical values to corresponding text labels.
belgium_data <- belgium_data %>%
mutate(polintr = recode(polintr,
`1` = "Very interested",
`2` = "Quite interested",
`3` = "Hardly interested",
`4` = "Not at all interested"))table(belgium_data$polintr)##
## Hardly interested Not at all interested Quite interested
## 5762 3292 6258
## Very interested
## 1623library(dplyr)
# Mapping numerical values to corresponding text labels.
belgium_data <- belgium_data %>%
mutate(clsprty = recode(clsprty,
`1` = "Yes",
`2` = "No"))
table(belgium_data$clsprty)##
## No Yes
## 8309 8626table(belgium_data$stfgov)##
## 0 1 2 3 4 5 6 7 8 9 10
## 972 588 1322 1943 2133 3545 3052 2357 798 139 86Here, y remains numeric #### Plot analytical graph.
library(ggplot2)
# Draw bar plot for polintr
polintr_bar <- ggplot(belgium_data, aes(x = polintr, fill = polintr)) +
geom_bar() +
labs(title = "Interest in Politics",
x = "Interest Level",
y = "Count") +
scale_fill_brewer(palette = "Set3")
# Draw bar plot for clsprty
clsprty_bar <- ggplot(belgium_data, aes(x = clsprty, fill = clsprty)) +
geom_bar() +
labs(title = "Party Membership",
x = "Membership Status",
y = "Count") +
scale_fill_brewer(palette = "Set3")
# Draw histogram for stfgov
stfgov_hist <- ggplot(belgium_data, aes(x = stfgov)) +
geom_histogram(binwidth = 1, fill = "skyblue", color = "black") +
labs(title = "Appendix 1. Distribution of Satisfaction with Government",
x = "Satisfaction Level",
y = "Frequency")
# Display plots
polintr_bar
clsprty_bar
stfgov_hist
library(ggplot2)
# Calculate percentages for each category
polintr_percent <- prop.table(table(belgium_data$polintr)) * 100
clsprty_percent <- prop.table(table(belgium_data$clsprty)) * 100
stfgov_percent <- prop.table(table(belgium_data$stfgov)) * 100
# Create data frames
polintr_df <- data.frame(polintr = names(polintr_percent), percent = as.numeric(polintr_percent))
clsprty_df <- data.frame(clsprty = names(clsprty_percent), percent = as.numeric(clsprty_percent))
stfgov_df <- data.frame(stfgov = names(stfgov_percent), percent = as.numeric(stfgov_percent))
# Draw pie chart for polintr
polintr_pie <- ggplot(polintr_df, aes(x = "", y = percent, fill = polintr)) +
geom_bar(width = 1, stat = "identity") +
coord_polar(theta = "y") +
geom_text(aes(label = paste0(round(percent), "%")), position = position_stack(vjust = 0.5)) +
labs(title = "Interest in Politics") +
scale_fill_brewer(palette = "Set3") +
theme_void()
# Draw pie chart for clsprty
clsprty_pie <- ggplot(clsprty_df, aes(x = "", y = percent, fill = clsprty)) +
geom_bar(width = 1, stat = "identity") +
coord_polar(theta = "y") +
geom_text(aes(label = paste0(round(percent), "%")), position = position_stack(vjust = 0.5)) +
labs(title = "Party Membership") +
scale_fill_brewer(palette = "Set3") +
theme_void()
# Draw pie chart for stfgov
stfgov_pie <- ggplot(stfgov_df, aes(x = "", y = percent, fill = stfgov)) +
geom_bar(width = 1, stat = "identity") +
coord_polar(theta = "y") +
geom_text(aes(label = paste0(round(percent), "%")), position = position_stack(vjust = 0.5)) +
labs(title = "Satisfaction with Government") +
scale_fill_brewer(palette = "Set3") +
theme_void()
# Display pie charts
polintr_pie
clsprty_pie
stfgov_pie
table1 <- datasummary_skim(belgium_data %>% dplyr::select(polintr, clsprty), type = "categorical",
title = "Table 1. Descriptive statistics for main variables", output = "flextable")
table1
|
| N | % |
|---|---|---|---|
polintr | Hardly interested | 5762 | 34.0 |
Not at all interested | 3292 | 19.4 | |
Quite interested | 6258 | 37.0 | |
Very interested | 1623 | 9.6 | |
clsprty | No | 8309 | 49.1 |
Yes | 8626 | 50.9 |
table1b <- datasummary_skim(belgium_data %>% dplyr::select(stfgov), title = "Table 2. 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.table1b
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
stfgov | 11 | 0 | 4.6 | 2.1 | 0.0 | 5.0 | 10.0 |
Perform hypothesis test to determine independence
# List of packages
packages <- c("tidyverse", "infer", "fst") # 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)## Warning: 程辑包'infer'是用R版本4.3.3 来建造的## [[1]]
## [1] "fstcore" "modelsummary" "fst" "lubridate" "forcats"
## [6] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [11] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "infer" "fstcore" "modelsummary" "fst" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[3]]
## [1] "infer" "fstcore" "modelsummary" "fst" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"## Not necessary, but always good to rename when doing operations so you can backtrack to clean dataset
tab_dat <- belgium_data
# Create a crosstable and save it into the object mytab
mytab <- table(belgium_data$stfgov, belgium_data$polintr)
# Calculate row-wise sums of counts in table
rsums <- rowSums(mytab)
# Calculate column-wise sums of counts in table
csums <- colSums(mytab)
# Get the total count of the table (N)
N <- sum(mytab)
# Generate table of expected proportions
ptab <- tcrossprod(rsums/N, csums/N)
cat("Table of Expected Proportions:\n")## Table of Expected Proportions:print(ptab)## [,1] [,2] [,3] [,4]
## [1,] 0.019528510 0.0111572121 0.021209548 0.0055006547
## [2,] 0.011813543 0.0067494246 0.012830468 0.0033275565
## [3,] 0.026560381 0.0151747268 0.028846731 0.0074813431
## [4,] 0.039036930 0.0223029456 0.042397276 0.0109956503
## [5,] 0.042854231 0.0244838821 0.046543176 0.0120708811
## [6,] 0.071222808 0.0406916841 0.077353754 0.0200615441
## [7,] 0.061317915 0.0350327277 0.066596236 0.0172716030
## [8,] 0.047354629 0.0270550915 0.051430973 0.0133385217
## [9,] 0.016032666 0.0091599334 0.017412777 0.0045159696
## [10,] 0.002792657 0.0015955272 0.003033053 0.0007866163
## [11,] 0.001727831 0.0009871607 0.001876565 0.0004866834ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2] [,3] [,4]
## [1,] 330.71532 188.94739 359.18370 93.153587
## [2,] 200.06236 114.30151 217.28397 56.352170
## [3,] 449.80006 256.98400 488.51940 126.696546
## [4,] 661.09040 377.70038 717.99787 186.211337
## [5,] 725.73640 414.63454 788.20868 204.420372
## [6,] 1206.15825 689.11367 1309.98583 339.742250
## [7,] 1038.41890 593.27924 1127.80726 292.494597
## [8,] 801.95063 458.17797 870.98353 225.887865
## [9,] 271.51320 155.12347 294.88539 76.477945
## [10,] 47.29365 27.02025 51.36475 13.321346
## [11,] 29.26082 16.71757 31.77963 8.241984#(11 - 1) * (4 - 1) = 10 * 3 = 30
# Critical Value for Independence:
alpha <- 0.05 # significance level
c_val <- qchisq(alpha, df = 30, lower.tail = FALSE)
cat("Critical Value:", round(c_val, 3), "\n")## Critical Value: 43.773# Pearson's chi-squared statistic:
test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")## Pearson's X^2: 707.9037p_val <- pchisq(test_stat, df = 30, lower.tail = F)
cat("p-value:", round(p_val, 4), "\n")## p-value: 0## Not necessary, but always good to rename when doing operations so you can backtrack to clean dataset
tab_dat <- belgium_data
# Create a crosstable and save it into the object mytab
mytab <- table(belgium_data$stfgov, belgium_data$clsprty)
# Calculate row-wise sums of counts in table
rsums <- rowSums(mytab)
# Calculate column-wise sums of counts in table
csums <- colSums(mytab)
# Get the total count of the table (N)
N <- sum(mytab)
# Generate table of expected proportions
ptab <- tcrossprod(rsums/N, csums/N)
cat("Table of Expected Proportions:\n")## Table of Expected Proportions:print(ptab)## [,1] [,2]
## [1,] 0.028160776 0.029235149
## [2,] 0.017035531 0.017685461
## [3,] 0.038300973 0.039762209
## [4,] 0.056292581 0.058440221
## [5,] 0.061797259 0.064154911
## [6,] 0.102705712 0.106624079
## [7,] 0.088422520 0.091795963
## [8,] 0.068286985 0.070892229
## [9,] 0.023119650 0.024001697
## [10,] 0.004027107 0.004180747
## [11,] 0.002491591 0.002586649ftab <- N * ptab
cat("Table of Expected Frequencies:\n")## Table of Expected Frequencies:print(ftab)## [,1] [,2]
## [1,] 476.90275 495.09725
## [2,] 288.49672 299.50328
## [3,] 648.62699 673.37301
## [4,] 953.31485 989.68515
## [5,] 1046.53658 1086.46342
## [6,] 1739.32123 1805.67877
## [7,] 1497.43537 1554.56463
## [8,] 1156.44009 1200.55991
## [9,] 391.53127 406.46873
## [10,] 68.19906 70.80094
## [11,] 42.19510 43.80490#(11 - 1) * (2 - 1) = 10 * 1 = 10
# Critical Value for Independence:
alpha <- 0.05 # significance level
c_val <- qchisq(alpha, df = 10, lower.tail = FALSE)
cat("Critical Value:", round(c_val, 3), "\n")## Critical Value: 18.307# Pearson's chi-squared statistic:
test_stat <- sum((mytab - ftab)^2 / ftab)
cat("Pearson's X^2:", round(test_stat, 4), "\n")## Pearson's X^2: 135.4551p_val <- pchisq(test_stat, df = 10, lower.tail = F)
cat("p-value:", round(p_val, 4), "\n")## p-value: 0# 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)## Warning: 程辑包'effects'是用R版本4.3.3 来建造的## 载入需要的程辑包：carData## Warning: 程辑包'carData'是用R版本4.3.3 来建造的## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Warning: 程辑包'survey'是用R版本4.3.3 来建造的## 载入需要的程辑包：grid## 载入需要的程辑包：Matrix##
## 载入程辑包：'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## 载入需要的程辑包：survival##
## 载入程辑包：'survey'## The following object is masked from 'package:graphics':
##
## dotchart## Warning: 程辑包'performance'是用R版本4.3.3 来建造的## Warning: 程辑包'flextable'是用R版本4.3.3 来建造的##
## 载入程辑包：'flextable'## The following object is masked from 'package:purrr':
##
## compose##
## 载入程辑包：'scales'## The following object is masked from 'package:purrr':
##
## discard## The following object is masked from 'package:readr':
##
## col_factor## Warning: 程辑包'ggeffects'是用R版本4.3.3 来建造的## Warning: 程辑包'marginaleffects'是用R版本4.3.3 来建造的## [[1]]
## [1] "infer" "fstcore" "modelsummary" "fst" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[2]]
## [1] "infer" "fstcore" "modelsummary" "fst" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[3]]
## [1] "infer" "fstcore" "modelsummary" "fst" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[4]]
## [1] "infer" "fstcore" "modelsummary" "fst" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[5]]
## [1] "effects" "carData" "infer" "fstcore" "modelsummary"
## [6] "fst" "lubridate" "forcats" "stringr" "dplyr"
## [11] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [16] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [21] "datasets" "methods" "base"
##
## [[6]]
## [1] "survey" "survival" "Matrix" "grid" "effects"
## [6] "carData" "infer" "fstcore" "modelsummary" "fst"
## [11] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [16] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"
##
## [[7]]
## [1] "performance" "survey" "survival" "Matrix" "grid"
## [6] "effects" "carData" "infer" "fstcore" "modelsummary"
## [11] "fst" "lubridate" "forcats" "stringr" "dplyr"
## [16] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [21] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [26] "datasets" "methods" "base"
##
## [[8]]
## [1] "flextable" "performance" "survey" "survival" "Matrix"
## [6] "grid" "effects" "carData" "infer" "fstcore"
## [11] "modelsummary" "fst" "lubridate" "forcats" "stringr"
## [16] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [21] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [26] "utils" "datasets" "methods" "base"
##
## [[9]]
## [1] "broom" "flextable" "performance" "survey" "survival"
## [6] "Matrix" "grid" "effects" "carData" "infer"
## [11] "fstcore" "modelsummary" "fst" "lubridate" "forcats"
## [16] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [21] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [26] "grDevices" "utils" "datasets" "methods" "base"
##
## [[10]]
## [1] "scales" "broom" "flextable" "performance" "survey"
## [6] "survival" "Matrix" "grid" "effects" "carData"
## [11] "infer" "fstcore" "modelsummary" "fst" "lubridate"
## [16] "forcats" "stringr" "dplyr" "purrr" "readr"
## [21] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [26] "graphics" "grDevices" "utils" "datasets" "methods"
## [31] "base"
##
## [[11]]
## [1] "ggeffects" "scales" "broom" "flextable" "performance"
## [6] "survey" "survival" "Matrix" "grid" "effects"
## [11] "carData" "infer" "fstcore" "modelsummary" "fst"
## [16] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [21] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [26] "stats" "graphics" "grDevices" "utils" "datasets"
## [31] "methods" "base"
##
## [[12]]
## [1] "marginaleffects" "ggeffects" "scales" "broom"
## [5] "flextable" "performance" "survey" "survival"
## [9] "Matrix" "grid" "effects" "carData"
## [13] "infer" "fstcore" "modelsummary" "fst"
## [17] "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble"
## [25] "ggplot2" "tidyverse" "stats" "graphics"
## [29] "grDevices" "utils" "datasets" "methods"
## [33] "base"test_stat <- belgium_data %>%
specify(explanatory = polintr, # change variable name for explanatory variable
response = stfgov) %>% # change variable name for outcome of interest
hypothesize(null = "independence") %>%
calculate(stat = "f")
print(test_stat$stat)## [1] 111.9702null_dist <- belgium_data %>%
specify(response = stfgov, explanatory = polintr) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "f")
null_dist## Response: stfgov (numeric)
## Explanatory: polintr (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.121
## 2 2 1.15
## 3 3 0.171
## 4 4 0.628
## 5 5 0.395
## 6 6 0.442
## 7 7 0.191
## 8 8 0.953
## 9 9 0.0653
## 10 10 1.17
## # ℹ 990 more rowsp_val <- null_dist %>% #
get_pvalue(obs_stat = test_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_dist%>%
get_confidence_interval(level = 0.95, type = "percentile")
null_dist %>%
visualize(data, bins = 10, method = "simulation", dens_color = "black") +
shade_p_value(obs_stat = test_stat, direction = "greater") +
shade_confidence_interval(endpoints = conf_int)+
labs(title = "Figure 1. Simulation-Based Null Distribution polintr~stfgov")## Warning in min(diff(unique_loc)): min里所有的参数都不存在; 回覆Inf
test_stat <- 111.9702
conf_int <- null_dist %>%
get_confidence_interval(level = 0.95, type = "percentile")
# Visualize the null distribution with manual x-axis limits
null_dist %>%
visualize(bins = 15, method = "simulation", dens_color = "black") +
geom_vline(aes(xintercept = test_stat), color = "red", linetype = "dashed", size = 1) +
shade_p_value(obs_stat = test_stat, direction = "greater") +
shade_confidence_interval(endpoints = conf_int) +
coord_cartesian(xlim = c(0, max(null_dist$stat) * 1.1)) + # Set manual x-axis limits
labs(title = "Appendix 2. Simulation-Based Null Distribution polintr~stfgov",
x = "F",
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)): min里所有的参数都不存在; 回覆Inf
Enlarge the graphs for easier observation.
test_stat <- belgium_data %>%
specify(explanatory = clsprty,
response = stfgov) %>%
hypothesize(null = "independence") %>%
calculate(stat = "t", order = c("No", "Yes"))
print(test_stat$stat)## t
## -8.364563null_dist <- belgium_data %>%
specify(response = stfgov, explanatory = clsprty) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "t", order = c("No", "Yes"))
null_dist## Response: stfgov (numeric)
## Explanatory: clsprty (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.366
## 2 2 -1.62
## 3 3 0.632
## 4 4 0.243
## 5 5 -1.18
## 6 6 0.481
## 7 7 0.387
## 8 8 -0.0225
## 9 9 -2.00
## 10 10 -0.210
## # ℹ 990 more rows# Perform a t-test using the t.test() function and obtain the p-value.
t_result <- t.test(stfgov ~ clsprty, data = belgium_data)
p_value <- t_result$p.value
print(p_value)## [1] 6.499193e-17p_val <- null_dist %>% # replace name here if you assigned something other than null_distribution above
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 0null_dist %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided")
null_dist## Response: stfgov (numeric)
## Explanatory: clsprty (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 0.366
## 2 2 -1.62
## 3 3 0.632
## 4 4 0.243
## 5 5 -1.18
## 6 6 0.481
## 7 7 0.387
## 8 8 -0.0225
## 9 9 -2.00
## 10 10 -0.210
## # ℹ 990 more rowsconf_int <- null_dist %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_dist %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int)+
labs(title = "Figure 2. Simulation-Based Null Distribution clsprty~stfgov")
test_stat <- belgium_data %>%
specify(explanatory = clsprty,
response = stfgov) %>%
hypothesize(null = "independence") %>%
calculate(stat = "t", order = c("No", "Yes"))
print(test_stat$stat)## t
## -8.364563Univariate regression model
Using political interest (polintr) as a single predictor variable to predict government satisfaction (stfgov).
model_single1 <- lm(stfgov ~ polintr, data = belgium_data)
summary(model_single1)##
## Call:
## lm(formula = stfgov ~ polintr, data = belgium_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8725 -1.6786 0.3214 1.3214 5.9408
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.67858 0.02789 167.757 < 2e-16 ***
## polintrNot at all interested -0.61935 0.04625 -13.391 < 2e-16 ***
## polintrQuite interested 0.19390 0.03865 5.017 5.31e-07 ***
## polintrVery interested 0.13411 0.05949 2.254 0.0242 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.117 on 16931 degrees of freedom
## Multiple R-squared: 0.01945, Adjusted R-squared: 0.01928
## F-statistic: 112 on 3 and 16931 DF, p-value: < 2.2e-16model_single2 <- lm(stfgov ~ clsprty, data = belgium_data)
summary(model_single2)##
## Call:
## lm(formula = stfgov ~ clsprty, data = belgium_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7773 -1.5029 0.2227 1.4971 5.4971
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.50295 0.02340 192.400 <2e-16 ***
## clsprtyYes 0.27435 0.03279 8.366 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.133 on 16933 degrees of freedom
## Multiple R-squared: 0.004117, Adjusted R-squared: 0.004058
## F-statistic: 69.99 on 1 and 16933 DF, p-value: < 2.2e-16coef(model_single2)## (Intercept) clsprtyYes
## 4.5029486 0.2743526# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "interactions") # 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)## Warning: 程辑包'interactions'是用R版本4.3.3 来建造的##
## 载入程辑包：'interactions'## The following object is masked from 'package:ggeffects':
##
## johnson_neyman## [[1]]
## [1] "marginaleffects" "ggeffects" "scales" "broom"
## [5] "flextable" "performance" "survey" "survival"
## [9] "Matrix" "grid" "effects" "carData"
## [13] "infer" "fstcore" "modelsummary" "fst"
## [17] "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble"
## [25] "ggplot2" "tidyverse" "stats" "graphics"
## [29] "grDevices" "utils" "datasets" "methods"
## [33] "base"
##
## [[2]]
## [1] "marginaleffects" "ggeffects" "scales" "broom"
## [5] "flextable" "performance" "survey" "survival"
## [9] "Matrix" "grid" "effects" "carData"
## [13] "infer" "fstcore" "modelsummary" "fst"
## [17] "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble"
## [25] "ggplot2" "tidyverse" "stats" "graphics"
## [29] "grDevices" "utils" "datasets" "methods"
## [33] "base"
##
## [[3]]
## [1] "marginaleffects" "ggeffects" "scales" "broom"
## [5] "flextable" "performance" "survey" "survival"
## [9] "Matrix" "grid" "effects" "carData"
## [13] "infer" "fstcore" "modelsummary" "fst"
## [17] "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble"
## [25] "ggplot2" "tidyverse" "stats" "graphics"
## [29] "grDevices" "utils" "datasets" "methods"
## [33] "base"
##
## [[4]]
## [1] "marginaleffects" "ggeffects" "scales" "broom"
## [5] "flextable" "performance" "survey" "survival"
## [9] "Matrix" "grid" "effects" "carData"
## [13] "infer" "fstcore" "modelsummary" "fst"
## [17] "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble"
## [25] "ggplot2" "tidyverse" "stats" "graphics"
## [29] "grDevices" "utils" "datasets" "methods"
## [33] "base"
##
## [[5]]
## [1] "marginaleffects" "ggeffects" "scales" "broom"
## [5] "flextable" "performance" "survey" "survival"
## [9] "Matrix" "grid" "effects" "carData"
## [13] "infer" "fstcore" "modelsummary" "fst"
## [17] "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble"
## [25] "ggplot2" "tidyverse" "stats" "graphics"
## [29] "grDevices" "utils" "datasets" "methods"
## [33] "base"
##
## [[6]]
## [1] "marginaleffects" "ggeffects" "scales" "broom"
## [5] "flextable" "performance" "survey" "survival"
## [9] "Matrix" "grid" "effects" "carData"
## [13] "infer" "fstcore" "modelsummary" "fst"
## [17] "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble"
## [25] "ggplot2" "tidyverse" "stats" "graphics"
## [29] "grDevices" "utils" "datasets" "methods"
## [33] "base"
##
## [[7]]
## [1] "interactions" "marginaleffects" "ggeffects" "scales"
## [5] "broom" "flextable" "performance" "survey"
## [9] "survival" "Matrix" "grid" "effects"
## [13] "carData" "infer" "fstcore" "modelsummary"
## [17] "fst" "lubridate" "forcats" "stringr"
## [21] "dplyr" "purrr" "readr" "tidyr"
## [25] "tibble" "ggplot2" "tidyverse" "stats"
## [29] "graphics" "grDevices" "utils" "datasets"
## [33] "methods" "base"models <- list(
"Model 1" = lm(stfgov ~ polintr, data = belgium_data),
"Model 2" = lm(stfgov ~ clsprty, data = belgium_data))
modelsummary(models, fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL)| Model 1 | Model 2 | |
|---|---|---|
| (Intercept) | 4.7 (0.0)*** | 4.5 (0.0)*** |
| polintrNot at all interested | -0.6 (0.0)*** | |
| polintrQuite interested | 0.2 (0.0)*** | |
| polintrVery interested | 0.1 (0.1)* | |
| clsprtyYes | 0.3 (0.0)*** | |
| Num.Obs. | 16935 | 16935 |
| R2 | 0.019 | 0.004 |
| R2 Adj. | 0.019 | 0.004 |
| AIC | 73467.9 | 73726.8 |
| BIC | 73506.6 | 73750.0 |
| Log.Lik. | -36728.975 | -36860.396 |
| RMSE | 2.12 | 2.13 |
modelplot(models, coef_omit = 'Interc')
Bivariate regression mode
Using political interest (polintr) and party affiliation (clsprty) as predictor variables to simultaneously predict government satisfaction (stfgov).
model_double <- lm(stfgov ~ polintr + clsprty, data = belgium_data)
summary(model_double)##
## Call:
## lm(formula = stfgov ~ polintr + clsprty, data = belgium_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9224 -1.6209 0.2474 1.3791 5.9781
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.62086 0.03165 146.006 < 2e-16 ***
## polintrNot at all interested -0.59900 0.04653 -12.873 < 2e-16 ***
## polintrQuite interested 0.16979 0.03914 4.338 1.45e-05 ***
## polintrVery interested 0.08803 0.06066 1.451 0.146719
## clsprtyYes 0.13171 0.03419 3.853 0.000117 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.116 on 16930 degrees of freedom
## Multiple R-squared: 0.02031, Adjusted R-squared: 0.02008
## F-statistic: 87.76 on 4 and 16930 DF, p-value: < 2.2e-16Interaction regression model.
Building upon the bivariate model, adding the interaction term between political interest and party affiliation to examine the interaction effect on government satisfaction.
model_interaction <- lm(stfgov ~ polintr * clsprty, data = belgium_data)
summary(model_interaction)##
## Call:
## lm(formula = stfgov ~ polintr * clsprty, data = belgium_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9576 -1.6577 0.2947 1.3423 5.9521
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.657708 0.037189 125.246 <2e-16
## polintrNot at all interested -0.609786 0.057284 -10.645 <2e-16
## polintrQuite interested 0.075204 0.057201 1.315 0.1886
## polintrVery interested -0.076312 0.119986 -0.636 0.5248
## clsprtyYes 0.047639 0.056178 0.848 0.3965
## polintrNot at all interested:clsprtyYes -0.007766 0.099235 -0.078 0.9376
## polintrQuite interested:clsprtyYes 0.177012 0.078716 2.249 0.0245
## polintrVery interested:clsprtyYes 0.245868 0.140249 1.753 0.0796
##
## (Intercept) ***
## polintrNot at all interested ***
## polintrQuite interested
## polintrVery interested
## clsprtyYes
## polintrNot at all interested:clsprtyYes
## polintrQuite interested:clsprtyYes *
## polintrVery interested:clsprtyYes .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.116 on 16927 degrees of freedom
## Multiple R-squared: 0.02077, Adjusted R-squared: 0.02037
## F-statistic: 51.29 on 7 and 16927 DF, p-value: < 2.2e-16mode666 <- lm(stfgov ~ +polintr+clsprty+polintr * clsprty, data = belgium_data)
summary(mode666)##
## Call:
## lm(formula = stfgov ~ +polintr + clsprty + polintr * clsprty,
## data = belgium_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9576 -1.6577 0.2947 1.3423 5.9521
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.657708 0.037189 125.246 <2e-16
## polintrNot at all interested -0.609786 0.057284 -10.645 <2e-16
## polintrQuite interested 0.075204 0.057201 1.315 0.1886
## polintrVery interested -0.076312 0.119986 -0.636 0.5248
## clsprtyYes 0.047639 0.056178 0.848 0.3965
## polintrNot at all interested:clsprtyYes -0.007766 0.099235 -0.078 0.9376
## polintrQuite interested:clsprtyYes 0.177012 0.078716 2.249 0.0245
## polintrVery interested:clsprtyYes 0.245868 0.140249 1.753 0.0796
##
## (Intercept) ***
## polintrNot at all interested ***
## polintrQuite interested
## polintrVery interested
## clsprtyYes
## polintrNot at all interested:clsprtyYes
## polintrQuite interested:clsprtyYes *
## polintrVery interested:clsprtyYes .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.116 on 16927 degrees of freedom
## Multiple R-squared: 0.02077, Adjusted R-squared: 0.02037
## F-statistic: 51.29 on 7 and 16927 DF, p-value: < 2.2e-16modelsummary(
list(model_double, mode666),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL)| (1) | (2) | |
|---|---|---|
| (Intercept) | 4.6 (0.0)*** | 4.7 (0.0)*** |
| polintrNot at all interested | -0.6 (0.0)*** | -0.6 (0.1)*** |
| polintrQuite interested | 0.2 (0.0)*** | 0.1 (0.1) |
| polintrVery interested | 0.1 (0.1) | -0.1 (0.1) |
| clsprtyYes | 0.1 (0.0)*** | 0.0 (0.1) |
| polintrNot at all interested × clsprtyYes | 0.0 (0.1) | |
| polintrQuite interested × clsprtyYes | 0.2 (0.1)* | |
| polintrVery interested × clsprtyYes | 0.2 (0.1)+ | |
| Num.Obs. | 16935 | 16935 |
| R2 | 0.020 | 0.021 |
| R2 Adj. | 0.020 | 0.020 |
| AIC | 73455.1 | 73453.2 |
| BIC | 73501.5 | 73522.8 |
| Log.Lik. | -36721.555 | -36717.591 |
| RMSE | 2.12 | 2.12 |
interaction_plot <- effect("polintr*clsprty", mode666)
plot(interaction_plot,
main = "Figure 3. Interaction Effect on Government Satisfaction",
xlab = "Political Interest",
ylab = "Government Satisfaction")
mode666 <- lm(stfgov ~ +polintr+clsprty+polintr * clsprty, data = belgium_data)
interaction.plot(belgium_data$clsprty, belgium_data$polintr, belgium_data$stfgov)
cat_plot(mode666, pred = polintr, modx = clsprty, jnplot = TRUE)
check_model(model_single1)
check_model(model_single2)
check_model(model_double)
check_model(mode666)