Study 4 analyses

In this report, we reproduce the analyses testing between group differences in Study 4.

prep data

First, we load the relevant packages, define functions and plotting aesthetics, and load and tidy the data.

load packages

if(!require('pacman')) {
    install.packages('pacman')
}

pacman::p_load(tidyverse, purrr, fs, knitr, lmerTest, ggeffects, parameters, kableExtra, boot, devtools, EMAtools, patchwork, install = TRUE)

define functions

# parameter estimate plotting function for purrr models
plot_model = function(model_data, palette, size = .75, facet = NULL, sharing_type = FALSE) {
  mod_renamed = model_data %>%
    mutate(term = gsub("msg_", "", term),
           term = gsub("_", " ", term),
           term = gsub(":", " x ", term),
           term = gsub(" z", "", term),
           term = gsub("rel self x topichealth", "self-relevance x\ntopic (health)", term),
           term = gsub("topichealth x rel social", "social relevance x\ntopic (health)", term),
           term = gsub("topichealth", "topic (health)", term),
           term = gsub("rel self", "self\nrelevance", term),
           term = gsub("rel social", "social\nrelevance", term),
           term = gsub(" within", "\nwithin", term),
           term = gsub(" between", "\nbetween", term),
           term = gsub("article condother", "other > control", term),
           term = gsub("article condself", "self > control", term),
           term = gsub("condition", "group", term),
           term = gsub("grouptimed", "group (timed)", term),
           term = gsub("groupuntimed", "group (untimed)", term),
           term = gsub("x topic \\(health\\)", "x\ntopic (health)", term)) 
  
  if (isTRUE(sharing_type)) {
    mod_renamed = mod_renamed %>%
      mutate(sharing_type = recode(sharing_type, "msg_share_broad" = "broadcast sharing",
                               "msg_share_narrow" = "narrowcast sharing"))
  }
  
  mod = mod_renamed %>%
    ggplot(aes(x = term, y = estimate)) +
    geom_pointrange(aes( ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = size, linewidth = size) +
    geom_hline(yintercept = 0, color = "grey50", linetype = "dotted") +
    coord_flip() +
    scale_fill_manual(name = "", values = palette) +
    scale_color_manual(name = "", values = palette) +
    labs(x = "", y = "\nstandardized  regression coefficient\n") +
    plot_aes
  
  if (!is.null(facet)) {
    mod + 
      facet_grid(~ get(facet))
  } else {
    mod
  }
}

# plot model predictions function
# parameter estimate plotting function for purrr models
plot_predicted = function(predicted_data, palette, size = .75, facet = NULL, sharing_type = FALSE) {

  if (isTRUE(sharing_type)) {
    predicted_data = predicted_data %>%
      mutate(sharing_type = recode(sharing_type, "msg_share_broad" = "broadcast sharing",
                               "msg_share_narrow" = "narrowcast sharing"))
  }
  
  mod = predicted_data %>%
    ggplot(aes(x = group, y = predicted)) +
    geom_pointrange(aes( ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = size, linewidth = size) +
    geom_hline(yintercept = 0, color = "grey50", linetype = "dotted") +
    coord_flip() +
    scale_fill_manual(name = "", values = palette) +
    scale_color_manual(name = "", values = palette) +
    labs(x = "", y = "\nstandardized  regression coefficient\n") +
    plot_aes
  
  if (!is.null(facet)) {
    mod + 
      facet_grid(facet)
  } else {
    mod
  }
}


# MLM results table function
table_model = function(model_data, sharing_type = FALSE, intercept = FALSE, spread = FALSE) {
  
  mod = model_data %>%
    rename("SE" = std.error,
           "t" = statistic,
           "p" = p.value) %>%
    select(-group, -effect) %>%
    mutate_at(vars(-contains("term"), -contains("value"), -contains("sharing_type"), -p), round, 2) %>%
    mutate(term = gsub("msg_", "", term),
           term = gsub("_", " ", term),
           term = gsub(":", " x ", term),
           term = gsub("z", "", term),
           term = gsub("topichealth", "topic (health)", term),
           term = gsub("rel self", "self-relevance", term),
           term = gsub("rel social", "social relevance", term),
           term = gsub("within", "within", term),
           term = gsub("between", "between", term),
           term = gsub("sharing type", "sharing type (narrowcast)", term),
           term = ifelse(grepl("between x ", term), "sharing type (narrowcast) x social relevance between", term),
           term = gsub("article condself", "self - control", term),
           term = gsub("article condother", "other - control", term),
           term = gsub("conditiontimed", "group (timed)", term),
           term = gsub("conditionuntimed", "group (untimed)", term),
           term = gsub("\\(Intercept\\)", "control", term),
           p = ifelse(p < .001, "< .001",
                      ifelse(p == 1, "1.000", gsub("0.(.*)", ".\\1", sprintf("%.3f", p)))),
           `b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", estimate, conf.low, conf.high))
  
  if (isTRUE(intercept)) {
    mod = mod %>%
      mutate(term = recode(term, "control" = "intercept"))
  }
  
    if (isTRUE(sharing_type)) {
    mod = mod %>%
      mutate(sharing_type = recode(sharing_type, "msg_share_broad" = "broadcast sharing",
                               "msg_share_narrow" = "narrowcast sharing")) %>%
      select(sharing_type, term, `b [95% CI]`, df, t, p)
      
    } else {
      
    mod = mod %>%
      select(term, `b [95% CI]`, df, t, p)
    }
  
  if (isTRUE(spread)) {
    mod %>%
      select(-df, -t, -p) %>%
      spread(`b [95% CI]`) %>%
      kable() %>%
      kableExtra::kable_styling()
    
  } else {
    mod %>%
      kable() %>%
      kableExtra::kable_styling()
  }
}

define aesthetics

palette_condition = c("self" = "#ee9b00",
                      "control" = "#bb3e03",
                      "other" = "#005f73")
palette_group = c("comment" = "#005f73",
                  "timed" = "#5F0F40",
                  "untimed" = "#D295BF")
palette_group_con = c("timed > comment" = "#5F0F40",
                      "untimed > comment" = "#D295BF")


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 data

data = read.csv("../data/study4_data.csv", stringsAsFactors = FALSE) %>%
  select(-sharing_type) %>%
  rename("sharing_type" = sharing_type_key) %>%
  group_by(sharing_type) %>%
  mutate(article_cond = ifelse(article_cond == "social", "other", article_cond),
         msg_share_z = scale(msg_share, center = TRUE, scale = 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)) %>%
  rename("condition" = group) %>%
  ungroup()

sub_conditions = data %>%
  select(SID, condition) %>%
  unique()

descriptives

group ns

Sample size by group

data %>%
  select(condition, SID) %>%
  unique() %>%
  group_by(condition) %>%
  summarize(n = n()) %>%
  kable() %>%
  kable_styling()

condition	n
comment	124
timed	157
untimed	167

ratings

Summarize means and SDs

data %>%
  select(-msg_share_z) %>%
  spread(sharing_type, msg_share) %>%
  gather(variable, value, msg_share_broad, msg_share_narrow, msg_rel_self, msg_rel_social) %>%
  group_by(variable) %>%
  summarize(M = mean(value, na.rm = TRUE),
            SD = sd(value, na.rm = TRUE)) %>%
  mutate(variable = recode(variable, "msg_rel_self" ="self-relevance",
                           "msg_rel_social" = "social relevance",
                           "msg_share_broad" = "broadcast intention",
                           "msg_share_narrow" = "narrowcast intention")) %>%
  kable(digits = 2) %>%
  kableExtra::kable_styling()

variable	M	SD
self-relevance	49.33	34.19
social relevance	54.47	33.04
broadcast intention	27.56	31.40
narrowcast intention	35.65	33.96

H1: sharing ~ self-relvance + social relevance

Greater (a) self-relevance and (b) social relevance ratings will be associated with stronger news sharing intentions and behavior.

run models

fit_mod = function(data){
  mod = lmerTest::lmer(msg_share_z ~ 1 + msg_rel_self_z * condition + msg_rel_social_z * condition +
                         (1 + msg_rel_self_z + msg_rel_social_z | SID) +
                         (1 | article_number), data = data,
                       control = lmerControl(optimizer = "bobyqa"))
  return(mod)
}
model_lmer = data  %>%
  group_by(sharing_type) %>%
  nest() %>%
  mutate(test = map(data, fit_mod))

model_data = model_lmer %>% 
  mutate(tidied = map(test, broom.mixed::tidy, conf.int = TRUE)) %>%
  select(-data, -test) %>%
  unnest(cols = tidied) %>%
  filter(effect == "fixed") %>%
  ungroup()

predicted_data = model_lmer %>% 
  mutate(`self-relevance` = map(test, ggeffects::ggpredict,
                                terms = c("msg_rel_self_z [-3,-2,-1,0,1,2,3]", "condition")),
         `social relevance` = map(test, ggeffects::ggpredict,
                                  terms = c("msg_rel_social_z [-3,-2,-1,0,1,2,3]", "condition"))) %>%
  select(-data, -test) %>%
  gather(variable, predicted, contains("relevance")) %>%
  unnest(cols = predicted) %>%
  mutate(sharing_type = recode(sharing_type, "msg_share_broad" = "broadcast sharing",
                               "msg_share_narrow" = "narrowcast sharing"))

predicted_random_data = model_lmer %>% 
  mutate(`self-relevance` = map(test, ggeffects::ggpredict,
                                terms = c("msg_rel_self_z [-3,-2,-1,0,1,2,3]", "condition", "SID"), type = "random"),
         `social relevance` = map(test, ggeffects::ggpredict,
                                  terms = c("msg_rel_social_z [-3,-2,-1,0,1,2,3]", "condition", "SID"), type = "random")) %>%
  select(-data, -test) %>%
  gather(variable, predicted, contains("relevance")) %>%
  unnest(cols = predicted) %>%
  mutate(sharing_type = recode(sharing_type, "msg_share_broad" = "broadcast sharing",
                               "msg_share_narrow" = "narrowcast sharing"))

model summary table

table_model(model_data, sharing_type = TRUE, intercept = TRUE)

sharing_type	term	b [95% CI]	df	t	p
broadcast sharing	intercept	-0.11 [-0.23, 0.01]	458.36	-1.88	.061
broadcast sharing	self-relevance	0.16 [0.10, 0.22]	331.72	5.29	< .001
broadcast sharing	group (timed)	0.05 [-0.11, 0.21]	442.85	0.64	.525
broadcast sharing	group (untimed)	0.08 [-0.07, 0.23]	443.26	1.02	.308
broadcast sharing	social relevance	0.23 [0.17, 0.29]	347.90	7.17	< .001
broadcast sharing	self-relevance x group (timed)	0.03 [-0.05, 0.11]	338.44	0.73	.463
broadcast sharing	self-relevance x group (untimed)	0.09 [0.01, 0.17]	334.40	2.32	.021
broadcast sharing	group (timed) x social relevance	0.03 [-0.05, 0.11]	347.04	0.70	.487
broadcast sharing	group (untimed) x social relevance	-0.03 [-0.11, 0.06]	345.11	-0.62	.537
narrowcast sharing	intercept	-0.02 [-0.12, 0.08]	451.88	-0.46	.649
narrowcast sharing	self-relevance	0.05 [-0.01, 0.12]	282.28	1.59	.113
narrowcast sharing	group (timed)	-0.03 [-0.16, 0.10]	434.90	-0.42	.675
narrowcast sharing	group (untimed)	-0.05 [-0.18, 0.07]	436.13	-0.85	.398
narrowcast sharing	social relevance	0.56 [0.49, 0.64]	360.20	15.16	< .001
narrowcast sharing	self-relevance x group (timed)	0.04 [-0.04, 0.13]	290.35	0.97	.332
narrowcast sharing	self-relevance x group (untimed)	0.11 [0.02, 0.20]	285.87	2.53	.012
narrowcast sharing	group (timed) x social relevance	-0.07 [-0.17, 0.03]	364.22	-1.46	.145
narrowcast sharing	group (untimed) x social relevance	-0.15 [-0.25, -0.05]	359.77	-3.07	.002

plot continuous relationships

(plot_a = predicted_data %>%
  ggplot(aes(x, predicted)) +
  stat_smooth(data = predicted_random_data, aes(group = interaction(facet, group), color = group),
              geom = "line", method = "lm", alpha = .15, linewidth = .1, se = FALSE) +
  geom_ribbon(aes(fill = group, ymin = conf.low, ymax = conf.high), alpha = .2) +
  geom_line(aes(color = group), size = 1) + 
  facet_grid(variable ~ sharing_type) +
  scale_fill_manual(name = "", values = palette_group) +
  scale_color_manual(name = "", values = palette_group) +
  #coord_cartesian(ylim = c(-2, 2.5)) +
  labs(x = "\nrelevance rating (SD)", y = "predicted sharing intention (SD)\n") + 
  plot_aes)

plot coefficients

plot_model(model_data, palette, facet = "sharing_type", sharing_type = TRUE, size = 1)

H2: relevance ~ intervention condition

Compared to the control condition, the (a) self-focused condition will increase self-relevance ratings, and (b) other-focused condition will increase social relevance ratings.

H2a self-relevance

run models

fit_mod = function(data){
  mod = lmerTest::lmer(msg_rel_self_z ~ 1 + article_cond * condition +
                         (1 | SID) +
                         (1 | article_number), data = data,
                       control = lmerControl(optimizer = "bobyqa"))
  return(mod)
}

model_lmer = data  %>%
  nest() %>%
  mutate(test = map(data, fit_mod))

model_data_self = model_lmer %>% 
  mutate(tidied = map(test, broom.mixed::tidy, conf.int = TRUE)) %>%
  select(-data, -test) %>%
  unnest(cols = tidied) %>%
  filter(effect == "fixed") %>%
  ungroup()

predicted_data_self = model_lmer %>% 
  mutate(predicted = map(test, ggeffects::ggpredict, terms = c("article_cond", "condition"))) %>%
  select(-data, -test) %>%
  unnest(cols = predicted)

model summary table

table_model(model_data_self)

term	b [95% CI]	df	t	p
control	0.01 [-0.12, 0.14]	607.40	0.20	.844
other - control	0.12 [0.05, 0.18]	10109.18	3.43	< .001
self - control	0.44 [0.37, 0.51]	10108.40	12.85	< .001
group (timed)	-0.10 [-0.25, 0.05]	562.50	-1.26	.208
group (untimed)	-0.12 [-0.27, 0.03]	563.72	-1.52	.130
other - control x group (timed)	-0.16 [-0.25, -0.07]	10105.42	-3.60	< .001
self - control x group (timed)	-0.31 [-0.40, -0.22]	10102.49	-6.82	< .001
other - control x group (untimed)	-0.13 [-0.22, -0.04]	10106.53	-2.89	.004
self - control x group (untimed)	-0.38 [-0.46, -0.29]	10108.99	-8.43	< .001

plot coefficients

plot_predicted(predicted_data_self, facet = "~ x", palette, size = 1)

H2b social relevance

run models

fit_mod = function(data){
  mod = lmerTest::lmer(msg_rel_social_z ~ 1 + article_cond * condition +
                         (1 | SID) +
                         (1 | article_number), data = data,
                       control = lmerControl(optimizer = "bobyqa"))
  return(mod)
}

model_lmer = data  %>%
  nest() %>%
  mutate(test = map(data, fit_mod))

model_data_social = model_lmer %>% 
  mutate(tidied = map(test, broom.mixed::tidy, conf.int = TRUE)) %>%
  select(-data, -test) %>%
  unnest(cols = tidied) %>%
  filter(effect == "fixed") %>%
  ungroup()

predicted_data_social = model_lmer %>% 
  mutate(predicted = map(test, ggeffects::ggpredict, terms = c("article_cond", "condition"))) %>%
  select(-data, -test) %>%
  unnest(cols = predicted)

model summary table

table_model(model_data_social)

term	b [95% CI]	df	t	p
control	-0.01 [-0.15, 0.12]	621.63	-0.21	.830
other - control	0.42 [0.36, 0.49]	10111.86	12.78	< .001
self - control	0.30 [0.23, 0.36]	10111.82	9.00	< .001
group (timed)	-0.10 [-0.26, 0.06]	539.49	-1.18	.237
group (untimed)	-0.13 [-0.29, 0.03]	540.45	-1.62	.106
other - control x group (timed)	-0.41 [-0.50, -0.32]	10109.14	-9.40	< .001
self - control x group (timed)	-0.21 [-0.30, -0.12]	10106.04	-4.81	< .001
other - control x group (untimed)	-0.31 [-0.39, -0.22]	10110.70	-7.15	< .001
self - control x group (untimed)	-0.23 [-0.31, -0.14]	10113.01	-5.30	< .001

plot coefficients

plot_predicted(predicted_data_social, facet = "~ x", palette, size = 1)

H3: sharing ~ intervention condition

Compared to the control condition, the (a) self-focused and (b) other-focused conditions will increase news sharing intentions and behavior.

run models

fit_mod = function(data){
  mod = lmerTest::lmer(msg_share_z ~ 1 + article_cond * condition +
                         (1 | SID) +
                         (1 | article_number), data = data,
                       control = lmerControl(optimizer = "bobyqa"))
  return(mod)
}

model_lmer = data  %>%
  group_by(sharing_type) %>%
  nest() %>%
  mutate(test = map(data, fit_mod))

model_data_share = model_lmer %>% 
  mutate(tidied = map(test, broom.mixed::tidy, conf.int = TRUE)) %>%
  select(-data, -test) %>%
  unnest(cols = tidied) %>%
  filter(effect == "fixed") %>%
  ungroup()

predicted_data_share = model_lmer %>% 
  mutate(predicted = map(test, ggeffects::ggpredict, terms = c("article_cond", "condition"))) %>%
  select(-data, -test) %>%
  unnest(cols = predicted)

model summary table

table_model(model_data_share, sharing_type = TRUE)

sharing_type	term	b [95% CI]	df	t	p
broadcast sharing	control	-0.09 [-0.25, 0.06]	545.83	-1.17	.244
broadcast sharing	other - control	0.21 [0.13, 0.28]	4812.68	5.55	< .001
broadcast sharing	self - control	0.18 [0.11, 0.26]	4814.83	4.90	< .001
broadcast sharing	group (timed)	0.06 [-0.15, 0.26]	520.52	0.54	.591
broadcast sharing	group (untimed)	0.06 [-0.14, 0.26]	521.24	0.59	.556
broadcast sharing	other - control x group (timed)	-0.15 [-0.25, -0.05]	4812.93	-3.06	.002
broadcast sharing	self - control x group (timed)	-0.18 [-0.27, -0.08]	4809.25	-3.60	< .001
broadcast sharing	other - control x group (untimed)	-0.18 [-0.27, -0.08]	4815.83	-3.60	< .001
broadcast sharing	self - control x group (untimed)	-0.16 [-0.25, -0.06]	4817.33	-3.21	.001
narrowcast sharing	control	-0.08 [-0.22, 0.07]	624.88	-1.05	.292
narrowcast sharing	other - control	0.42 [0.33, 0.51]	4809.39	9.60	< .001
narrowcast sharing	self - control	0.28 [0.19, 0.36]	4810.82	6.34	< .001
narrowcast sharing	group (timed)	0.02 [-0.17, 0.21]	575.78	0.18	.858
narrowcast sharing	group (untimed)	-0.01 [-0.19, 0.18]	577.00	-0.08	.937
narrowcast sharing	other - control x group (timed)	-0.39 [-0.51, -0.28]	4808.50	-6.81	< .001
narrowcast sharing	self - control x group (timed)	-0.25 [-0.36, -0.14]	4804.93	-4.32	< .001
narrowcast sharing	other - control x group (untimed)	-0.38 [-0.49, -0.27]	4811.05	-6.66	< .001
narrowcast sharing	self - control x group (untimed)	-0.28 [-0.39, -0.17]	4812.81	-4.92	< .001

plot coefficients

plot_predicted(predicted_data_share, facet = "x ~ sharing_type", palette, size = 1, sharing_type = TRUE)

plot contrasts

(plot_b = model_data_share %>%
  mutate(model = recode(sharing_type, "msg_share_broad" = "broadcast\nsharing",
                        "msg_share_narrow" = "narrowcast\nsharing")) %>%
  bind_rows(model_data_self %>% mutate(model = "self\nrelevance")) %>%
  bind_rows(model_data_social %>% mutate(model = "social\nrelevance")) %>%
  filter(!term == "(Intercept)") %>%
  mutate(term = recode(term, "article_condother" = "other > control",
                       "article_condself" = "self > control",
                       "article_condother:conditiontimed" = "other > control x group (timed > comment)",
                       "article_condself:conditiontimed" = "self > control x group (timed > comment)",
                       "article_condother:conditionuntimed" = "other > control x group (untimed > comment)",
                       "article_condself:conditionuntimed" = "self > control x group (untimed > comment)",                      
                       "conditiontimed" = "group (timed > comment)",
                       "conditionuntimed" = "group (untimed > comment)")) %>%
  filter(grepl("x", term)) %>%
  extract(term, c("term", "group"), "(.*) x group \\((.*)\\)") %>%
  mutate(term = factor(term, levels = c("self > control", "other > control"))) %>%
  ggplot(aes(x = model, y = estimate, color = group)) +
    geom_pointrange(aes( ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = 1, linewidth = 1) +
    geom_hline(yintercept = 0, color = "grey50", linetype = "dotted") +
    coord_flip() +
    facet_grid(~ term) +
    scale_fill_manual(name = "", values = palette_group_con) +
    scale_color_manual(name = "", values = palette_group_con) +
    labs(x = "", y = "\nstandardized  regression coefficient\n") +
    plot_aes)

combined H1-3 coefficient plot

(plot_a + labs(title = "H1: Correlational relationships") + plot_b + labs(title = "H2-3: Intervention effects")) +
  patchwork::plot_layout(ncol = 2) & theme(legend.position = 'top') & patchwork::plot_annotation(tag_levels = 'A')

cite packages

report::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>.
##   - 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>.
##   - Kleiman E (2021). _EMAtools: Data Management Tools for Real-Time Monitoring/Ecological Momentary Assessment Data_. R package version 0.1.4, <https://CRAN.R-project.org/package=EMAtools>.
##   - 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>.
##   - Lüdecke D, Ben-Shachar M, Patil I, Makowski D (2020). "Extracting, Computing and Exploring the Parameters of Statistical Models using R." _Journal of Open Source Software_, *5*(53), 2445. doi:10.21105/joss.02445 <https://doi.org/10.21105/joss.02445>.
##   - Müller K, Wickham H (2023). _tibble: Simple Data Frames_. R package version 3.2.1, <https://CRAN.R-project.org/package=tibble>.
##   - Pedersen T (2024). _patchwork: The Composer of Plots_. R package version 1.2.0, <https://CRAN.R-project.org/package=patchwork>.
##   - 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>.