Research Report
2024-04-09
# List of packages
packages <- c("tidyverse", "infer", "fst", "modelsummary", "effects", "survey", "MASS", "aod", "interactions", "kableExtra", "flextable", "scales", "parsnip") # 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: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select## Warning: package 'aod' was built under R version 4.3.3##
## Attaching package: 'aod'
##
## The following object is masked from 'package:survival':
##
## rats## Warning: package 'interactions' was built under R version 4.3.3##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows## Warning: package 'flextable' was built under R version 4.3.3##
## 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## Warning: package 'parsnip' 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] "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] "flextable" "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] "scales" "flextable" "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] "parsnip" "scales" "flextable" "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"Context
my_data <- read_fst("C:\\Users\\eesha\\Documents\\Report\\sweden_data.fst")# clean variables
sweden_data <- my_data %>%
mutate(
trstlgl = ifelse(trstlgl %in% c(77, 88, 99), NA, trstlgl),
)
table(sweden_data$trstlgl)##
## 0 1 2 3 4 5 6 7 8 9 10
## 345 311 669 1148 1385 2540 2275 3449 3582 1521 612sweden_data$year <- NA
replacements <- c(2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020)
for(i in 1:10){
sweden_data$year[sweden_data$essround == i] <- replacements[i]
}
#filter for 2020
sweden_data_2020 <- sweden_data %>%
filter(year == 2020)
# Check the table of the 'year' column in the filtered data
table(sweden_data_2020$year)##
## 2020
## 2287library(dplyr)
# Filter data for the year 2020
sweden_data_2020 <- sweden_data %>%
filter(year == 2020)
# Calculate the mean trust for 2020
trust_by_year_2020 <- sweden_data_2020 %>%
summarize(mean_trust = mean(trstlgl, na.rm = TRUE))
# Print the mean trust for 2020
print(trust_by_year_2020)## mean_trust
## 1 5.824229sweden_data_2020 <- sweden_data %>%
filter(year == 2020)
# Calculate the standard deviation (sd) for trust for 2020
sd_trust_2020 <- sweden_data_2020 %>%
summarize(sd_trust = sd(trstlgl, na.rm = TRUE))
# Print the standard deviation for trust in 2020
print(sd_trust_2020)## sd_trust
## 1 2.577906sweden_data_2020 <- sweden_data %>%
filter(year == 2020)
#group lrscale
sweden_data_2020 <- sweden_data_2020 %>%
mutate(
lrscale = case_when(
lrscale %in% 0:1 ~ "1", # far left
lrscale %in% 2:5 ~ "2", # center left
lrscale %in% c(6,7,8,9,10,77, 88, 99) ~ NA_character_ # NA for other values
)
)
# Calculate counts for each level of 'lrscale' in 2020 data
count_1 <- sum(sweden_data_2020$lrscale == "1", na.rm = TRUE)
count_2 <- sum(sweden_data_2020$lrscale == "2", na.rm = TRUE)
# Create a data frame to display the counts
counts_df <- data.frame(Level = c("1", "2"), Count = c(count_1, count_2))
# Print the counts data frame
print(counts_df)## Level Count
## 1 1 203
## 2 2 1077# Mutate the 'educ.ba' variable
sweden_data_2020 <- sweden_data_2020 %>%
mutate(
# Recoding education - only edulvlb recoded as edulvla not used in 2020
educ.ba = case_when(
essround >= 5 & edulvlb > 600 ~ "2", # BA or more
TRUE ~ "1" # no BA
),
# Handle NAs for education levels
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
# Convert 'educ.ba' to a factor with specified levels
educ.ba = factor(educ.ba, levels = c("1", "2"))
)
# Create a table of counts for 'educ.ba' with labeled columns
educ_table <- table(sweden_data_2020$educ.ba)
dimnames(educ_table) <- list(Education_Level = c("1", "2"))
print(educ_table)## Education_Level
## 1 2
## 1549 738For the next step (datasummary) this was the code used but there were some errors
summary_legal <- datasummary(legal_data) summary_trust <- datasummary(trust_by_year) summary_gender <- datasummary(gender_data) summary_edu <- datasummary(edu_data)
summary_table <- datasummary(summary_legal, summary_trust, summary_gender, summary_edu) summary_table
Instead I renamed the data in hopes of it working.
This was the next code used; it did not work as the argument “data” was missing.
summary_table <- datasummary( “Legal Data” = legal_data, “Trust by Year” = trust_by_year, “Gender Data” = gender_data, “Education Data” = edu_data )
summary_table
This was the next code I tried:
library(modelsummary) library(dplyr)
cleaned_data <- sweden_data_2020 %>% filter(!is.na(lrscale) & !is.na(educ.ba) & !lrscale %in% c(77, 88, 99)) %>% mutate( lrscale = as.numeric(as.character(lrscale)), educ.ba = as.numeric(as.character(educ.ba)) )
summary_skim <- cleaned_data %>% datasummary_skim(statistics = list(mean, sd))
print(summary_skim)
This also did not work. I kept trying multiple different codes and eventually landed on this one:
library(modelsummary)
library(dplyr)
# change variables to numeric to fit datasummary model
subset_data <- sweden_data_2020 %>%
dplyr::select("trstlgl", "lrscale", "educ.ba") %>%
mutate(
lrscale = as.numeric(as.character(lrscale)),
educ.ba = as.numeric(as.character(educ.ba))
)
# filter NA values
subset_data <- subset_data %>%
filter(!is.na(lrscale) & !lrscale %in% c(77, 88, 99) & !is.na(trstlgl) &!trstlgl %in% c(77, 88, 99) & !is.na(educ.ba) & !educ.ba %in% c(77, 88, 99)) %>%
mutate(
lrscale = as.numeric(as.character(lrscale)),
educ.ba = as.numeric(as.character(educ.ba)),
)
summary_skim <- subset_data %>%
datasummary_skim(
statistics = list(mean, sd, nobs) ,
title = "Table 1"
)
print(summary_skim)## <table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
## <caption>Table 1</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> </th>
## <th style="text-align:right;"> Unique (#) </th>
## <th style="text-align:right;"> Missing (%) </th>
## <th style="text-align:right;"> Mean </th>
## <th style="text-align:right;"> SD </th>
## <th style="text-align:right;"> Min </th>
## <th style="text-align:right;"> Median </th>
## <th style="text-align:right;"> Max </th>
## <th style="text-align:right;"> </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> trstlgl </td>
## <td style="text-align:right;"> 11 </td>
## <td style="text-align:right;"> 0 </td>
## <td style="text-align:right;"> 6.0 </td>
## <td style="text-align:right;"> 2.5 </td>
## <td style="text-align:right;"> 0.0 </td>
## <td style="text-align:right;"> 7.0 </td>
## <td style="text-align:right;"> 10.0 </td>
## <td style="text-align:right;"> <svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" class="svglite" width="48.00pt" height="12.00pt" viewBox="0 0 48.00 12.00"><defs><style type="text/css">
## .svglite line, .svglite polyline, .svglite polygon, .svglite path, .svglite rect, .svglite circle {
## fill: none;
## stroke: #000000;
## stroke-linecap: round;
## stroke-linejoin: round;
## stroke-miterlimit: 10.00;
## }
## .svglite text {
## white-space: pre;
## }
## </style></defs><rect width="100%" height="100%" style="stroke: none; fill: none;"></rect><defs><clipPath id="cpMC4wMHw0OC4wMHwwLjAwfDEyLjAw"><rect x="0.00" y="0.00" width="48.00" height="12.00"></rect></clipPath></defs><g clip-path="url(#cpMC4wMHw0OC4wMHwwLjAwfDEyLjAw)">
## </g><defs><clipPath id="cpMC4wMHw0OC4wMHwyLjg4fDEyLjAw"><rect x="0.00" y="2.88" width="48.00" height="9.12"></rect></clipPath></defs><g clip-path="url(#cpMC4wMHw0OC4wMHwyLjg4fDEyLjAw)"><rect x="1.78" y="9.31" width="4.44" height="2.35" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="6.22" y="9.18" width="4.44" height="2.49" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="10.67" y="8.36" width="4.44" height="3.30" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="15.11" y="9.38" width="4.44" height="2.28" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="19.56" y="5.43" width="4.44" height="6.23" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="24.00" y="6.79" width="4.44" height="4.87" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="28.44" y="4.68" width="4.44" height="6.98" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="32.89" y="3.22" width="4.44" height="8.44" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="37.33" y="7.64" width="4.44" height="4.02" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="41.78" y="9.35" width="4.44" height="2.32" style="stroke-width: 0.38; fill: #000000;"></rect></g></svg>
## </td>
## </tr>
## <tr>
## <td style="text-align:left;"> lrscale </td>
## <td style="text-align:right;"> 2 </td>
## <td style="text-align:right;"> 0 </td>
## <td style="text-align:right;"> 1.8 </td>
## <td style="text-align:right;"> 0.4 </td>
## <td style="text-align:right;"> 1.0 </td>
## <td style="text-align:right;"> 2.0 </td>
## <td style="text-align:right;"> 2.0 </td>
## <td style="text-align:right;"> <svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" class="svglite" width="48.00pt" height="12.00pt" viewBox="0 0 48.00 12.00"><defs><style type="text/css">
## .svglite line, .svglite polyline, .svglite polygon, .svglite path, .svglite rect, .svglite circle {
## fill: none;
## stroke: #000000;
## stroke-linecap: round;
## stroke-linejoin: round;
## stroke-miterlimit: 10.00;
## }
## .svglite text {
## white-space: pre;
## }
## </style></defs><rect width="100%" height="100%" style="stroke: none; fill: none;"></rect><defs><clipPath id="cpMC4wMHw0OC4wMHwwLjAwfDEyLjAw"><rect x="0.00" y="0.00" width="48.00" height="12.00"></rect></clipPath></defs><g clip-path="url(#cpMC4wMHw0OC4wMHwwLjAwfDEyLjAw)">
## </g><defs><clipPath id="cpMC4wMHw0OC4wMHwyLjg4fDEyLjAw"><rect x="0.00" y="2.88" width="48.00" height="9.12"></rect></clipPath></defs><g clip-path="url(#cpMC4wMHw0OC4wMHwyLjg4fDEyLjAw)"><rect x="1.78" y="10.08" width="4.44" height="1.59" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="6.22" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="10.67" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="15.11" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="19.56" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="24.00" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="28.44" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="32.89" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="37.33" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="41.78" y="3.22" width="4.44" height="8.44" style="stroke-width: 0.38; fill: #000000;"></rect></g></svg>
## </td>
## </tr>
## <tr>
## <td style="text-align:left;"> educ.ba </td>
## <td style="text-align:right;"> 2 </td>
## <td style="text-align:right;"> 0 </td>
## <td style="text-align:right;"> 1.3 </td>
## <td style="text-align:right;"> 0.5 </td>
## <td style="text-align:right;"> 1.0 </td>
## <td style="text-align:right;"> 1.0 </td>
## <td style="text-align:right;"> 2.0 </td>
## <td style="text-align:right;"> <svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" class="svglite" width="48.00pt" height="12.00pt" viewBox="0 0 48.00 12.00"><defs><style type="text/css">
## .svglite line, .svglite polyline, .svglite polygon, .svglite path, .svglite rect, .svglite circle {
## fill: none;
## stroke: #000000;
## stroke-linecap: round;
## stroke-linejoin: round;
## stroke-miterlimit: 10.00;
## }
## .svglite text {
## white-space: pre;
## }
## </style></defs><rect width="100%" height="100%" style="stroke: none; fill: none;"></rect><defs><clipPath id="cpMC4wMHw0OC4wMHwwLjAwfDEyLjAw"><rect x="0.00" y="0.00" width="48.00" height="12.00"></rect></clipPath></defs><g clip-path="url(#cpMC4wMHw0OC4wMHwwLjAwfDEyLjAw)">
## </g><defs><clipPath id="cpMC4wMHw0OC4wMHwyLjg4fDEyLjAw"><rect x="0.00" y="2.88" width="48.00" height="9.12"></rect></clipPath></defs><g clip-path="url(#cpMC4wMHw0OC4wMHwyLjg4fDEyLjAw)"><rect x="1.78" y="3.22" width="4.44" height="8.44" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="6.22" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="10.67" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="15.11" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="19.56" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="24.00" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="28.44" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="32.89" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="37.33" y="11.66" width="4.44" height="0.00" style="stroke-width: 0.38; fill: #000000;"></rect><rect x="41.78" y="7.49" width="4.44" height="4.17" style="stroke-width: 0.38; fill: #000000;"></rect></g></svg>
## </td>
## </tr>
## </tbody>
## </table>write(summary_skim, "skimtable.pdf") # export as pdfFindings
#print data
table_trstlgl <- table(subset_data$trstlgl)
print(table_trstlgl)##
## 0 1 2 3 4 5 6 7 8 9 10
## 42 27 73 97 67 183 143 205 248 118 68table_lrscale <- table(subset_data$lrscale)
print(table_lrscale)##
## 1 2
## 201 1070table_educ.ba <- table(subset_data$educ.ba)
print(table_educ.ba)##
## 1 2
## 851 420This was the first code I tried for the next step.
library(infer)
library(dplyr)
subset_data <- subset_data %>%
filter(!is.na(trstlgl), !is.na(lrscale), !is.na(educ.ba)) %>%
mutate(
trstlgl = factor(trstlgl),
lrscale = factor(lrscale),
educ.ba = factor(educ.ba)
)
test_stat <- subset_data %>%
specify(explanatory = lrscale + educ.ba,
response = trstlgl) %>%
hypothesize(null = "independence") %>%
fit(stat = "chisq") # fit instead of calculate because there are 2 explanatory variables
print(test_stat$stat)Error read: The independence null hypothesis does not inform
calculation of the observed fit and will be ignored.Error in
fit(): ! infer does not support fitting models for
categorical response variables with more than two levels.
Error said to use multinom_reg() instead; here is the subsequent code:
library(infer)
library(dplyr)
library(parsnip)
# Assuming subset_data is your dataframe
# Ensure columns are factors and remove NA values
subset_data <- subset_data %>%
mutate(
trstlgl = factor(trstlgl),
lrscale = factor(lrscale),
educ.ba = factor(educ.ba)
) %>%
filter(!is.na(trstlgl), !is.na(lrscale), !is.na(educ.ba))
# Fit multinomial logistic regression model
model <- multinom_reg() %>%
set_engine("nnet", trace = FALSE) %>%
set_mode("classification") %>%
fit(trstlgl ~ lrscale + educ.ba, data = subset_data)
# Perform hypothesis test
test_stat <- subset_data %>%
specify(
formula = trstlgl ~ lrscale + educ.ba,
success = "5"
) %>%
hypothesize(null = "independence") %>%
fit(model)
print(test_stat$stat)
This did not work as I had not changed the response to 2 levels. So i added that with 5 and up being confident and under 5 being not confident and ran a different way
These methods did not work and I realized the null has to be run separately not interaction for both variables so I did that.
table_educ.ba <- table(subset_data$educ.ba)
print(table_educ.ba)##
## 1 2
## 851 420table_trstlgl <- table(subset_data$trstlgl)
print(table_trstlgl)##
## 0 1 2 3 4 5 6 7 8 9 10
## 42 27 73 97 67 183 143 205 248 118 68#rename variables back to original groups
subset_clean <- subset_data %>%
mutate(
ID = case_when(
educ.ba %in% 1 ~ "BA or More",
educ.ba %in% 2 ~ "No BA",
educ.ba %in% c(77, 88, 99) ~ NA_character_
)
)
table(subset_clean$ID)##
## BA or More No BA
## 851 420library(infer)
library(dplyr)
# filter for NA
subset_data <- subset_clean %>%
filter(!is.na(educ.ba), !is.na(trstlgl))
# factor
subset_data <- subset_data %>%
mutate(
educ.ba = factor(educ.ba),
trstlgl = factor(trstlgl)
)
# get test statistic
test_stat <- subset_data %>%
specify(explanatory = educ.ba,
response = trstlgl) %>%
hypothesize(null = "independence") %>%
calculate(stat = "chisq")
print(test_stat$stat)## X-squared
## 48.93853Chi-squared statistic is 48.93853 This statistic is the deviation of data to be expected if the null hypothesis were true. We will use this number to test the null.
# set variables and test stat
null_distribution <- subset_data %>%
specify(explanatory = educ.ba,
response = trstlgl) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "chisq")
# get p value
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.p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0p value is read as 0. Because this was gained from an approximation the p value is not truly 0.000 but this gives enough data to know that it is under the alpha threshold of 0.05. The real p value will be seen in our regression model.
# null distribution
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided")## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
null_distribution## Response: trstlgl (factor)
## Explanatory: educ.ba (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 7.36
## 2 2 8.32
## 3 3 15.4
## 4 4 6.37
## 5 5 8.23
## 6 6 6.41
## 7 7 3.87
## 8 8 5.49
## 9 9 5.42
## 10 10 11.2
## # ℹ 990 more rows# confidence intervals
conf_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)## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
null_distribution## Response: trstlgl (factor)
## Explanatory: educ.ba (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 7.36
## 2 2 8.32
## 3 3 15.4
## 4 4 6.37
## 5 5 8.23
## 6 6 6.41
## 7 7 3.87
## 8 8 5.49
## 9 9 5.42
## 10 10 11.2
## # ℹ 990 more rowsThis visual shows that the test statistic is far from the null distribution; therefore we reject the null.
Now let us test using the other variable: lrscale
table_trstlgl <- table(subset_data$trstlgl)
print(table_trstlgl)##
## 0 1 2 3 4 5 6 7 8 9 10
## 42 27 73 97 67 183 143 205 248 118 68table_lrscale <- table(subset_data$lrscale)
print(table_lrscale)##
## 1 2
## 201 1070subset_clean <- subset_data %>%
mutate(
ID = case_when(
lrscale %in% 0:1 ~ "Far Left",
lrscale %in% 2:5 ~ "Center Left",
lrscale %in% c(77, 88, 99) ~ NA_character_
)
)
table(subset_clean$ID)##
## Center Left Far Left
## 1070 201library(infer)
library(dplyr)
subset_data <- subset_clean %>%
filter(!is.na(lrscale), !is.na(trstlgl))
subset_data <- subset_data %>%
mutate(
lrscale = factor(lrscale),
trstlgl = factor(trstlgl)
)
test_stat <- subset_data %>%
specify(explanatory = lrscale,
response = trstlgl) %>%
hypothesize(null = "independence") %>%
calculate(stat = "chisq")
print(test_stat$stat)## X-squared
## 19.53151Chi-square test is 19.53151. This will be our test statistic
null_distribution <- subset_data %>%
specify(explanatory = lrscale,
response = trstlgl) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "chisq")
p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "two-sided")
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.066p value is 0.062; over the alpha threshold 0.05
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided")## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
null_distribution## Response: trstlgl (factor)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 9.40
## 2 2 11.3
## 3 3 9.48
## 4 4 7.42
## 5 5 12.5
## 6 6 7.73
## 7 7 4.73
## 8 8 14.9
## 9 9 7.75
## 10 10 8.23
## # ℹ 990 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)## Warning: Chi-Square usually corresponds to right-tailed tests. Proceed with
## caution.
null_distribution## Response: trstlgl (factor)
## Explanatory: lrscale (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 9.40
## 2 2 11.3
## 3 3 9.48
## 4 4 7.42
## 5 5 12.5
## 6 6 7.73
## 7 7 4.73
## 8 8 14.9
## 9 9 7.75
## 10 10 8.23
## # ℹ 990 more rowsThis time the test statistic lies within the null distribution; therefore we cannot reject the null.
Linear Regression Models
table(subset_data$trstlgl)##
## 0 1 2 3 4 5 6 7 8 9 10
## 42 27 73 97 67 183 143 205 248 118 68table(subset_data$educ.ba)##
## 1 2
## 851 420table(subset_data$lrscale)##
## 1 2
## 201 1070# Convert 'trstlgl' to numeric
subset_data$trstlgl <- as.numeric(as.character(subset_data$trstlgl))
# regression model for main predictor
model1 <- lm(trstlgl ~ educ.ba, data = subset_data)
#multivariate model
model2 <- lm(trstlgl ~ educ.ba + lrscale, data = subset_data)
#interaction model
model3 <- lm(trstlgl ~ educ.ba + lrscale + educ.ba*lrscale, data = subset_data)
# Summarize the models
summary(model1)##
## Call:
## lm(formula = trstlgl ~ educ.ba, data = subset_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.6667 -1.6667 0.3333 2.2820 4.2820
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.71798 0.08483 67.403 < 2e-16 ***
## educ.ba2 0.94869 0.14757 6.429 1.82e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.475 on 1269 degrees of freedom
## Multiple R-squared: 0.03154, Adjusted R-squared: 0.03078
## F-statistic: 41.33 on 1 and 1269 DF, p-value: 1.819e-10summary(model2)##
## Call:
## lm(formula = trstlgl ~ educ.ba + lrscale, data = subset_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.8679 -1.6325 0.3208 2.0853 4.3208
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.9147 0.1802 32.828 < 2e-16 ***
## educ.ba2 0.9532 0.1476 6.459 1.5e-10 ***
## lrscale2 -0.2354 0.1903 -1.237 0.216
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.474 on 1268 degrees of freedom
## Multiple R-squared: 0.03271, Adjusted R-squared: 0.03118
## F-statistic: 21.44 on 2 and 1268 DF, p-value: 6.978e-10summary(model3)##
## Call:
## lm(formula = trstlgl ~ educ.ba + lrscale + educ.ba * lrscale,
## data = subset_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1148 -1.5905 0.2996 2.1929 4.2996
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.8071 0.2091 27.771 < 2e-16 ***
## educ.ba2 1.3076 0.3796 3.445 0.00059 ***
## lrscale2 -0.1067 0.2288 -0.467 0.64094
## educ.ba2:lrscale2 -0.4175 0.4120 -1.013 0.31108
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.474 on 1267 degrees of freedom
## Multiple R-squared: 0.03349, Adjusted R-squared: 0.0312
## F-statistic: 14.63 on 3 and 1267 DF, p-value: 2.258e-09# models list
models_list <- list(model1, model2, model3)
# generate summary and title
modelsummary(
models_list,
fmt = 1,
estimate = c("{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}",
"{estimate} ({std.error}){stars}"),
statistic = NULL
)| (1) | (2) | (3) | |
|---|---|---|---|
| (Intercept) | 5.7 (0.1)*** | 5.9 (0.2)*** | 5.8 (0.2)*** |
| educ.ba2 | 0.9 (0.1)*** | 1.0 (0.1)*** | 1.3 (0.4)*** |
| lrscale2 | −0.2 (0.2) | −0.1 (0.2) | |
| educ.ba2 × lrscale2 | −0.4 (0.4) | ||
| Num.Obs. | 1271 | 1271 | 1271 |
| R2 | 0.032 | 0.033 | 0.033 |
| R2 Adj. | 0.031 | 0.031 | 0.031 |
| AIC | 5914.3 | 5914.8 | 5915.8 |
| BIC | 5929.8 | 5935.4 | 5941.5 |
| Log.Lik. | −2954.157 | −2953.390 | −2952.875 |
| RMSE | 2.47 | 2.47 | 2.47 |
Interaction effect
interaction_plot <- effect("educ.ba * lrscale", model3, na.rm = TRUE)
plot(interaction_plot,
main = "Interaction effect",
xlab = "Education",
ylab = "Trust in legal system")
Lines are almost parallel indicating low statistical interaction.
```