This is the Prolific study with Veterans. Data collected in April/May
of 2026.
Load libraries
Clean data
T1_raw <- read_csv("data/Time 1.csv") |> slice(-c(1:2))
## New names:
## Rows: 32 Columns: 46
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (43): StartDate, EndDate, Status, IPAddress, Progress, Duration (in seco... lgl
## (3): ...44, ...45, ...46
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...44`
## • `` -> `...45`
## • `` -> `...46`
T2_raw <- read_csv("data/Time 2.csv") |> slice(-c(1:2))
## Rows: 26 Columns: 47
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (47): StartDate, EndDate, Status, IPAddress, Progress, Duration (in seco...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Clean helper
clean_wave <- function(df, suffix) {
df |>
filter(
Finished %in% c(TRUE, "TRUE", 1, "1"),
!is.na(prolific_ID),
prolific_ID != "{{%PROLIFIC_PID%}}"
) |>
mutate(
across(
c(PHQ_4_1, PHQ_4_2, PHQ_4_3, PHQ_4_4),
~ as.numeric(str_extract(as.character(.), "\\d+(\\.\\d+)?"))
),
# PHQ-4 items appear coded 1-4 in Qualtrics, so convert to 0-3
PHQ_4_1_recoded = PHQ_4_1 - 1,
PHQ_4_2_recoded = PHQ_4_2 - 1,
PHQ_4_3_recoded = PHQ_4_3 - 1,
PHQ_4_4_recoded = PHQ_4_4 - 1,
# PHQ-2 depression = first two PHQ-4 items
PHQ2_depression = PHQ_4_1_recoded + PHQ_4_2_recoded,
# GAD-2 anxiety = second two PHQ-4 items
GAD2_anxiety = PHQ_4_3_recoded + PHQ_4_4_recoded
) |>
select(
prolific_ID,
PHQ2_depression,
GAD2_anxiety,
everything()
) |>
rename_with(~ paste0(.x, "_", suffix), -prolific_ID)
}
T1 <- clean_wave(T1_raw, "T1")
T2 <- clean_wave(T2_raw, "T2")
merged <- inner_join(T1, T2, by = "prolific_ID")
Sample and retention
sample_flow <- tibble(
stage = c("Completed Time 1", "Completed Time 2", "Matched T1-T2"),
n = c(
n_distinct(T1$prolific_ID),
n_distinct(T2$prolific_ID),
n_distinct(merged$prolific_ID)
)
) |>
mutate(
percent_of_T1 = round(n / first(n) * 100, 1)
)
sample_flow |>
knitr::kable(digits = 1)
| Completed Time 1 |
30 |
100 |
| Completed Time 2 |
24 |
80 |
| Matched T1-T2 |
24 |
80 |
missing_T2 <- anti_join(
T1 |> select(prolific_ID),
T2 |> select(prolific_ID),
by = "prolific_ID"
)
missing_T2
## # A tibble: 6 × 1
## prolific_ID
## <chr>
## 1 600c167e28de48080e812884
## 2 69aa1de4c116eb95cf2e345b
## 3 69ebc155d394b6d05839fae7
## 4 5d227f901e99f80018c60e1f
## 5 69db0d18f4e828f3974d5586
## 6 69e4380a2de8698c069f7a3a
Demographics
demographics <- merged |>
transmute(
prolific_ID,
Age = as.numeric(Age_T1),
Sex = Sex_T1,
Gender = Gender_T1,
RaceEthnicity = RaceEthnicity_T1,
Education = Education_T1,
Therapy_before = Therapy_before_T1,
medication_before = medication_before_T1,
Device = Device_T1
) |>
distinct()
demographics |>
summarise(
n = n(),
mean_age = mean(Age, na.rm = TRUE),
sd_age = sd(Age, na.rm = TRUE),
min_age = min(Age, na.rm = TRUE),
max_age = max(Age, na.rm = TRUE)
) |>
knitr::kable(digits = 2)
make_percent_table <- function(data, var) {
data |>
count({{ var }}, name = "n") |>
mutate(percent = round(n / sum(n) * 100, 1)) |>
arrange(desc(n))
}
make_percent_table(demographics, Sex) |> knitr::kable()
| Male |
14 |
58.3 |
| Female |
10 |
41.7 |
make_percent_table(demographics, Gender) |> knitr::kable()
| Man |
15 |
62.5 |
| Woman |
8 |
33.3 |
| Genderqueer/Gender non-conforming,Gender
non-binary |
1 |
4.2 |
make_percent_table(demographics, RaceEthnicity) |> knitr::kable()
| White/Caucasian/European American |
17 |
70.8 |
| Black or African American |
3 |
12.5 |
| American Indian or Alaskan
Native,White/Caucasian/European American |
2 |
8.3 |
| Black or African American,Hispanic or Latina/e/o/x |
1 |
4.2 |
| Hispanic or Latina/e/o/x |
1 |
4.2 |
make_percent_table(demographics, Education) |> knitr::kable()
| Bachelor’s |
11 |
45.8 |
| Other (please specify): |
6 |
25.0 |
| Master’s |
5 |
20.8 |
| Associate’s |
1 |
4.2 |
| Non-degree student |
1 |
4.2 |
make_percent_table(demographics, Therapy_before) |> knitr::kable()
| Yes, in the past |
14 |
58.3 |
| Yes, currently |
6 |
25.0 |
| No |
4 |
16.7 |
make_percent_table(demographics, medication_before) |> knitr::kable()
| Yes, in the past |
12 |
50 |
| No |
6 |
25 |
| Yes, currently |
6 |
25 |
make_percent_table(demographics, Device) |> knitr::kable()
| android |
13 |
54.2 |
| ios |
11 |
45.8 |
Pre-post outcome summaries
outcomes <- c(
"PHQ2_depression",
"GAD2_anxiety"
)
prepost_summary <- map_dfr(outcomes, function(v) {
t1 <- merged[[paste0(v, "_T1")]]
t2 <- merged[[paste0(v, "_T2")]]
diff <- t2 - t1
test <- t.test(t2, t1, paired = TRUE)
d <- effectsize::cohens_d(t2, t1, paired = TRUE)
tibble(
outcome = v,
n = sum(!is.na(t1) & !is.na(t2)),
mean_T1 = mean(t1, na.rm = TRUE),
sd_T1 = sd(t1, na.rm = TRUE),
mean_T2 = mean(t2, na.rm = TRUE),
sd_T2 = sd(t2, na.rm = TRUE),
mean_change = mean(diff, na.rm = TRUE),
change_CI_low = test$conf.int[1],
change_CI_high = test$conf.int[2],
t = unname(test$statistic),
p = test$p.value,
cohens_d_paired = d$Cohens_d
)
})
## For paired samples, 'repeated_measures_d()' provides more options.
## For paired samples, 'repeated_measures_d()' provides more options.
prepost_summary_clean <- prepost_summary |>
mutate(
outcome = recode(
outcome,
"PHQ2_depression" = "PHQ-2 Depression",
"GAD2_anxiety" = "GAD-2 Anxiety"
),
outcome = factor(
outcome,
levels = c("PHQ-2 Depression", "GAD-2 Anxiety")
),
across(
c(mean_T1, sd_T1, mean_T2, sd_T2, mean_change, change_CI_low, change_CI_high, t, p, cohens_d_paired),
~ round(.x, 3)
)
) |>
arrange(outcome)
prepost_summary_clean |>
knitr::kable(
caption = "Pre-post changes among matched Time 1-Time 2 completers"
)
Pre-post changes among matched Time 1-Time 2
completers
| PHQ-2 Depression |
24 |
2.333 |
1.993 |
1.792 |
2.021 |
-0.542 |
-1.070 |
-0.014 |
-2.122 |
0.045 |
-0.433 |
| GAD-2 Anxiety |
24 |
2.708 |
2.255 |
2.208 |
2.085 |
-0.500 |
-1.028 |
0.028 |
-1.958 |
0.062 |
-0.400 |
Mixed-effects models
long_outcomes <- merged |>
dplyr::select(
prolific_ID,
PHQ2_depression_T1, PHQ2_depression_T2,
GAD2_anxiety_T1, GAD2_anxiety_T2
) |>
pivot_longer(
cols = -prolific_ID,
names_to = c("outcome", "time"),
names_pattern = "(.*)_T(1|2)",
values_to = "score"
) |>
mutate(
outcome = factor(
outcome,
levels = c("PHQ2_depression", "GAD2_anxiety")
),
time = as.numeric(time),
time_c = time - 1
)
model_phq2 <- lmer(
score ~ time_c + (1 | prolific_ID),
data = filter(long_outcomes, outcome == "PHQ2_depression")
)
model_gad2 <- lmer(
score ~ time_c + (1 | prolific_ID),
data = filter(long_outcomes, outcome == "GAD2_anxiety")
)
summary(model_phq2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c + (1 | prolific_ID)
## Data: filter(long_outcomes, outcome == "PHQ2_depression")
##
## REML criterion at convergence: 176.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72028 -0.55700 0.05565 0.30214 2.19125
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 3.2464 1.8018
## Residual 0.7817 0.8841
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.3333 0.4097 27.8867 5.696 4.23e-06 ***
## time_c -0.5417 0.2552 23.0000 -2.122 0.0448 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## time_c -0.311
summary(model_gad2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c + (1 | prolific_ID)
## Data: filter(long_outcomes, outcome == "GAD2_anxiety")
##
## REML criterion at convergence: 180.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02506 -0.48291 0.05541 0.47865 1.93130
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 3.9330 1.9832
## Residual 0.7826 0.8847
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.7083 0.4433 27.1287 6.110 1.55e-06 ***
## time_c -0.5000 0.2554 23.0000 -1.958 0.0625 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## time_c -0.288
Outcome plots
long_outcomes_plot <- long_outcomes |>
mutate(
outcome = factor(
outcome,
levels = c("PHQ2_depression", "GAD2_anxiety")
),
outcome_label = recode(
as.character(outcome),
"PHQ2_depression" = "PHQ-2 Depression",
"GAD2_anxiety" = "GAD-2 Anxiety"
),
outcome_label = factor(
outcome_label,
levels = c("PHQ-2 Depression", "GAD-2 Anxiety")
),
time_label = factor(
time,
levels = c(1, 2),
labels = c("Time 1", "Time 2")
)
)
Individual trajectories
ggplot(long_outcomes_plot, aes(x = time_label, y = score, group = prolific_ID)) +
geom_line(alpha = 0.20) +
geom_point(alpha = 0.35) +
stat_summary(
aes(group = 1),
fun = mean,
geom = "line",
linewidth = 1.2
) +
stat_summary(
aes(group = 1),
fun = mean,
geom = "point",
size = 3
) +
stat_summary(
aes(group = 1),
fun.data = mean_se,
geom = "errorbar",
width = 0.08,
linewidth = 0.8
) +
facet_wrap(~ outcome_label, scales = "free_y") +
labs(
x = NULL,
y = "Score",
title = "Pre-post symptom change among veterans using Flourish"
) +
theme_minimal(base_size = 13)
Mean change only
ggplot(long_outcomes_plot, aes(x = time_label, y = score, group = outcome_label)) +
stat_summary(fun = mean, geom = "line", linewidth = 1.1) +
stat_summary(fun = mean, geom = "point", size = 3) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.08) +
facet_wrap(~ outcome_label, scales = "free_y") +
labs(
x = NULL,
y = "Mean score",
title = "Mean pre-post symptom change"
) +
theme_minimal(base_size = 13)
Mean change plot
ggplot(long_outcomes_plot,
aes(x = time, y = score, color = outcome_label, group = outcome_label)) +
geom_jitter(
aes(group = outcome_label),
width = 0.08,
height = 0,
alpha = 0.25,
show.legend = FALSE
) +
geom_errorbar(
stat = "summary",
fun.data = mean_se,
width = 0.12,
linewidth = 0.8,
na.rm = TRUE
) +
geom_line(
stat = "summary",
fun = mean,
linewidth = 1,
na.rm = TRUE
) +
geom_point(
stat = "summary",
fun = mean,
size = 2.25,
na.rm = TRUE
) +
facet_grid(. ~ outcome_label) +
labs(
x = "Time",
y = "Symptom score",
color = "Outcome"
) +
scale_x_continuous(
breaks = c(1, 2),
labels = c("Time 1", "Time 2")
) +
scale_y_continuous(breaks = waiver()) +
scale_color_manual(
values = c(
"PHQ-2 Depression" = "#3A828E",
"GAD-2 Anxiety" = "#E9B12E"
)
) +
coord_cartesian(ylim = c(0, 6)) +
theme(
panel.background = element_rect(fill = "white"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "#7D7D7D", fill = NA, linewidth = 0.8),
axis.line = element_blank(),
axis.text = element_text(size = 13),
axis.text.x = element_text(margin = margin(t = 8)),
axis.text.y = element_text(margin = margin(r = 8)),
axis.title = element_text(size = 16, face = "bold"),
axis.title.x = element_text(margin = margin(t = 10)),
axis.title.y = element_text(margin = margin(r = 10)),
strip.background = element_blank(),
strip.text = element_text(size = 14, face = "bold"),
legend.position = "none"
)
Change scores
change_scores <- merged |>
transmute(
prolific_ID,
PHQ2_depression_change = PHQ2_depression_T2 - PHQ2_depression_T1,
GAD2_anxiety_change = GAD2_anxiety_T2 - GAD2_anxiety_T1
)
change_summary <- change_scores |>
summarise(
across(
ends_with("_change"),
list(
mean = ~ mean(.x, na.rm = TRUE),
sd = ~ sd(.x, na.rm = TRUE),
median = ~ median(.x, na.rm = TRUE),
min = ~ min(.x, na.rm = TRUE),
max = ~ max(.x, na.rm = TRUE)
)
)
) |>
pivot_longer(everything()) |>
separate(name, into = c("outcome", "stat"), sep = "_change_") |>
pivot_wider(names_from = stat, values_from = value) |>
mutate(
outcome = recode(
outcome,
"PHQ2_depression" = "PHQ-2 Depression",
"GAD2_anxiety" = "GAD-2 Anxiety"
),
outcome = factor(
outcome,
levels = c("PHQ-2 Depression", "GAD-2 Anxiety")
)
) |>
arrange(outcome)
change_summary |>
knitr::kable(digits = 2)
| PHQ-2 Depression |
-0.54 |
1.25 |
0 |
-4 |
1 |
| GAD-2 Anxiety |
-0.50 |
1.25 |
0 |
-4 |
1 |
change_scores_long <- change_scores |>
pivot_longer(
cols = -prolific_ID,
names_to = "outcome",
values_to = "change"
) |>
mutate(
outcome = recode(
outcome,
"PHQ2_depression_change" = "PHQ-2 Depression",
"GAD2_anxiety_change" = "GAD-2 Anxiety"
),
outcome = factor(
outcome,
levels = c("PHQ-2 Depression", "GAD-2 Anxiety")
)
)
ggplot(change_scores_long, aes(x = change)) +
geom_histogram(bins = 10, alpha = 0.75) +
geom_vline(xintercept = 0, linetype = "dashed") +
facet_wrap(~ outcome, scales = "free_x") +
labs(
x = "Change from Time 1 to Time 2",
y = "Number of participants",
title = "Distribution of individual symptom change scores",
subtitle = "Negative change indicates improvement"
) +
theme_minimal(base_size = 13)
Improvement / worsening counts
improvement_flags <- change_scores |>
mutate(
PHQ2_depression_status = case_when(
PHQ2_depression_change < 0 ~ "Improved",
PHQ2_depression_change == 0 ~ "No change",
PHQ2_depression_change > 0 ~ "Worsened",
TRUE ~ NA_character_
),
GAD2_anxiety_status = case_when(
GAD2_anxiety_change < 0 ~ "Improved",
GAD2_anxiety_change == 0 ~ "No change",
GAD2_anxiety_change > 0 ~ "Worsened",
TRUE ~ NA_character_
)
)
improvement_summary <- improvement_flags |>
select(ends_with("_status")) |>
pivot_longer(
everything(),
names_to = "outcome",
values_to = "status"
) |>
mutate(
outcome = str_remove(outcome, "_status"),
outcome = recode(
outcome,
"PHQ2_depression" = "PHQ-2 Depression",
"GAD2_anxiety" = "GAD-2 Anxiety"
),
outcome = factor(
outcome,
levels = c("PHQ-2 Depression", "GAD-2 Anxiety")
)
) |>
count(outcome, status) |>
group_by(outcome) |>
mutate(percent = round(n / sum(n) * 100, 1)) |>
ungroup() |>
arrange(outcome)
improvement_summary |>
knitr::kable()
| PHQ-2 Depression |
Improved |
9 |
37.5 |
| PHQ-2 Depression |
No change |
11 |
45.8 |
| PHQ-2 Depression |
Worsened |
4 |
16.7 |
| GAD-2 Anxiety |
Improved |
9 |
37.5 |
| GAD-2 Anxiety |
No change |
11 |
45.8 |
| GAD-2 Anxiety |
Worsened |
4 |
16.7 |
ggplot(improvement_summary, aes(x = outcome, y = percent, fill = status)) +
geom_col(position = "stack") +
coord_flip() +
labs(
x = NULL,
y = "Percent of participants",
fill = NULL,
title = "Percent improved, unchanged, or worsened"
) +
theme_minimal(base_size = 13)
Elevated baseline symptoms
For PHQ-2 and GAD-2, a commonly used elevated screening cutoff is 3
or higher.
cutoff_summary <- merged |>
transmute(
prolific_ID,
PHQ2_elevated_T1 = PHQ2_depression_T1 >= 3,
PHQ2_elevated_T2 = PHQ2_depression_T2 >= 3,
GAD2_elevated_T1 = GAD2_anxiety_T1 >= 3,
GAD2_elevated_T2 = GAD2_anxiety_T2 >= 3
)
cutoff_summary |>
summarise(
n = n(),
PHQ2_elevated_T1_n = sum(PHQ2_elevated_T1, na.rm = TRUE),
PHQ2_elevated_T1_percent = mean(PHQ2_elevated_T1, na.rm = TRUE) * 100,
PHQ2_elevated_T2_n = sum(PHQ2_elevated_T2, na.rm = TRUE),
PHQ2_elevated_T2_percent = mean(PHQ2_elevated_T2, na.rm = TRUE) * 100,
GAD2_elevated_T1_n = sum(GAD2_elevated_T1, na.rm = TRUE),
GAD2_elevated_T1_percent = mean(GAD2_elevated_T1, na.rm = TRUE) * 100,
GAD2_elevated_T2_n = sum(GAD2_elevated_T2, na.rm = TRUE),
GAD2_elevated_T2_percent = mean(GAD2_elevated_T2, na.rm = TRUE) * 100
) |>
knitr::kable(digits = 1)
| 24 |
8 |
33.3 |
5 |
20.8 |
11 |
45.8 |
10 |
41.7 |
cutoff_improvement <- cutoff_summary |>
summarise(
PHQ2_baseline_elevated_n = sum(PHQ2_elevated_T1, na.rm = TRUE),
PHQ2_moved_below_cutoff_n = sum(PHQ2_elevated_T1 == TRUE & PHQ2_elevated_T2 == FALSE, na.rm = TRUE),
PHQ2_moved_below_cutoff_percent = PHQ2_moved_below_cutoff_n / PHQ2_baseline_elevated_n * 100,
GAD2_baseline_elevated_n = sum(GAD2_elevated_T1, na.rm = TRUE),
GAD2_moved_below_cutoff_n = sum(GAD2_elevated_T1 == TRUE & GAD2_elevated_T2 == FALSE, na.rm = TRUE),
GAD2_moved_below_cutoff_percent = GAD2_moved_below_cutoff_n / GAD2_baseline_elevated_n * 100
)
cutoff_improvement |>
knitr::kable(digits = 1)
Meaningful symptom change
This section uses a reduction of 2 or more points as a meaningful
change threshold for PHQ-2 and GAD-2.
meaningful_change <- merged |>
transmute(
prolific_ID,
PHQ2_depression_T1,
PHQ2_depression_T2,
GAD2_anxiety_T1,
GAD2_anxiety_T2,
PHQ2_change = PHQ2_depression_T2 - PHQ2_depression_T1,
GAD2_change = GAD2_anxiety_T2 - GAD2_anxiety_T1,
PHQ2_reduced_2plus = PHQ2_change <= -2,
GAD2_reduced_2plus = GAD2_change <= -2
)
meaningful_change |>
summarise(
n = n(),
PHQ2_reduced_2plus_n = sum(PHQ2_reduced_2plus, na.rm = TRUE),
PHQ2_reduced_2plus_percent = mean(PHQ2_reduced_2plus, na.rm = TRUE) * 100,
GAD2_reduced_2plus_n = sum(GAD2_reduced_2plus, na.rm = TRUE),
GAD2_reduced_2plus_percent = mean(GAD2_reduced_2plus, na.rm = TRUE) * 100
) |>
knitr::kable(digits = 1)
meaningful_change |>
summarise(
PHQ2_elevated_baseline_n = sum(PHQ2_depression_T1 >= 3, na.rm = TRUE),
PHQ2_elevated_reduced_2plus_n = sum(PHQ2_depression_T1 >= 3 & PHQ2_reduced_2plus, na.rm = TRUE),
PHQ2_elevated_reduced_2plus_percent = PHQ2_elevated_reduced_2plus_n / PHQ2_elevated_baseline_n * 100,
GAD2_elevated_baseline_n = sum(GAD2_anxiety_T1 >= 3, na.rm = TRUE),
GAD2_elevated_reduced_2plus_n = sum(GAD2_anxiety_T1 >= 3 & GAD2_reduced_2plus, na.rm = TRUE),
GAD2_elevated_reduced_2plus_percent = GAD2_elevated_reduced_2plus_n / GAD2_elevated_baseline_n * 100
) |>
knitr::kable(digits = 1)
Acceptability ratings
acceptability <- merged |>
transmute(
prolific_ID,
useful = as.numeric(str_extract(as.character(useful_flourish_T2), "\\d+(\\.\\d+)?")),
engaging = as.numeric(str_extract(as.character(engaging_flourish_T2), "\\d+(\\.\\d+)?")),
chat = as.numeric(str_extract(as.character(chat_flourish_T2), "\\d+(\\.\\d+)?")),
bargain = as.numeric(str_extract(as.character(bargain_flourish_1_T2), "\\d+(\\.\\d+)?")),
expensive_but_wtp = as.numeric(str_extract(as.character(wtp_flourish_1_T2), "\\d+(\\.\\d+)?"))
)
acceptability_summary <- acceptability |>
summarise(
across(
-prolific_ID,
list(
n = ~ sum(!is.na(.x)),
mean = ~ mean(.x, na.rm = TRUE),
sd = ~ sd(.x, na.rm = TRUE),
median = ~ median(.x, na.rm = TRUE),
min = ~ min(.x, na.rm = TRUE),
max = ~ max(.x, na.rm = TRUE)
)
)
) |>
pivot_longer(everything()) |>
separate(name, into = c("variable", "stat"), sep = "_(?=[^_]+$)") |>
pivot_wider(names_from = stat, values_from = value)
acceptability_summary |>
knitr::kable(digits = 2)
| useful |
23 |
3.20 |
1.25 |
3 |
1 |
5 |
| engaging |
22 |
3.86 |
1.35 |
4 |
1 |
5 |
| chat |
22 |
3.66 |
1.37 |
4 |
1 |
5 |
| bargain |
24 |
4.83 |
3.16 |
5 |
0 |
10 |
| expensive_but_wtp |
24 |
10.79 |
7.03 |
10 |
0 |
30 |
acceptability_long <- acceptability |>
select(prolific_ID, useful, engaging, chat) |>
pivot_longer(
cols = -prolific_ID,
names_to = "rating_type",
values_to = "score"
) |>
mutate(
rating_type = recode(
rating_type,
useful = "Useful",
engaging = "Engaging",
chat = "Chat with Sunnie"
)
)
ggplot(acceptability_long, aes(x = rating_type, y = score)) +
geom_jitter(width = 0.08, alpha = 0.45, height = 0) +
stat_summary(fun = mean, geom = "point", size = 3) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.12) +
coord_cartesian(ylim = c(1, 5)) +
labs(
x = NULL,
y = "Rating",
title = "Acceptability ratings for Flourish"
) +
theme_minimal(base_size = 13)
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 5 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
Engagement metrics
engagement <- merged |>
transmute(
prolific_ID,
days_flourished = as.numeric(str_extract(as.character(metrics_flourish_1_T2), "\\d+(\\.\\d+)?")),
chats = as.numeric(str_extract(as.character(metrics_flourish_2_T2), "\\d+(\\.\\d+)?")),
weeklies = as.numeric(str_extract(as.character(metrics_flourish_3_T2), "\\d+(\\.\\d+)?")),
activities = as.numeric(str_extract(as.character(metrics_flourish_6_T2), "\\d+(\\.\\d+)?")),
badges = as.numeric(str_extract(as.character(metrics_flourish_7_T2), "\\d+(\\.\\d+)?")),
journal = as.numeric(str_extract(as.character(metrics_flourish_8_T2), "\\d+(\\.\\d+)?"))
)
engagement_summary <- engagement |>
summarise(
across(
-prolific_ID,
list(
n = ~ sum(!is.na(.x)),
mean = ~ mean(.x, na.rm = TRUE),
sd = ~ sd(.x, na.rm = TRUE),
median = ~ median(.x, na.rm = TRUE),
min = ~ min(.x, na.rm = TRUE),
max = ~ max(.x, na.rm = TRUE)
)
)
) |>
pivot_longer(everything()) |>
separate(name, into = c("metric", "stat"), sep = "_(?=[^_]+$)") |>
pivot_wider(names_from = stat, values_from = value)
engagement_summary |>
knitr::kable(digits = 2)
| days_flourished |
24 |
2.96 |
0.86 |
3.0 |
2 |
5 |
| chats |
24 |
1.58 |
1.56 |
1.0 |
0 |
6 |
| weeklies |
24 |
0.00 |
0.00 |
0.0 |
0 |
0 |
| activities |
24 |
1.42 |
1.82 |
0.5 |
0 |
6 |
| badges |
24 |
5.38 |
2.63 |
6.0 |
1 |
12 |
| journal |
24 |
2.12 |
1.60 |
2.0 |
1 |
8 |
engagement_long <- engagement |>
pivot_longer(
cols = -prolific_ID,
names_to = "metric",
values_to = "value"
) |>
mutate(
metric = recode(
metric,
days_flourished = "Days flourished",
chats = "Chats",
weeklies = "Weeklies",
activities = "Activities",
badges = "Badges",
journal = "Journal entries"
)
)
ggplot(engagement_long, aes(x = metric, y = value)) +
geom_boxplot(outlier.shape = NA, alpha = 0.6) +
geom_jitter(width = 0.12, alpha = 0.45) +
coord_flip() +
labs(
x = NULL,
y = "Count",
title = "Self-reported Flourish engagement"
) +
theme_minimal(base_size = 13)
Engagement and symptom change
This section tests whether specific engagement metrics predict
greater symptom improvement. For PHQ-2 and GAD-2, negative change scores
indicate improvement.
engagement_change <- merged |>
transmute(
prolific_ID,
PHQ2_depression_change = PHQ2_depression_T2 - PHQ2_depression_T1,
GAD2_anxiety_change = GAD2_anxiety_T2 - GAD2_anxiety_T1,
days_flourished = as.numeric(str_extract(as.character(metrics_flourish_1_T2), "\\d+(\\.\\d+)?")),
chats = as.numeric(str_extract(as.character(metrics_flourish_2_T2), "\\d+(\\.\\d+)?")),
activities = as.numeric(str_extract(as.character(metrics_flourish_6_T2), "\\d+(\\.\\d+)?"))
)
engagement_change |>
select(
days_flourished,
chats,
activities,
PHQ2_depression_change,
GAD2_anxiety_change
) |>
cor(use = "pairwise.complete.obs") |>
round(2)
## days_flourished chats activities PHQ2_depression_change
## days_flourished 1.00 0.41 0.54 -0.14
## chats 0.41 1.00 -0.24 0.21
## activities 0.54 -0.24 1.00 -0.22
## PHQ2_depression_change -0.14 0.21 -0.22 1.00
## GAD2_anxiety_change -0.10 0.20 -0.06 0.57
## GAD2_anxiety_change
## days_flourished -0.10
## chats 0.20
## activities -0.06
## PHQ2_depression_change 0.57
## GAD2_anxiety_change 1.00
# Change-score models
# Negative coefficients mean higher engagement is associated with greater symptom improvement.
lm(PHQ2_depression_change ~ days_flourished, data = engagement_change) |> summary()
##
## Call:
## lm(formula = PHQ2_depression_change ~ days_flourished, data = engagement_change)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4496 -0.7518 0.4459 0.6026 1.5504
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07617 0.94491 0.081 0.936
## days_flourished -0.20885 0.30724 -0.680 0.504
##
## Residual standard error: 1.265 on 22 degrees of freedom
## Multiple R-squared: 0.02057, Adjusted R-squared: -0.02395
## F-statistic: 0.462 on 1 and 22 DF, p-value: 0.5038
lm(PHQ2_depression_change ~ chats, data = engagement_change) |> summary()
##
## Call:
## lm(formula = PHQ2_depression_change ~ chats, data = engagement_change)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3582 -0.5653 0.2985 0.6418 1.8134
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.8134 0.3674 -2.214 0.0375 *
## chats 0.1716 0.1671 1.027 0.3156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.249 on 22 degrees of freedom
## Multiple R-squared: 0.04574, Adjusted R-squared: 0.002369
## F-statistic: 1.055 on 1 and 22 DF, p-value: 0.3156
lm(PHQ2_depression_change ~ activities, data = engagement_change) |> summary()
##
## Call:
## lm(formula = PHQ2_depression_change ~ activities, data = engagement_change)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5220 -0.6956 0.3253 0.6690 1.9363
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.3253 0.3254 -1.000 0.328
## activities -0.1527 0.1432 -1.067 0.298
##
## Residual standard error: 1.247 on 22 degrees of freedom
## Multiple R-squared: 0.0492, Adjusted R-squared: 0.005987
## F-statistic: 1.139 on 1 and 22 DF, p-value: 0.2975
lm(GAD2_anxiety_change ~ days_flourished, data = engagement_change) |> summary()
##
## Call:
## lm(formula = GAD2_anxiety_change ~ days_flourished, data = engagement_change)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4939 -0.5307 0.4324 0.5799 1.5061
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.06388 0.95043 -0.067 0.947
## days_flourished -0.14742 0.30904 -0.477 0.638
##
## Residual standard error: 1.273 on 22 degrees of freedom
## Multiple R-squared: 0.01024, Adjusted R-squared: -0.03475
## F-statistic: 0.2276 on 1 and 22 DF, p-value: 0.638
lm(GAD2_anxiety_change ~ chats, data = engagement_change) |> summary()
##
## Call:
## lm(formula = GAD2_anxiety_change ~ chats, data = engagement_change)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4060 -0.4463 0.4328 0.7552 1.7552
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.7552 0.3687 -2.048 0.0527 .
## chats 0.1612 0.1677 0.961 0.3469
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.253 on 22 degrees of freedom
## Multiple R-squared: 0.0403, Adjusted R-squared: -0.003324
## F-statistic: 0.9238 on 1 and 22 DF, p-value: 0.3469
lm(GAD2_anxiety_change ~ activities, data = engagement_change) |> summary()
##
## Call:
## lm(formula = GAD2_anxiety_change ~ activities, data = engagement_change)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5165 -0.5264 0.4440 0.5725 1.5626
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.44396 0.33335 -1.332 0.197
## activities -0.03956 0.14665 -0.270 0.790
##
## Residual standard error: 1.277 on 22 degrees of freedom
## Multiple R-squared: 0.003297, Adjusted R-squared: -0.04201
## F-statistic: 0.07277 on 1 and 22 DF, p-value: 0.7899
# Compact table of engagement-change associations
engagement_predictors <- c("days_flourished", "chats", "activities")
change_outcomes <- c("PHQ2_depression_change", "GAD2_anxiety_change")
engagement_change_models <- crossing(
outcome = change_outcomes,
predictor = engagement_predictors
) |>
mutate(
model = map2(outcome, predictor, ~ lm(as.formula(paste(.x, "~", .y)), data = engagement_change)),
summary = map(model, broom::tidy)
) |>
unnest(summary) |>
filter(term != "(Intercept)") |>
mutate(
outcome = recode(
outcome,
"PHQ2_depression_change" = "PHQ-2 Depression Change",
"GAD2_anxiety_change" = "GAD-2 Anxiety Change"
),
outcome = factor(
outcome,
levels = c("PHQ-2 Depression Change", "GAD-2 Anxiety Change")
),
predictor = recode(
predictor,
"days_flourished" = "Days flourished",
"chats" = "Chats",
"activities" = "Activities"
),
predictor = factor(
predictor,
levels = c("Days flourished", "Chats", "Activities")
),
across(c(estimate, std.error, statistic, p.value), ~ round(.x, 3))
) |>
arrange(outcome, predictor) |>
select(outcome, predictor, estimate, std.error, statistic, p.value)
engagement_change_models |>
knitr::kable(
caption = "Engagement metrics predicting symptom change"
)
Engagement metrics predicting symptom change
| PHQ-2 Depression Change |
Days flourished |
-0.209 |
0.307 |
-0.680 |
0.504 |
| PHQ-2 Depression Change |
Chats |
0.172 |
0.167 |
1.027 |
0.316 |
| PHQ-2 Depression Change |
Activities |
-0.153 |
0.143 |
-1.067 |
0.298 |
| GAD-2 Anxiety Change |
Days flourished |
-0.147 |
0.309 |
-0.477 |
0.638 |
| GAD-2 Anxiety Change |
Chats |
0.161 |
0.168 |
0.961 |
0.347 |
| GAD-2 Anxiety Change |
Activities |
-0.040 |
0.147 |
-0.270 |
0.790 |
# Mixed-effects moderation models
# These test whether each engagement metric moderates pre-post change.
# The key term is time_c:engagement_metric.
engagement_moderation_long <- long_outcomes |>
left_join(
engagement_change |>
select(prolific_ID, days_flourished, chats, activities),
by = "prolific_ID"
) |>
mutate(
days_flourished_z = as.numeric(scale(days_flourished)),
chats_z = as.numeric(scale(chats)),
activities_z = as.numeric(scale(activities))
)
model_phq2_days <- lmer(
score ~ time_c * days_flourished_z + (1 | prolific_ID),
data = filter(engagement_moderation_long, outcome == "PHQ2_depression")
)
model_phq2_chats <- lmer(
score ~ time_c * chats_z + (1 | prolific_ID),
data = filter(engagement_moderation_long, outcome == "PHQ2_depression")
)
model_phq2_activities <- lmer(
score ~ time_c * activities_z + (1 | prolific_ID),
data = filter(engagement_moderation_long, outcome == "PHQ2_depression")
)
model_gad2_days <- lmer(
score ~ time_c * days_flourished_z + (1 | prolific_ID),
data = filter(engagement_moderation_long, outcome == "GAD2_anxiety")
)
model_gad2_chats <- lmer(
score ~ time_c * chats_z + (1 | prolific_ID),
data = filter(engagement_moderation_long, outcome == "GAD2_anxiety")
)
model_gad2_activities <- lmer(
score ~ time_c * activities_z + (1 | prolific_ID),
data = filter(engagement_moderation_long, outcome == "GAD2_anxiety")
)
summary(model_phq2_days)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c * days_flourished_z + (1 | prolific_ID)
## Data: filter(engagement_moderation_long, outcome == "PHQ2_depression")
##
## REML criterion at convergence: 175.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.67994 -0.46495 -0.05313 0.29578 2.17584
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 3.0440 1.7447
## Residual 0.8004 0.8947
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.3333 0.4002 27.0446 5.830 3.29e-06 ***
## time_c -0.5417 0.2583 22.0000 -2.097 0.0477 *
## days_flourished_z 0.6644 0.4023 27.0446 1.651 0.1102
## time_c:days_flourished_z -0.1765 0.2596 22.0000 -0.680 0.5038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time_c dys_f_
## time_c -0.323
## dys_flrshd_ 0.000 0.000
## tm_c:dys_f_ 0.000 0.000 -0.323
summary(model_phq2_chats)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c * chats_z + (1 | prolific_ID)
## Data: filter(engagement_moderation_long, outcome == "PHQ2_depression")
##
## REML criterion at convergence: 176.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.66710 -0.53386 -0.07378 0.40437 2.13568
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 3.3821 1.8390
## Residual 0.7799 0.8831
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.33333 0.41643 26.50033 5.603 6.46e-06 ***
## time_c -0.54167 0.25493 22.00000 -2.125 0.0451 *
## chats_z 0.03662 0.41862 26.50033 0.087 0.9310
## time_c:chats_z 0.26317 0.25626 22.00000 1.027 0.3156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time_c chts_z
## time_c -0.306
## chats_z 0.000 0.000
## tm_c:chts_z 0.000 0.000 -0.306
summary(model_phq2_activities)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c * activities_z + (1 | prolific_ID)
## Data: filter(engagement_moderation_long, outcome == "PHQ2_depression")
##
## REML criterion at convergence: 175.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.75137 -0.49601 -0.06856 0.26362 2.24412
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 3.243 1.8009
## Residual 0.777 0.8815
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.3333 0.4093 26.6536 5.701 4.89e-06 ***
## time_c -0.5417 0.2545 22.0000 -2.129 0.0447 *
## activities_z 0.5341 0.4114 26.6536 1.298 0.2054
## time_c:activities_z -0.2729 0.2558 22.0000 -1.067 0.2975
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time_c actvt_
## time_c -0.311
## activitis_z 0.000 0.000
## tm_c:ctvts_ 0.000 0.000 -0.311
summary(model_gad2_days)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c * days_flourished_z + (1 | prolific_ID)
## Data: filter(engagement_moderation_long, outcome == "GAD2_anxiety")
##
## REML criterion at convergence: 178.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.99606 -0.38934 0.01001 0.48954 1.88647
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 3.6304 1.9054
## Residual 0.8098 0.8999
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.7083 0.4301 26.3709 6.297 1.08e-06 ***
## time_c -0.5000 0.2598 22.0000 -1.925 0.0673 .
## days_flourished_z 0.7329 0.4324 26.3709 1.695 0.1019
## time_c:days_flourished_z -0.1246 0.2611 22.0000 -0.477 0.6380
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time_c dys_f_
## time_c -0.302
## dys_flrshd_ 0.000 0.000
## tm_c:dys_f_ 0.000 0.000 -0.302
summary(model_gad2_chats)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c * chats_z + (1 | prolific_ID)
## Data: filter(engagement_moderation_long, outcome == "GAD2_anxiety")
##
## REML criterion at convergence: 178.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94493 -0.45567 0.06347 0.43907 1.89876
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 3.6754 1.9171
## Residual 0.7852 0.8861
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.7083 0.4311 26.2072 6.282 1.15e-06 ***
## time_c -0.5000 0.2558 22.0000 -1.955 0.0635 .
## chats_z 0.5241 0.4334 26.2072 1.209 0.2374
## time_c:chats_z 0.2472 0.2571 22.0000 0.961 0.3469
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time_c chts_z
## time_c -0.297
## chats_z 0.000 0.000
## tm_c:chts_z 0.000 0.000 -0.297
summary(model_gad2_activities)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c * activities_z + (1 | prolific_ID)
## Data: filter(engagement_moderation_long, outcome == "GAD2_anxiety")
##
## REML criterion at convergence: 181.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9900 -0.4800 0.0356 0.4607 1.9041
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 4.0955 2.024
## Residual 0.8155 0.903
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.70833 0.45235 25.95153 5.987 2.56e-06 ***
## time_c -0.50000 0.26069 22.00000 -1.918 0.0682 .
## activities_z 0.16298 0.45473 25.95153 0.358 0.7229
## time_c:activities_z -0.07069 0.26205 22.00000 -0.270 0.7899
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time_c actvt_
## time_c -0.288
## activitis_z 0.000 0.000
## tm_c:ctvts_ 0.000 0.000 -0.288
# Compact table of mixed-effects moderation results
moderation_models <- list(
"PHQ-2 Depression: Days flourished" = model_phq2_days,
"PHQ-2 Depression: Chats" = model_phq2_chats,
"PHQ-2 Depression: Activities" = model_phq2_activities,
"GAD-2 Anxiety: Days flourished" = model_gad2_days,
"GAD-2 Anxiety: Chats" = model_gad2_chats,
"GAD-2 Anxiety: Activities" = model_gad2_activities
)
extract_lmer_fixed <- function(model, model_name) {
coef_table <- summary(model)$coefficients
as.data.frame(coef_table) |>
rownames_to_column("term") |>
as_tibble() |>
rename(
estimate = Estimate,
std.error = `Std. Error`,
statistic = `t value`,
p.value = `Pr(>|t|)`
) |>
mutate(model = model_name)
}
moderation_summary <- imap_dfr(
moderation_models,
~ extract_lmer_fixed(.x, .y)
) |>
filter(str_detect(term, "time_c:")) |>
mutate(
model = factor(
model,
levels = names(moderation_models)
),
across(c(estimate, std.error, statistic, p.value), ~ round(.x, 3))
) |>
arrange(model) |>
select(model, term, estimate, std.error, statistic, p.value)
moderation_summary |>
knitr::kable(
caption = "Engagement moderation of pre-post symptom change"
)
Engagement moderation of pre-post symptom change
| PHQ-2 Depression: Days flourished |
time_c:days_flourished_z |
-0.176 |
0.260 |
-0.680 |
0.504 |
| PHQ-2 Depression: Chats |
time_c:chats_z |
0.263 |
0.256 |
1.027 |
0.316 |
| PHQ-2 Depression: Activities |
time_c:activities_z |
-0.273 |
0.256 |
-1.067 |
0.298 |
| GAD-2 Anxiety: Days flourished |
time_c:days_flourished_z |
-0.125 |
0.261 |
-0.477 |
0.638 |
| GAD-2 Anxiety: Chats |
time_c:chats_z |
0.247 |
0.257 |
0.961 |
0.347 |
| GAD-2 Anxiety: Activities |
time_c:activities_z |
-0.071 |
0.262 |
-0.270 |
0.790 |
# Plots: engagement metrics and PHQ-2 depression change
ggplot(engagement_change, aes(x = days_flourished, y = PHQ2_depression_change)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
x = "Days flourished",
y = "PHQ-2 depression change",
title = "Days flourished and PHQ-2 depression change",
subtitle = "Negative change indicates improvement"
) +
theme_minimal(base_size = 13)
## `geom_smooth()` using formula = 'y ~ x'
ggplot(engagement_change, aes(x = chats, y = PHQ2_depression_change)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
x = "Chats",
y = "PHQ-2 depression change",
title = "Chats and PHQ-2 depression change",
subtitle = "Negative change indicates improvement"
) +
theme_minimal(base_size = 13)
## `geom_smooth()` using formula = 'y ~ x'
ggplot(engagement_change, aes(x = activities, y = PHQ2_depression_change)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
x = "Activities",
y = "PHQ-2 depression change",
title = "Activities and PHQ-2 depression change",
subtitle = "Negative change indicates improvement"
) +
theme_minimal(base_size = 13)
## `geom_smooth()` using formula = 'y ~ x'
# Plots: engagement metrics and GAD-2 anxiety change
ggplot(engagement_change, aes(x = days_flourished, y = GAD2_anxiety_change)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
x = "Days flourished",
y = "GAD-2 anxiety change",
title = "Days flourished and GAD-2 anxiety change",
subtitle = "Negative change indicates improvement"
) +
theme_minimal(base_size = 13)
## `geom_smooth()` using formula = 'y ~ x'
ggplot(engagement_change, aes(x = chats, y = GAD2_anxiety_change)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
x = "Chats",
y = "GAD-2 anxiety change",
title = "Chats and GAD-2 anxiety change",
subtitle = "Negative change indicates improvement"
) +
theme_minimal(base_size = 13)
## `geom_smooth()` using formula = 'y ~ x'
ggplot(engagement_change, aes(x = activities, y = GAD2_anxiety_change)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
x = "Activities",
y = "GAD-2 anxiety change",
title = "Activities and GAD-2 anxiety change",
subtitle = "Negative change indicates improvement"
) +
theme_minimal(base_size = 13)
## `geom_smooth()` using formula = 'y ~ x'
Therapy and medication history as moderators
This section tests whether prior/current therapy or medication
history moderates pre-post symptom change.
therapy_med_moderation_long <- long_outcomes |>
left_join(
merged |>
select(prolific_ID, Therapy_before_T1, medication_before_T1),
by = "prolific_ID"
) |>
mutate(
Therapy_before_T1 = factor(
Therapy_before_T1,
levels = c("No", "Yes, in the past", "Yes, currently", "Prefer not to say")
),
medication_before_T1 = factor(
medication_before_T1,
levels = c("No", "Yes, in the past", "Yes, currently", "Prefer not to say")
)
)
# Therapy history moderation models
model_phq2_therapy <- lmer(
score ~ time_c * Therapy_before_T1 + (1 | prolific_ID),
data = filter(therapy_med_moderation_long, outcome == "PHQ2_depression")
)
model_gad2_therapy <- lmer(
score ~ time_c * Therapy_before_T1 + (1 | prolific_ID),
data = filter(therapy_med_moderation_long, outcome == "GAD2_anxiety")
)
summary(model_phq2_therapy)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c * Therapy_before_T1 + (1 | prolific_ID)
## Data: filter(therapy_med_moderation_long, outcome == "PHQ2_depression")
##
## REML criterion at convergence: 162
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4375 -0.4747 -0.1595 0.4853 2.1417
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 2.6349 1.6232
## Residual 0.7364 0.8581
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 1.500e+00 9.181e-01 2.607e+01
## time_c 6.410e-16 6.068e-01 2.100e+01
## Therapy_before_T1Yes, in the past 5.000e-01 1.041e+00 2.607e+01
## Therapy_before_T1Yes, currently 2.167e+00 1.185e+00 2.607e+01
## time_c:Therapy_before_T1Yes, in the past -9.286e-01 6.880e-01 2.100e+01
## time_c:Therapy_before_T1Yes, currently -7.401e-16 7.834e-01 2.100e+01
## t value Pr(>|t|)
## (Intercept) 1.634 0.114
## time_c 0.000 1.000
## Therapy_before_T1Yes, in the past 0.480 0.635
## Therapy_before_T1Yes, currently 1.828 0.079 .
## time_c:Therapy_before_T1Yes, in the past -1.350 0.192
## time_c:Therapy_before_T1Yes, currently 0.000 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time_c T__itp T__T1c t_:itp
## time_c -0.330
## Th__T1Y,itp -0.882 0.291
## Thrp__T1Y,c -0.775 0.256 0.683
## t_:T__T1itp 0.291 -0.882 -0.330 -0.226
## t_:T__T1Y,c 0.256 -0.775 -0.226 -0.330 0.683
summary(model_gad2_therapy)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c * Therapy_before_T1 + (1 | prolific_ID)
## Data: filter(therapy_med_moderation_long, outcome == "GAD2_anxiety")
##
## REML criterion at convergence: 169.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.80385 -0.36393 -0.02765 0.44770 1.81194
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 3.5306 1.879
## Residual 0.7902 0.889
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 2.000e+00 1.039e+00 2.518e+01
## time_c -7.051e-16 6.286e-01 2.100e+01
## Therapy_before_T1Yes, in the past 3.571e-01 1.178e+00 2.518e+01
## Therapy_before_T1Yes, currently 2.000e+00 1.342e+00 2.518e+01
## time_c:Therapy_before_T1Yes, in the past -7.857e-01 7.128e-01 2.100e+01
## time_c:Therapy_before_T1Yes, currently -1.667e-01 8.115e-01 2.100e+01
## t value Pr(>|t|)
## (Intercept) 1.924 0.0657 .
## time_c 0.000 1.0000
## Therapy_before_T1Yes, in the past 0.303 0.7643
## Therapy_before_T1Yes, currently 1.491 0.1485
## time_c:Therapy_before_T1Yes, in the past -1.102 0.2828
## time_c:Therapy_before_T1Yes, currently -0.205 0.8393
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time_c T__itp T__T1c t_:itp
## time_c -0.302
## Th__T1Y,itp -0.882 0.267
## Thrp__T1Y,c -0.775 0.234 0.683
## t_:T__T1itp 0.267 -0.882 -0.302 -0.207
## t_:T__T1Y,c 0.234 -0.775 -0.207 -0.302 0.683
# Medication history moderation models
model_phq2_medication <- lmer(
score ~ time_c * medication_before_T1 + (1 | prolific_ID),
data = filter(therapy_med_moderation_long, outcome == "PHQ2_depression")
)
model_gad2_medication <- lmer(
score ~ time_c * medication_before_T1 + (1 | prolific_ID),
data = filter(therapy_med_moderation_long, outcome == "GAD2_anxiety")
)
summary(model_phq2_medication)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c * medication_before_T1 + (1 | prolific_ID)
## Data: filter(therapy_med_moderation_long, outcome == "PHQ2_depression")
##
## REML criterion at convergence: 169.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.52784 -0.51270 -0.07328 0.35116 1.92875
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 3.1984 1.7884
## Residual 0.8393 0.9161
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 1.33333 0.82034 25.80674
## time_c -0.50000 0.52893 21.00000
## medication_before_T1Yes, in the past 1.25000 1.00470 25.80674
## medication_before_T1Yes, currently 1.50000 1.16013 25.80674
## time_c:medication_before_T1Yes, in the past 0.08333 0.64780 21.00000
## time_c:medication_before_T1Yes, currently -0.33333 0.74801 21.00000
## t value Pr(>|t|)
## (Intercept) 1.625 0.116
## time_c -0.945 0.355
## medication_before_T1Yes, in the past 1.244 0.225
## medication_before_T1Yes, currently 1.293 0.207
## time_c:medication_before_T1Yes, in the past 0.129 0.899
## time_c:medication_before_T1Yes, currently -0.446 0.660
##
## Correlation of Fixed Effects:
## (Intr) time_c m__itp m__T1c t_:itp
## time_c -0.322
## md__T1Y,itp -0.816 0.263
## mdct__T1Y,c -0.707 0.228 0.577
## t_:__T1Yitp 0.263 -0.816 -0.322 -0.186
## tm_:__T1Y,c 0.228 -0.707 -0.186 -0.322 0.577
summary(model_gad2_medication)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: score ~ time_c * medication_before_T1 + (1 | prolific_ID)
## Data: filter(therapy_med_moderation_long, outcome == "GAD2_anxiety")
##
## REML criterion at convergence: 172.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9777 -0.3958 -0.1047 0.4412 1.9267
##
## Random effects:
## Groups Name Variance Std.Dev.
## prolific_ID (Intercept) 4.0000 2.0000
## Residual 0.8036 0.8964
## Number of obs: 48, groups: prolific_ID, 24
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 1.667e+00 8.948e-01 2.480e+01
## time_c 4.487e-16 5.175e-01 2.100e+01
## medication_before_T1Yes, in the past 1.667e+00 1.096e+00 2.480e+01
## medication_before_T1Yes, currently 8.333e-01 1.265e+00 2.480e+01
## time_c:medication_before_T1Yes, in the past -7.500e-01 6.339e-01 2.100e+01
## time_c:medication_before_T1Yes, currently -5.000e-01 7.319e-01 2.100e+01
## t value Pr(>|t|)
## (Intercept) 1.863 0.0744 .
## time_c 0.000 1.0000
## medication_before_T1Yes, in the past 1.521 0.1409
## medication_before_T1Yes, currently 0.659 0.5162
## time_c:medication_before_T1Yes, in the past -1.183 0.2499
## time_c:medication_before_T1Yes, currently -0.683 0.5020
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time_c m__itp m__T1c t_:itp
## time_c -0.289
## md__T1Y,itp -0.816 0.236
## mdct__T1Y,c -0.707 0.205 0.577
## t_:__T1Yitp 0.236 -0.816 -0.289 -0.167
## tm_:__T1Y,c 0.205 -0.707 -0.167 -0.289 0.577
Compact therapy/medication moderation table
therapy_med_models <- list(
"PHQ-2 Depression: Therapy history" = model_phq2_therapy,
"PHQ-2 Depression: Medication history" = model_phq2_medication,
"GAD-2 Anxiety: Therapy history" = model_gad2_therapy,
"GAD-2 Anxiety: Medication history" = model_gad2_medication
)
extract_lmer_fixed <- function(model, model_name) {
coef_table <- summary(model)$coefficients
as.data.frame(coef_table) |>
tibble::rownames_to_column("term") |>
as_tibble() |>
rename(
estimate = Estimate,
std.error = `Std. Error`,
statistic = `t value`,
p.value = `Pr(>|t|)`
) |>
mutate(model = model_name)
}
therapy_med_moderation_summary <- imap_dfr(
therapy_med_models,
~ extract_lmer_fixed(.x, .y)
) |>
filter(str_detect(term, "time_c:")) |>
mutate(
model = factor(
model,
levels = names(therapy_med_models)
),
across(c(estimate, std.error, statistic, p.value), ~ round(.x, 3))
) |>
arrange(model, term) |>
select(model, term, estimate, std.error, statistic, p.value)
therapy_med_moderation_summary |>
knitr::kable(
caption = "Therapy and medication history moderation of pre-post symptom change"
)
Therapy and medication history moderation of pre-post symptom
change
| PHQ-2 Depression: Therapy history |
time_c:Therapy_before_T1Yes, currently |
0.000 |
0.783 |
0.000 |
1.000 |
| PHQ-2 Depression: Therapy history |
time_c:Therapy_before_T1Yes, in the past |
-0.929 |
0.688 |
-1.350 |
0.192 |
| PHQ-2 Depression: Medication history |
time_c:medication_before_T1Yes, currently |
-0.333 |
0.748 |
-0.446 |
0.660 |
| PHQ-2 Depression: Medication history |
time_c:medication_before_T1Yes, in the past |
0.083 |
0.648 |
0.129 |
0.899 |
| GAD-2 Anxiety: Therapy history |
time_c:Therapy_before_T1Yes, currently |
-0.167 |
0.812 |
-0.205 |
0.839 |
| GAD-2 Anxiety: Therapy history |
time_c:Therapy_before_T1Yes, in the past |
-0.786 |
0.713 |
-1.102 |
0.283 |
| GAD-2 Anxiety: Medication history |
time_c:medication_before_T1Yes, currently |
-0.500 |
0.732 |
-0.683 |
0.502 |
| GAD-2 Anxiety: Medication history |
time_c:medication_before_T1Yes, in the past |
-0.750 |
0.634 |
-1.183 |
0.250 |
Willingness to pay
wtp_data <- acceptability |>
select(prolific_ID, bargain, expensive_but_wtp)
wtp_data |>
summarise(
n_bargain = sum(!is.na(bargain)),
mean_bargain = mean(bargain, na.rm = TRUE),
sd_bargain = sd(bargain, na.rm = TRUE),
median_bargain = median(bargain, na.rm = TRUE),
n_expensive_but_wtp = sum(!is.na(expensive_but_wtp)),
mean_expensive_but_wtp = mean(expensive_but_wtp, na.rm = TRUE),
sd_expensive_but_wtp = sd(expensive_but_wtp, na.rm = TRUE),
median_expensive_but_wtp = median(expensive_but_wtp, na.rm = TRUE)
) |>
knitr::kable(digits = 2)
| 24 |
4.83 |
3.16 |
5 |
24 |
10.79 |
7.03 |
10 |
wtp_long <- wtp_data |>
pivot_longer(
cols = c(bargain, expensive_but_wtp),
names_to = "item",
values_to = "amount"
) |>
mutate(
item = recode(
item,
bargain = "Bargain / great deal",
expensive_but_wtp = "Expensive but willing to pay"
)
)
ggplot(wtp_long, aes(x = item, y = amount)) +
geom_boxplot(outlier.shape = NA, alpha = 0.6) +
geom_jitter(width = 0.10, alpha = 0.50) +
labs(
x = NULL,
y = "Monthly amount ($)",
title = "Monthly price perceptions for Flourish"
) +
theme_minimal(base_size = 13)
