Study 3 sensitivity analyses
Dani Cosme
2025-01-31
In this report, we reproduce the sensitivity analyses testing H4-6 in Study 3 in control and value ROIs.
prep data
First, we load the relevant packages, define functions and plotting aesthetics, and load and tidy the data.
load packages
define functions
source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")
# MLM results table function
table_model = function(model_data, print = TRUE) {
table = model_data %>%
broom.mixed::tidy(conf.int = TRUE) %>%
filter(effect == "fixed") %>%
rename("SE" = std.error,
"t" = statistic,
"p" = p.value) %>%
select(-group, -effect) %>%
mutate_at(vars(-contains("term"), -contains("p")), round, 2) %>%
mutate(term = gsub("article_cond", "", term),
term = gsub("\\(Intercept\\)", "intercept", term),
term = gsub("other$", "other - control", term),
term = gsub("self$", "self - control", term),
term = gsub("siteUSA", "sample (USA)", term),
term = gsub("value", "value", term),
term = gsub("value", "value", term),
term = gsub("msg_rel_self_z", "self-relevance", term),
term = gsub("msg_rel_social_z", "social relevance", term),
term = gsub("topichealth", "topic (health)", term),
term = gsub(":", " x ", term),
term = gsub("sample \\(USA\\) x social relevance", "social relevance x sample (USA)", term),
p = ifelse(p < .001, "< .001",
ifelse(p > .999, "1.000", gsub("0.(.*)", ".\\1", sprintf("%.3f", p)))),
`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", estimate, conf.low, conf.high)) %>%
select(term, `b [95% CI]`, df, t, p) %>%
arrange(term)
if (isTRUE(print)) {
table %>%
kable() %>%
kableExtra::kable_styling()
} else {
table
}
}
simple_slopes = function(model, var, moderator, continuous = TRUE) {
if (isTRUE(continuous)) {
emmeans::emtrends(model, as.formula(paste("~", moderator)), var = var) %>%
data.frame() %>%
rename("trend" = 2) %>%
mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", trend, asymp.LCL, asymp.UCL)) %>%
select(!!moderator, `b [95% CI]`) %>%
kable() %>%
kableExtra::kable_styling()
} else {
confint(emmeans::contrast(emmeans::emmeans(model, as.formula(paste("~", var, "|", moderator))), "revpairwise", by = moderator, adjust = "none")) %>%
data.frame() %>%
filter(grepl("control", contrast)) %>%
mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, asymp.LCL, asymp.UCL)) %>%
select(contrast, !!moderator, `b [95% CI]`) %>%
arrange(contrast) %>%
kable() %>%
kableExtra::kable_styling()
}
}define aesthetics
palette_condition = c("self" = "#ee9b00",
"control" = "#0a9396",
"other" = "#005f73")
palette_roi = c("value" = "#0a9396",
"auditory" = "grey50")
plot_aes = theme_minimal() +
theme(legend.position = "top",
legend.text = element_text(size = 12),
text = element_text(size = 16, family = "Futura Medium"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(color = "black"),
axis.line = element_line(colour = "black"),
axis.ticks.y = element_blank())load and tidy data
merged_all = read.csv("../data/study3_data_sensitivity.csv")
ratings_z = merged_all %>%
select(SID, trial, article_number, article_cond, msg_rel_self, msg_rel_social, msg_share) %>%
unique() %>%
mutate(msg_share_z = scale(msg_share, scale = TRUE, center = TRUE),
msg_rel_self_z = scale(msg_rel_self, center = TRUE, scale = TRUE),
msg_rel_social_z = scale(msg_rel_social, center = TRUE, scale = TRUE))
merged = merged_all %>%
filter(outlier == "no" | is.na(outlier)) %>%
group_by(SID, atlas) %>%
mutate(parameter_estimate_std = parameter_estimate / sd(parameter_estimate, na.rm = TRUE)) %>%
left_join(., ratings_z)
merged_wide = merged %>%
select(SID, site, trial, article_number, topic, article_cond, msg_share, msg_share_z,
msg_rel_self, msg_rel_self_z, msg_rel_social, msg_rel_social_z, atlas, parameter_estimate_std) %>%
spread(atlas, parameter_estimate_std)value ROI
H4: relevance ~ ROI activity
self-relevance
mod_h4a = lmer(msg_rel_self_z ~ value + (1 + value | SID),
data = merged_wide,
control = lmerControl(optimizer = "bobyqa"))model table
| term | b [95% CI] | df | t | p |
|---|---|---|---|---|
| intercept | 0.01 [-0.06, 0.09] | 84.57 | 0.3 | .764 |
| value | 0.04 [0.02, 0.07] | 83.01 | 3.6 | < .001 |
summary
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_rel_self_z ~ value + (1 + value | SID)
## Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 16577.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4223 -0.7039 0.1449 0.6835 2.4176
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## SID (Intercept) 0.1070261 0.32715
## value 0.0004833 0.02198 -0.83
## Residual 0.8910432 0.94395
## Number of obs: 6014, groups: SID, 85
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.01134 0.03765 84.57358 0.301 0.763913
## value 0.04440 0.01233 83.00829 3.602 0.000537 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## value -0.073plot
predicted_h4_value = ggeffects::ggpredict(mod_h4a, c("value [-4.5:5]")) %>%
data.frame() %>%
mutate(roi = "value",
variable = "self-relevance") %>%
bind_rows(ggeffects::ggpredict(mod_h4b, c("value [-4.5:5]")) %>%
data.frame() %>%
mutate(roi = "value",
variable = "social relevance"))
predicted_sub_h4_value = ggeffects::ggpredict(mod_h4a, terms = c("value [-4.5:5]", "SID"), type = "random") %>%
data.frame() %>%
mutate(roi = "value",
variable = "self-relevance") %>%
bind_rows(ggeffects::ggpredict(mod_h4b, c("value [-4.5:5]", "SID"), type = "random") %>%
data.frame() %>%
mutate(roi = "value",
variable = "social relevance"))
predicted_h4_value %>%
ggplot(aes(x, predicted)) +
stat_smooth(data = predicted_sub_h4_value, aes(group = group, color = roi), geom ='line', method = "lm", alpha = .1, size = 1, se = FALSE) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = roi), alpha = .5, color = NA) +
geom_line(aes(color = roi), size = 2) +
facet_grid(~variable) +
scale_color_manual(name = "", values = palette_roi, guide = FALSE) +
scale_fill_manual(name = "", values = palette_roi, guide = FALSE) +
labs(x = "\nROI activity (SD)", y = "predicted rating (SD)\n") +
plot_aes
H5: sharing ~ ROI activity
value ROI
mod_h5a = lmer(msg_share_z ~ value + (1 + value | SID),
data = merged_wide,
control = lmerControl(optimizer = "bobyqa"))model table
| term | b [95% CI] | df | t | p |
|---|---|---|---|---|
| intercept | 0.02 [-0.05, 0.10] | 83.50 | 0.59 | .560 |
| value | 0.07 [0.04, 0.10] | 83.33 | 4.71 | < .001 |
summary
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_share_z ~ value + (1 + value | SID)
## Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 16369
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5879 -0.7209 0.1236 0.7342 2.1967
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## SID (Intercept) 0.102154 0.31962
## value 0.006365 0.07978 -0.01
## Residual 0.887993 0.94233
## Number of obs: 5935, groups: SID, 85
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.02170 0.03707 83.49714 0.585 0.56
## value 0.07046 0.01495 83.32668 4.713 0.00000968 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## value 0.065plot
vals = seq(-4.5, 4.5, .1)
predicted_h5_value = ggeffects::ggpredict(mod_h5a, c("value [vals]")) %>%
data.frame() %>%
mutate(roi = "value",
variable = "narrowcast sharing")
predicted_sub_h5_value = ggeffects::ggpredict(mod_h5a, terms = c("value [vals]", "SID"), type = "random") %>%
data.frame() %>%
mutate(roi = "value",
variable = "narrowcast sharing")
predicted_h5_value %>%
ggplot(aes(x = x, y = predicted, color = roi, fill = roi)) +
stat_smooth(data = predicted_sub_h5_value, aes(group = group), geom ='line', method = "lm", alpha = .1, size = 1, se = FALSE) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2, color = NA) +
geom_line(size = 2) +
scale_color_manual(name = "", values = palette_roi) +
scale_fill_manual(name = "", values = palette_roi) +
labs(y = "predicted sharing intention (SD)\n", x = "\nROI activity (SD)") +
plot_aes +
theme(legend.position = "none")
H6 ROI activity ~ intervention condition
value ROI
mod_h6a = lmer(value ~ article_cond + (1 + article_cond | SID),
data = merged_wide,
control = lmerControl(optimizer = "bobyqa"))model table
| term | b [95% CI] | df | t | p |
|---|---|---|---|---|
| intercept | -0.26 [-0.38, -0.15] | 83.94 | -4.60 | < .001 |
| other - control | 0.02 [-0.06, 0.10] | 83.42 | 0.52 | .606 |
| self - control | -0.02 [-0.10, 0.06] | 83.69 | -0.52 | .603 |
summary
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: value ~ article_cond + (1 + article_cond | SID)
## Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 17300.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8328 -0.6379 0.0033 0.6476 3.7572
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## SID (Intercept) 0.23652 0.4863
## article_condother 0.04100 0.2025 -0.17
## article_condself 0.06259 0.2502 -0.06 0.57
## Residual 0.98337 0.9917
## Number of obs: 6014, groups: SID, 85
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.26309 0.05722 83.94166 -4.598 0.0000149 ***
## article_condother 0.01983 0.03826 83.42071 0.518 0.606
## article_condself -0.02162 0.04146 83.68796 -0.522 0.603
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) artcl_cndt
## artcl_cndth -0.314
## artcl_cndsl -0.242 0.522plot
predicted_h6_value = ggeffects::ggpredict(mod_h6a, c("article_cond")) %>%
data.frame() %>%
mutate(atlas = "value")
predicted_sub_h6_value = ggeffects::ggpredict(mod_h6a, terms = c("article_cond", "SID"), type = "random") %>%
data.frame() %>%
mutate(atlas = "value")
predicted_h6_value %>%
ggplot(aes(x = x, y = predicted)) +
stat_summary(data = predicted_sub_h6_value, aes(group = group), fun = "mean", geom = "line",
size = .1, color = "grey50") +
stat_summary(aes(group = group), fun = "mean", geom = "line", size = 1) +
geom_pointrange(aes(color = x, ymin = conf.low, ymax = conf.high), size = .75) +
scale_color_manual(name = "", values = palette_condition, guide = "none") +
scale_alpha_manual(name = "", values = c(1, .5)) +
labs(x = "", y = "predicted ROI activity (SD)\n") +
plot_aes
auditory ROI
H4: relevance ~ ROI activity
self-relevance
mod_h4a = lmer(msg_rel_self_z ~ auditory + (1 + auditory | SID),
data = merged_wide,
control = lmerControl(optimizer = "bobyqa"))model table
| term | b [95% CI] | df | t | p |
|---|---|---|---|---|
| auditory | -0.01 [-0.04, 0.02] | 83.27 | -0.77 | .443 |
| intercept | 0.00 [-0.07, 0.08] | 84.70 | 0.11 | .911 |
summary
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_rel_self_z ~ auditory + (1 + auditory | SID)
## Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 16592.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4090 -0.7083 0.1536 0.6690 2.3400
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## SID (Intercept) 0.107207 0.32743
## auditory 0.001519 0.03897 0.00
## Residual 0.892237 0.94458
## Number of obs: 6014, groups: SID, 85
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.004260 0.037952 84.701506 0.112 0.911
## auditory -0.009886 0.012832 83.270687 -0.770 0.443
##
## Correlation of Fixed Effects:
## (Intr)
## auditory -0.130plot
predicted_h4_auditory = ggeffects::ggpredict(mod_h4a, c("auditory [-4.5:5]")) %>%
data.frame() %>%
mutate(roi = "auditory",
variable = "self-relevance") %>%
bind_rows(ggeffects::ggpredict(mod_h4b, c("auditory [-4.5:5]")) %>%
data.frame() %>%
mutate(roi = "auditory",
variable = "social relevance"))
predicted_sub_h4_auditory = ggeffects::ggpredict(mod_h4a, terms = c("auditory [-4.5:5]", "SID"), type = "random") %>%
data.frame() %>%
mutate(roi = "auditory",
variable = "self-relevance") %>%
bind_rows(ggeffects::ggpredict(mod_h4b, c("auditory [-4.5:5]", "SID"), type = "random") %>%
data.frame() %>%
mutate(roi = "auditory",
variable = "social relevance"))
predicted_h4_auditory %>%
ggplot(aes(x, predicted)) +
stat_smooth(data = predicted_sub_h4_auditory, aes(group = group, color = roi), geom ='line', method = "lm", alpha = .1, size = 1, se = FALSE) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = roi), alpha = .5, color = NA) +
geom_line(aes(color = roi), size = 2) +
facet_grid(~variable) +
scale_color_manual(name = "", values = palette_roi, guide = FALSE) +
scale_fill_manual(name = "", values = palette_roi, guide = FALSE) +
labs(x = "\nROI activity (SD)", y = "predicted rating (SD)\n") +
plot_aes
H5: sharing ~ ROI activity
value ROI
mod_h5a = lmer(msg_share_z ~ auditory + (1 + auditory | SID),
data = merged_wide,
control = lmerControl(optimizer = "bobyqa"))model table
| term | b [95% CI] | df | t | p |
|---|---|---|---|---|
| auditory | -0.01 [-0.04, 0.01] | 81.68 | -1.00 | .321 |
| intercept | 0.01 [-0.07, 0.08] | 86.04 | 0.22 | .823 |
summary
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_share_z ~ auditory + (1 + auditory | SID)
## Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 16409.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5384 -0.7053 0.1130 0.7247 2.0442
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## SID (Intercept) 0.10084796 0.317566
## auditory 0.00007759 0.008809 0.36
## Residual 0.89944528 0.948391
## Number of obs: 5935, groups: SID, 85
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.008288 0.036963 86.038095 0.224 0.823
## auditory -0.012290 0.012306 81.682672 -0.999 0.321
##
## Correlation of Fixed Effects:
## (Intr)
## auditory -0.116plot
vals = seq(-4.5, 4.5, .1)
predicted_h5_auditory = ggeffects::ggpredict(mod_h5a, c("auditory [vals]")) %>%
data.frame() %>%
mutate(roi = "auditory",
variable = "narrowcast sharing")
predicted_sub_h5_auditory = ggeffects::ggpredict(mod_h5a, terms = c("auditory [vals]", "SID"), type = "random") %>%
data.frame() %>%
mutate(roi = "auditory",
variable = "narrowcast sharing")
predicted_h5_auditory %>%
ggplot(aes(x = x, y = predicted, color = roi, fill = roi)) +
stat_smooth(data = predicted_sub_h5_auditory, aes(group = group), geom ='line', method = "lm", alpha = .1, size = 1, se = FALSE) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2, color = NA) +
geom_line(size = 2) +
scale_color_manual(name = "", values = palette_roi) +
scale_fill_manual(name = "", values = palette_roi) +
labs(y = "predicted sharing intention (SD)\n", x = "\nROI activity (SD)") +
plot_aes +
theme(legend.position = "none")
H6 ROI activity ~ intervention condition
auditory ROI
mod_h6a = lmer(auditory ~ article_cond + (1 + article_cond | SID),
data = merged_wide,
control = lmerControl(optimizer = "bobyqa"))model table
| term | b [95% CI] | df | t | p |
|---|---|---|---|---|
| intercept | 0.45 [0.34, 0.56] | 84.23 | 8.14 | < .001 |
| other - control | -0.04 [-0.11, 0.02] | 83.94 | -1.28 | .205 |
| self - control | -0.02 [-0.09, 0.06] | 84.76 | -0.50 | .620 |
summary
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: auditory ~ article_cond + (1 + article_cond | SID)
## Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 17313.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8301 -0.6150 0.0111 0.6304 4.0809
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## SID (Intercept) 0.21736 0.4662
## article_condother 0.01199 0.1095 0.08
## article_condself 0.03259 0.1805 0.07 0.45
## Residual 0.99204 0.9960
## Number of obs: 6014, groups: SID, 85
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.44958 0.05525 84.22798 8.137 0.00000000000317 ***
## article_condother -0.04297 0.03363 83.93840 -1.278 0.205
## article_condself -0.01844 0.03707 84.75753 -0.497 0.620
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) artcl_cndt
## artcl_cndth -0.240
## artcl_cndsl -0.208 0.480plot
predicted_h6_auditory = ggeffects::ggpredict(mod_h6a, c("article_cond")) %>%
data.frame() %>%
mutate(atlas = "auditory")
predicted_sub_h6_auditory = ggeffects::ggpredict(mod_h6a, terms = c("article_cond", "SID"), type = "random") %>%
data.frame() %>%
mutate(atlas = "auditory")
predicted_h6_auditory %>%
ggplot(aes(x = x, y = predicted)) +
stat_summary(data = predicted_sub_h6_auditory, aes(group = group), fun = "mean", geom = "line",
size = .1, color = "grey50") +
stat_summary(aes(group = group), fun = "mean", geom = "line", size = 1) +
geom_pointrange(aes(color = x, ymin = conf.low, ymax = conf.high), size = .75) +
scale_color_manual(name = "", values = palette_condition, guide = "none") +
scale_alpha_manual(name = "", values = c(1, .5)) +
labs(x = "", y = "predicted ROI activity (SD)\n") +
plot_aes
combined figure
predicted = predicted_h4_value %>%
bind_rows(predicted_h5_value) %>%
bind_rows(predicted_h4_auditory) %>%
bind_rows(predicted_h5_auditory)
predicted_sub = predicted_sub_h4_value %>%
bind_rows(predicted_sub_h5_value) %>%
bind_rows(predicted_sub_h4_auditory) %>%
bind_rows(predicted_sub_h5_auditory)
predicted %>%
ggplot(aes(x, predicted, color = roi)) +
stat_smooth(data = predicted_sub, aes(group = interaction(group, roi)), geom ='line', method = "lm", alpha = .1, size = 1, se = FALSE) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = roi), alpha = .5, color = NA) +
geom_line(aes(color = roi), size = 2) +
facet_grid(~variable) +
scale_color_manual(name = "", values = palette_roi) +
scale_fill_manual(name = "", values = palette_roi) +
labs(x = "\nROI activity (SD)", y = "predicted rating (SD)\n") +
plot_aes
cite packages
## - Angelo Canty, B. D. Ripley (2024). _boot: Bootstrap R (S-Plus) Functions_. R package version 1.3-30. A. C. Davison, D. V. Hinkley (1997). _Bootstrap Methods and Their Applications_. Cambridge University Press, Cambridge. ISBN 0-521-57391-2, <doi:10.1017/CBO9780511802843>.
## - Bates D, Mächler M, Bolker B, Walker S (2015). "Fitting Linear Mixed-Effects Models Using lme4." _Journal of Statistical Software_, *67*(1), 1-48. doi:10.18637/jss.v067.i01 <https://doi.org/10.18637/jss.v067.i01>.
## - Bates D, Maechler M, Jagan M (2024). _Matrix: Sparse and Dense Matrix Classes and Methods_. R package version 1.7-0, <https://CRAN.R-project.org/package=Matrix>.
## - Bürkner P (2017). "brms: An R Package for Bayesian Multilevel Models Using Stan." _Journal of Statistical Software_, *80*(1), 1-28. doi:10.18637/jss.v080.i01 <https://doi.org/10.18637/jss.v080.i01>. Bürkner P (2018). "Advanced Bayesian Multilevel Modeling with the R Package brms." _The R Journal_, *10*(1), 395-411. doi:10.32614/RJ-2018-017 <https://doi.org/10.32614/RJ-2018-017>. Bürkner P (2021). "Bayesian Item Response Modeling in R with brms and Stan." _Journal of Statistical Software_, *100*(5), 1-54. doi:10.18637/jss.v100.i05 <https://doi.org/10.18637/jss.v100.i05>.
## - Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N, Ucar I, Bates D, Chambers J (2024). _Rcpp: Seamless R and C++ Integration_. R package version 1.0.13, <https://CRAN.R-project.org/package=Rcpp>. Eddelbuettel D, François R (2011). "Rcpp: Seamless R and C++ Integration." _Journal of Statistical Software_, *40*(8), 1-18. doi:10.18637/jss.v040.i08 <https://doi.org/10.18637/jss.v040.i08>. Eddelbuettel D (2013). _Seamless R and C++ Integration with Rcpp_. Springer, New York. doi:10.1007/978-1-4614-6868-4 <https://doi.org/10.1007/978-1-4614-6868-4>, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). "Extending R with C++: A Brief Introduction to Rcpp." _The American Statistician_, *72*(1), 28-36. doi:10.1080/00031305.2017.1375990 <https://doi.org/10.1080/00031305.2017.1375990>.
## - Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate." _Journal of Statistical Software_, *40*(3), 1-25. <https://www.jstatsoft.org/v40/i03/>.
## - Hester J, Wickham H, Csárdi G (2024). _fs: Cross-Platform File System Operations Based on 'libuv'_. R package version 1.6.4, <https://CRAN.R-project.org/package=fs>.
## - Kay M (2023). _tidybayes: Tidy Data and Geoms for Bayesian Models_. doi:10.5281/zenodo.1308151 <https://doi.org/10.5281/zenodo.1308151>, R package version 3.0.6, <http://mjskay.github.io/tidybayes/>.
## - Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest Package: Tests in Linear Mixed Effects Models." _Journal of Statistical Software_, *82*(13), 1-26. doi:10.18637/jss.v082.i13 <https://doi.org/10.18637/jss.v082.i13>.
## - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects from Regression Models." _Journal of Open Source Software_, *3*(26), 772. doi:10.21105/joss.00772 <https://doi.org/10.21105/joss.00772>.
## - Müller K, Wickham H (2023). _tibble: Simple Data Frames_. R package version 3.2.1, <https://CRAN.R-project.org/package=tibble>.
## - R Core Team (2024). _R: A Language and Environment for Statistical Computing_. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>.
## - Rinker TW, Kurkiewicz D (2018). _pacman: Package Management for R_. version 0.5.0, <http://github.com/trinker/pacman>.
## - Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_. Springer-Verlag New York. ISBN 978-3-319-24277-4, <https://ggplot2.tidyverse.org>.
## - Wickham H (2023). _forcats: Tools for Working with Categorical Variables (Factors)_. R package version 1.0.0, <https://CRAN.R-project.org/package=forcats>.
## - Wickham H (2023). _stringr: Simple, Consistent Wrappers for Common String Operations_. R package version 1.5.1, <https://CRAN.R-project.org/package=stringr>.
## - Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." _Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
## - Wickham H, Bryan J, Barrett M, Teucher A (2024). _usethis: Automate Package and Project Setup_. R package version 2.2.3, <https://CRAN.R-project.org/package=usethis>.
## - Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar of Data Manipulation_. R package version 1.1.4, <https://CRAN.R-project.org/package=dplyr>.
## - Wickham H, Henry L (2023). _purrr: Functional Programming Tools_. R package version 1.0.2, <https://CRAN.R-project.org/package=purrr>.
## - Wickham H, Hester J, Bryan J (2024). _readr: Read Rectangular Text Data_. R package version 2.1.5, <https://CRAN.R-project.org/package=readr>.
## - Wickham H, Hester J, Chang W, Bryan J (2022). _devtools: Tools to Make Developing R Packages Easier_. R package version 2.4.5, <https://CRAN.R-project.org/package=devtools>.
## - Wickham H, Vaughan D, Girlich M (2024). _tidyr: Tidy Messy Data_. R package version 1.3.1, <https://CRAN.R-project.org/package=tidyr>.
## - Xie Y (2024). _knitr: A General-Purpose Package for Dynamic Report Generation in R_. R package version 1.47, <https://yihui.org/knitr/>. Xie Y (2015). _Dynamic Documents with R and knitr_, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, <https://yihui.org/knitr/>. Xie Y (2014). "knitr: A Comprehensive Tool for Reproducible Research in R." In Stodden V, Leisch F, Peng RD (eds.), _Implementing Reproducible Computational Research_. Chapman and Hall/CRC. ISBN 978-1466561595.
## - Zhu H (2024). _kableExtra: Construct Complex Table with 'kable' and Pipe Syntax_. R package version 1.4.0, <https://CRAN.R-project.org/package=kableExtra>.
social
model table
summary