Final Report
Melanie Fidel Rodriguez
2024-03-25
# 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## Warning: package 'flextable' was built under R version 4.3.3##
## 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 '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"france_data <- read_fst("france_data.fst")table(france_data$polID)## < table of extent 0 >table(france_data$educ.ba)## < table of extent 0 >table(france_data$impenv)##
## 1 2 3 4 5 6 7 8 9
## 6507 5817 3394 2243 572 148 31 59 125france_data <- france_data %>%
# Recoding political ID and education
mutate(
polID = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
lrscale %in% 4:6 ~ "Moderate",
lrscale %in% c(77, 88, 99) ~ NA_character_
),
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
),
edulvla = ifelse(edulvla %in% c(77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = factor(educ.ba, levels = c("No BA", "BA or more"))
)Descriptive table
You will need your outcome and your two main predictors.
france_data <- france_data %>%
mutate(
impenv = impenv,
polID = case_when(
lrscale > 10 ~ NA_real_, # Set values above 10 to NA
TRUE ~ lrscale # Keep other values as is
),
educ.ba = ifelse(educ.ba %in% c(7, 8, 9), NA, educ.ba),
)
df <- france_data
df <- df %>% filter(!is.na(polID))
df <- df %>% filter(!is.na(educ.ba))
table1b <- datasummary_skim(df %>% dplyr::select(impenv, polID, educ.ba), 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.table1b
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
impenv | 10 | 1 | 2.2 | 1.3 | 1.0 | 2.0 | 9.0 |
polID | 11 | 0 | 4.9 | 2.4 | 0.0 | 5.0 | 10.0 |
educ.ba | 2 | 0 | 1.2 | 0.4 | 1.0 | 1.0 | 2.0 |
Null Distribution
france_task1 <- france_data %>%
mutate(
impenv = case_when (
impenv == 1 ~ "Very Much Like Me",
impenv == 2 ~ "Like Me",
impenv == 3 ~ "Somewhat Like Me",
impenv == 4 ~ "A Little Like Me",
impenv == 5 ~ "Not Like Me",
impenv == 6 ~ "Not Like Me At All",
impenv %in% c(7,8,9) ~ NA_character_,
TRUE ~ as.character(impenv)
),
polID = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
lrscale %in% 4:6 ~ "Moderate",
lrscale %in% c(77, 88, 99) ~ NA_character_)
) %>%
filter (!is.na(impenv), !is.na(polID)) mytab <- table(france_task1$impenv, france_task1$polID)
mytab##
## Left Moderate Right
## A Little Like Me 460 1091 530
## Like Me 1457 2735 1259
## Not Like Me 113 254 152
## Not Like Me At All 33 60 35
## Somewhat Like Me 793 1564 824
## Very Much Like Me 1788 2855 1432#step 1:
test_stat <- france_task1 %>%
specify(explanatory = polID,
response = impenv) %>%
hypothesize(null = "independence") %>%
calculate(stat = "chisq")
print(test_stat$stat)## X-squared
## 68.51947#step 2:
null_distribution <- france_task1 %>%
specify(response = impenv, explanatory = polID) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "Chisq")
null_distribution## Response: impenv (factor)
## Explanatory: polID (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 14.3
## 2 2 4.84
## 3 3 8.56
## 4 4 10.6
## 5 5 12.9
## 6 6 10.2
## 7 7 6.77
## 8 8 6.19
## 9 9 19.9
## 10 10 7.64
## # ℹ 990 more rows#step 3:
p_val <- null_distribution %>% #
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 0#step 4:
conf_int <- null_distribution %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution %>%
visualize(data, bins = 10, method = "simulation", dens_color = "black") +
shade_p_value(obs_stat = test_stat, direction = "greater") +
shade_confidence_interval(endpoints = conf_int)## Warning in min(diff(unique_loc)): no non-missing arguments to min; returning
## Inf
#regression models
table(france_data$impenv)##
## 1 2 3 4 5 6 7 8 9
## 6507 5817 3394 2243 572 148 31 59 125model_data <- france_data %>%
mutate(
impenv = case_when (
impenv == 1 ~ "Very Much Like Me",
impenv == 2 ~ "Like Me",
impenv == 3 ~ "Somewhat Like Me",
impenv == 4 ~ "A Little Like Me",
impenv == 5 ~ "Not Like Me",
impenv == 6 ~ "Not Like Me At All",
impenv %in% c(7, 8, 9) ~ NA_character_,
TRUE ~ as.character(impenv)
),
polID = case_when(
lrscale %in% 0:3 ~ "Left",
lrscale %in% 7:10 ~ "Right",
lrscale %in% 4:6 ~ "Moderate",
lrscale %in% c(77, 88, 99) ~ NA_character_)
) %>%
filter (!is.na(impenv), !is.na(polID)) table(model_data$impenv)##
## A Little Like Me Like Me Not Like Me Not Like Me At All
## 2081 5451 519 128
## Somewhat Like Me Very Much Like Me
## 3181 6075table(model_data$polID)##
## Left Moderate Right
## 4644 8559 4232model_1 <- lm(impenv ~ polID, data = france_data)
summary(model_1)##
## Call:
## lm(formula = impenv ~ polID, data = france_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3821 -1.2002 -0.2521 0.7479 6.8778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.122169 0.022792 93.109 < 2e-16 ***
## polID 0.025995 0.004195 6.196 5.92e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.309 on 17592 degrees of freedom
## (1444 observations deleted due to missingness)
## Multiple R-squared: 0.002178, Adjusted R-squared: 0.002121
## F-statistic: 38.39 on 1 and 17592 DF, p-value: 5.917e-10model_2 <- lm(impenv ~ polID + educ.ba, data = france_data)
summary(model_2)##
## Call:
## lm(formula = impenv ~ polID + educ.ba, data = france_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3918 -1.2063 -0.2658 0.7342 6.8945
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.199166 0.038014 57.851 < 2e-16 ***
## polID 0.025211 0.004206 5.994 2.09e-09 ***
## educ.ba -0.059439 0.023488 -2.531 0.0114 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.309 on 17591 degrees of freedom
## (1444 observations deleted due to missingness)
## Multiple R-squared: 0.002541, Adjusted R-squared: 0.002427
## F-statistic: 22.4 on 2 and 17591 DF, p-value: 1.917e-10model_3 <- lm(impenv ~ polID + educ.ba + polID * educ.ba, data = france_data)
summary(model_3)##
## Call:
## lm(formula = impenv ~ polID + educ.ba + polID * educ.ba, data = france_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5411 -1.2225 -0.2657 0.7343 7.0324
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.47919 0.07007 35.383 < 2e-16 ***
## polID -0.03406 0.01315 -2.590 0.00961 **
## educ.ba -0.28765 0.05341 -5.385 7.33e-08 ***
## polID:educ.ba 0.04889 0.01028 4.756 1.99e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.308 on 17590 degrees of freedom
## (1444 observations deleted due to missingness)
## Multiple R-squared: 0.003822, Adjusted R-squared: 0.003652
## F-statistic: 22.49 on 3 and 17590 DF, p-value: 1.57e-14models_list <- list(model_1, model_2, model_3)
# Generate the table summary with a 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) | 2.1 (0.0)*** | 2.2 (0.0)*** | 2.5 (0.1)*** |
| polID | 0.0 (0.0)*** | 0.0 (0.0)*** | 0.0 (0.0)** |
| educ.ba | −0.1 (0.0)* | −0.3 (0.1)*** | |
| polID × educ.ba | 0.0 (0.0)*** | ||
| Num.Obs. | 17594 | 17594 | 17594 |
| R2 | 0.002 | 0.003 | 0.004 |
| R2 Adj. | 0.002 | 0.002 | 0.004 |
| AIC | 59418.2 | 59413.8 | 59393.2 |
| BIC | 59441.5 | 59444.9 | 59432.1 |
| Log.Lik. | −29706.111 | −29702.909 | −29691.602 |
| RMSE | 1.31 | 1.31 | 1.31 |
#visual of choice
model4 <- glm(impenv ~ polID + educ.ba + polID*educ.ba, family = poisson, data = df, weights = round(weight)) # both predictors + interaction
mydf <- predict_response(model4, terms = c("polID", "educ.ba"))
print(mydf, n = Inf)## # Predicted counts of impenv
##
## educ.ba: 1
##
## polID | Predicted | 95% CI
## ------------------------------
## 0 | 2.20 | 2.18, 2.21
## 1 | 2.22 | 2.20, 2.23
## 2 | 2.24 | 2.23, 2.25
## 3 | 2.26 | 2.25, 2.27
## 4 | 2.28 | 2.28, 2.29
## 5 | 2.31 | 2.30, 2.32
## 6 | 2.33 | 2.32, 2.34
## 7 | 2.35 | 2.34, 2.36
## 8 | 2.38 | 2.36, 2.39
## 9 | 2.40 | 2.39, 2.42
## 10 | 2.43 | 2.41, 2.44
##
## educ.ba: 2
##
## polID | Predicted | 95% CI
## ------------------------------
## 0 | 2.12 | 2.09, 2.16
## 1 | 2.17 | 2.14, 2.20
## 2 | 2.21 | 2.19, 2.24
## 3 | 2.26 | 2.24, 2.28
## 4 | 2.31 | 2.29, 2.33
## 5 | 2.36 | 2.34, 2.38
## 6 | 2.41 | 2.39, 2.43
## 7 | 2.46 | 2.43, 2.48
## 8 | 2.51 | 2.48, 2.54
## 9 | 2.56 | 2.52, 2.60
## 10 | 2.62 | 2.57, 2.67ggplot(mydf, aes(x, predicted, colour = group)) + geom_line()