Hey all–
On our call two (three?!) weeks ago:
So what I want to do here is just clarify what I did, especially since I was just following the R models that Katie wrote and saved in dropbox.
First, we’re going to load all the cleaned data
# Initial settings --------------------------------------------------------
library(tidyverse)
library(estimatr)
library(metafor)
library(emmeans)
library(magrittr)
library(here)
library(broom)
# Read data ---------------------------------------------------------------
t_data <- tibble(
nms = str_c(
here(),
"/data/cleaned/"
) %>%
dir() %>%
str_subset("study ") %>%
str_extract(
"study \\d"
),
tabs = str_c(
here(),
"/data/cleaned/"
) %>%
dir(full.names = T) %>%
str_subset("study ") %>%
map(readRDS)
)
The object t_data
should look like this:
## # A tibble: 5 × 2
## nms tabs
## <chr> <list>
## 1 study 1 <tibble [1,827 × 220]>
## 2 study 2 <tibble [1,739 × 46]>
## 3 study 3 <tibble [3,207 × 52]>
## 4 study 4 <tibble [2,789 × 63]>
## 5 study 5 <tibble [1,983 × 69]>
Now we estimate all the models I could find, including many which we’re not going to plot
# Study 3 -----------------------------------------------------------------
t_s3 <- tribble(
~mod_name, ~frm, ~tab,
"perc_acc_ovr", "out2 ~ group + party + male + age + white + college", t_data$tabs[[3]],
"perc_acc_ov_nct", "out2 ~ group ", t_data$tabs[[3]],
"perc_acc_reps", "out2 ~ group + male + age + white + college", t_data$tabs[[3]] %>%
filter(party == "rep"),
"perc_acc_reps_nct", "out2 ~ group", t_data$tabs[[3]] %>%
filter(party == "rep"),
"perc_acc_dems", "out2 ~ group + male + age + white + college", t_data$tabs[[3]] %>%
filter(party == "dem"),
"perc_acc_dems_nct", "out2 ~ group", t_data$tabs[[3]] %>%
filter(party == "dem"),
"change_ovr", "change ~ group + male + age + white + college", t_data$tabs[[3]],
"change_reps", "change ~ group + male + age + white + college", t_data$tabs[[3]] %>%
filter(party == "rep"),
"change_dems", "change ~ group + male + age + white + college", t_data$tabs[[3]] %>%
filter(party == "dem"),
# certainty analysis, exploratory
"cert_ovr", "cert2 ~ group + party + male + age + white + college", t_data$tabs[[3]],
"cert_ovr_nct", "cert2 ~ group", t_data$tabs[[3]],
"cert_reps", "cert2 ~ group + male + age + white + college", t_data$tabs[[3]] %>%
filter(party == "rep"),
"cert_reps_nct", "cert2 ~ group", t_data$tabs[[3]] %>%
filter(party == "rep"),
"cert_dems", "cert2 ~ group + male + age + white + college", t_data$tabs[[3]] %>%
filter(party == "dem"),
"cert_dems_nct", "cert2 ~ group", t_data$tabs[[3]] %>%
filter(party == "dem")
) %>%
mutate(
mods = frm %>%
map2(
tab,
\(i, j)
i %>%
as.formula %>%
lm_robust(
data = j
)
)
)
# Study 4 -----------------------------------------------------------------
t_s4 <- tribble(
~mod_name, ~frm, ~tab,
"perc_acc_ovr", "out2 ~ group + party + male + age + white + college", t_data$tabs[[4]],
"perc_acc_ovr_nct", "out2 ~ group", t_data$tabs[[4]],
"perc_acc_reps", "out2 ~ group + male + age + white + college", t_data$tabs[[4]] %>%
filter(party == "rep"),
"perc_acc_reps_nct", "out2 ~ group", t_data$tabs[[4]] %>%
filter(party == "rep"),
"perc_acc_dems", "out2 ~ group + male + age + white + college", t_data$tabs[[4]] %>%
filter(party == "dem"),
"perc_acc_dems_nct", "out2 ~ group", t_data$tabs[[4]] %>%
filter(party == "dem")
) %>%
mutate(
mods = frm %>%
map2(
tab,
\(i, j)
i %>%
as.formula %>%
lm_robust(
data = j
)
)
)
# Study 5 -----------------------------------------------------------------
t_s5 <- tribble(
~mod_name, ~frm, ~tab,
"perc_acc_ovr", "out2 ~ group + party + male + age + white + college", t_data$tabs[[5]],
"perc_acc_ovr_nct", "out2 ~ group", t_data$tabs[[5]],
"perc_acc_reps", "out2 ~ group + male + age + white + college", t_data$tabs[[5]] %>%
filter(party == "rep"),
"perc_acc_reps_nct", "out2 ~ group", t_data$tabs[[5]] %>%
filter(party == "rep"),
"perc_acc_dems", "out2 ~ group + male + age + white + college", t_data$tabs[[5]] %>%
filter(party == "dem"),
"perc_acc_dems_nct", "out2 ~ group", t_data$tabs[[5]] %>%
filter(party == "dem")
) %>%
mutate(
mods = frm %>%
map2(
tab,
\(i, j)
i %>%
as.formula %>%
lm_robust(
data = j
)
)
)
I guess I should pause here to ensure that Katie is satisfied I didn’t screw anything up.
So the central thing that we plot in this figure are the conditional expected values (ie, the point ranges in the :
The quantities are extracted from the lm_robust
objects
here
list(
t_s3,
t_s4,
t_s5
) %>%
map2(
"study " %>%
str_c(3:5),
\(i, j)
i %>%
filter(
mod_name %>%
equals("perc_acc_ovr")
) %>%
pluck("mods", 1) %>%
emmeans(~group) %>%
tidy %>%
select(1:2) %>%
mutate(
study = j
)
) %>%
list_rbind %>%
spread(
study, estimate
) %>%
slice(
c(
2:4, 5, 1
)
)
## # A tibble: 5 × 4
## group `study 3` `study 4` `study 5`
## <chr> <dbl> <dbl> <dbl>
## 1 control 5.13 4.88 5.17
## 2 early 5.23 NA NA
## 3 late 5.33 NA NA
## 4 standard NA 5.01 5.14
## 5 best NA 2.86 3.11
And then, the differences in the right facet are just different comparisons of these values.
lemme know if we’re on the same page!