# =============================================================================
# STEP 1: DATA IMPORT AND VARIABLE CREATION
# Following Yarns et al. (2024) JAMA Network Open methods
# =============================================================================
# Load required packages
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(readxl)
library(psych)
##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(naniar)
library(readr)
library(dplyr)
options(scipen=999)
# -----------------------------------------------------------------------------
# 1.1 Import Data
# -----------------------------------------------------------------------------
setwd("/Users/melissalagunas/Desktop/Yarns Research Project")
library(readxl)
df_raw <- read_excel("EAET_vs_CBT_Study 2023-11-07.xlsx")
## New names:
## • `currentpain_b` -> `currentpain_b...37`
## • `currentpain_b` -> `currentpain_b...42`
# CRITICAL: Filter to only RANDOMIZED participants
# The dataset contains 150 rows, but only 126 were randomized
# Participants with missing treatment assignment (eaet = NA) were not randomized
df_raw <- df_raw %>%
filter(!is.na(eaet)) # Keep only those who were randomized
# Check dimensions and structure
cat("Dataset dimensions:", dim(df_raw)[1], "rows ×", dim(df_raw)[2], "columns\n")
## Dataset dimensions: 126 rows × 512 columns
cat("(After filtering to randomized participants only)\n\n")
## (After filtering to randomized participants only)
# View first few rows
head(df_raw)
## # A tibble: 6 × 512
## group id eaet dropout attend back_pain neck_pain leg_pain pelvic_pain
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 101 1 4 8 1 1 1 0
## 2 1 102 0 3 6 1 1 0 0
## 3 1 103 1 4 2 1 1 1 0
## 4 1 104 0 4 7 1 1 1 0
## 5 1 105 0 4 8 1 0 1 0
## 6 1 106 1 4 8 0 0 0 0
## # ℹ 503 more variables: tmj <dbl>, shoulder_pain <dbl>, arm_pain <dbl>,
## # fibro <dbl>, headache <dbl>, other_pain <dbl>, pain_duration <dbl>,
## # age <dbl>, female <dbl>, race_ethnicity <dbl>, marital_status <dbl>,
## # educ <dbl>, mmse <dbl>, med_cond <dbl>, psych_dx <dbl>, ptsd <dbl>,
## # medications <dbl>, opiates <dbl>, sat_therapy <dbl>, sat_therapist <dbl>,
## # sat_global <dbl>, wai <dbl>, pgic_p <dbl>, pgic_6m <dbl>,
## # worstpain_b <dbl>, leastpain_b <dbl>, avgpain_b <dbl>, …
# Check key variables exist
key_vars <- c("id", "group", "eaet", "dropout", "attend",
"race_ethnicity", "age", "female",
"wai", # Working Alliance Inventory - our MODERATOR
"meanpain_b", "meanpain_p", "meanpain_6m", # Primary outcome
"painint_b", "painint_p", "painint_6m", # Secondary outcome
"depress_b", "depress_p", "depress_6m") # Secondary outcome
missing_vars <- setdiff(key_vars, names(df_raw))
if(length(missing_vars) > 0) {
cat("WARNING: Missing variables:", paste(missing_vars, collapse = ", "), "\n")
} else {
cat("✓ All key variables present\n\n")
}
## ✓ All key variables present
# -----------------------------------------------------------------------------
# 1.2 Create Variables Following Paper's Methods
# -----------------------------------------------------------------------------
df <- df_raw %>%
mutate(
# ---------------------------------------------------------------------
# Treatment Variable
# ---------------------------------------------------------------------
# Original: eaet (0 = CBT, 1 = EAET)
# Create labeled factor for clarity
treatment = factor(eaet,
levels = c(0, 1),
labels = c("CBT", "EAET")),
# ---------------------------------------------------------------------
# Race/Ethnicity Variable - For Moderation Analysis
# ---------------------------------------------------------------------
# Based on Table 1 from paper:
# - 69 (55%) Black/African American
# - 39 (31%) White
# - 17 (14%) Multiracial or unknown
# - 8 (6%) Hispanic/Latino ethnicity
# From codebook: 1=Black, 2=White, 3=Hispanic, 4=Asian/PI, 5=Native Am, 6=Other
# White as reference group (standard practice)
race_white_vs_bipoc = case_when(
race_ethnicity == 2 ~ 0, # White = 0 (reference)
race_ethnicity %in% c(1, 3, 4, 5, 6) ~ 1, # BIPOC = 1
TRUE ~ NA_real_
),
# Create labels for reporting
race_group = factor(race_white_vs_bipoc,
levels = c(0, 1),
labels = c("White", "BIPOC")),
# ---------------------------------------------------------------------
# Completion Status
# ---------------------------------------------------------------------
# Based on dropout variable:
# 1 = dropped after individual, before group (baseline only)
# 2 = dropped during group (baseline only)
# 3 = dropped after group, before 6m (baseline + post)
# 4 = completed all assessments
completed_post = if_else(dropout >= 3, 1, 0), # Has post-treatment data
completed_6m = if_else(dropout == 4, 1, 0), # Has 6-month data
completer = if_else(dropout == 4, 1, 0), # Completed all
# ---------------------------------------------------------------------
# Attendance Variables
# ---------------------------------------------------------------------
# Paper reports: Mean 5.5 (SD 2.7) of 8 sessions attended
# attend already exists (0-8), create proportion
attend_prop = attend / 8, # Proportion of 8 sessions
# ---------------------------------------------------------------------
# Primary Outcome: Pain Severity (BPI)
# ---------------------------------------------------------------------
# Paper: "calculates the mean of worst, least, and average pain in the
# last week and current pain on scales ranging from 0 to 10"
# Note: meanpain_b, meanpain_p, meanpain_6m already calculated in dataset
# Calculate change scores (matching paper's approach)
pain_change_post = meanpain_p - meanpain_b,
pain_change_6m = meanpain_6m - meanpain_b,
# Calculate percent change
pain_percent_post = (pain_change_post / meanpain_b) * 100,
pain_percent_6m = (pain_change_6m / meanpain_b) * 100,
# Clinical improvement thresholds (from paper)
pain_improved_30pct_post = if_else(pain_percent_post <= -30, 1, 0),
pain_improved_50pct_post = if_else(pain_percent_post <= -50, 1, 0),
pain_improved_70pct_post = if_else(pain_percent_post <= -70, 1, 0),
pain_improved_30pct_6m = if_else(pain_percent_6m <= -30, 1, 0),
pain_improved_50pct_6m = if_else(pain_percent_6m <= -50, 1, 0),
pain_improved_70pct_6m = if_else(pain_percent_6m <= -70, 1, 0),
# ---------------------------------------------------------------------
# Secondary Outcomes: Pain Interference (PROMIS)
# ---------------------------------------------------------------------
painint_change_post = painint_p - painint_b,
painint_change_6m = painint_6m - painint_b,
# ---------------------------------------------------------------------
# Secondary Outcomes: Depression (PROMIS)
# ---------------------------------------------------------------------
depress_change_post = depress_p - depress_b,
depress_change_6m = depress_6m - depress_b,
# ---------------------------------------------------------------------
# Secondary Outcomes: Anxiety (PROMIS)
# ---------------------------------------------------------------------
anxiety_change_post = anxiety_p - anxiety_b,
anxiety_change_6m = anxiety_6m - anxiety_b,
# ---------------------------------------------------------------------
# Moderator: Working Alliance Inventory (WAI)
# ---------------------------------------------------------------------
# From codebook: "Overall score on the Working Alliance Inventory (WAI),
# obtained at post-treatment. Range 12-84 (sum of 12 items, each 1-7)"
# This is our KEY Moderator variable
# Note: Only measured at POST-treatment, not baseline
wai_present = if_else(!is.na(wai), 1, 0),
# Center WAI for interaction models (subtracts mean, keeps scale)
wai_c = as.numeric(scale(wai, center = TRUE, scale = FALSE)),
# Create interaction term for moderation analysis
wai_x_race = wai_c * race_white_vs_bipoc
)
# -----------------------------------------------------------------------------
# 1.3 Verify Variable Creation
# -----------------------------------------------------------------------------
cat("\n=== VARIABLE CREATION VERIFICATION ===\n\n")
##
## === VARIABLE CREATION VERIFICATION ===
# Check sample sizes match paper
cat("Total randomized:", nrow(df), "(paper reports 126)\n")
## Total randomized: 126 (paper reports 126)
cat("EAET:", sum(df$eaet == 1, na.rm=TRUE), "(paper reports 66)\n")
## EAET: 66 (paper reports 66)
cat("CBT:", sum(df$eaet == 0, na.rm=TRUE), "(paper reports 60)\n\n")
## CBT: 60 (paper reports 60)
# Check race/ethnicity distribution
cat("Race/Ethnicity Distribution:\n")
## Race/Ethnicity Distribution:
race_table <- table(df$race_group, useNA = "ifany")
print(race_table)
##
## White BIPOC
## 39 87
cat("\nWhite = 0 (reference), BIPOC = 1\n")
##
## White = 0 (reference), BIPOC = 1
cat("Paper reports: 69 (55%) Black, 39 (31%) White, 17 (14%) Multiracial/Unknown\n")
## Paper reports: 69 (55%) Black, 39 (31%) White, 17 (14%) Multiracial/Unknown
cat("Total with race data:", sum(race_table[1:2]), "\n")
## Total with race data: 126
if(sum(is.na(df$race_group)) > 0) {
cat("WARNING:", sum(is.na(df$race_group)), "participants missing race/ethnicity data\n")
}
cat("\n")
# Check completion rates
cat("Completion Rates:\n")
## Completion Rates:
cat("Post-treatment:", sum(df$completed_post, na.rm=T),
sprintf("(%.1f%%, paper reports 111, 88%%)\n",
mean(df$completed_post, na.rm=T)*100))
## Post-treatment: 111 (88.1%, paper reports 111, 88%)
cat("6-month follow-up:", sum(df$completed_6m, na.rm=T),
sprintf("(%.1f%%, paper reports 104, 82%%)\n\n",
mean(df$completed_6m, na.rm=T)*100))
## 6-month follow-up: 104 (82.5%, paper reports 104, 82%)
# Check attendance
cat("Attendance:\n")
## Attendance:
cat("Mean (SD):", sprintf("%.2f (%.2f) sessions\n",
mean(df$attend, na.rm=T),
sd(df$attend, na.rm=T)))
## Mean (SD): 5.50 (2.68) sessions
cat("Paper reports: 5.5 (2.7) sessions\n\n")
## Paper reports: 5.5 (2.7) sessions
# Check baseline pain severity
cat("Baseline Pain Severity (BPI):\n")
## Baseline Pain Severity (BPI):
df %>%
group_by(treatment) %>%
summarise(
n = n(),
Mean = mean(meanpain_b, na.rm=T),
SD = sd(meanpain_b, na.rm=T),
.groups = "drop"
) %>%
print()
## # A tibble: 2 × 4
## treatment n Mean SD
## <fct> <int> <dbl> <dbl>
## 1 CBT 60 6.23 1.63
## 2 EAET 66 5.97 1.98
cat("Paper reports: EAET 5.97(1.98), CBT 6.23(1.63)\n\n")
## Paper reports: EAET 5.97(1.98), CBT 6.23(1.63)
# Check WAI availability (critical for mediation!)
cat("Working Alliance Inventory (WAI) - Moderator VARIABLE:\n")
## Working Alliance Inventory (WAI) - Moderator VARIABLE:
cat("Available at post-treatment:", sum(!is.na(df$wai)), "participants\n")
## Available at post-treatment: 110 participants
cat("Missing:", sum(is.na(df$wai)), "participants\n")
## Missing: 16 participants
cat("This is expected - only measured post-treatment, not baseline\n\n")
## This is expected - only measured post-treatment, not baseline
# -----------------------------------------------------------------------------
# 1.4 Create Analysis-Ready Datasets
# -----------------------------------------------------------------------------
# Following paper's ITT approach: include all randomized participants
# But track who has data at each timepoint
# Full ITT sample (N=126)
df_itt <- df
# Sample with post-treatment data (for mediation analysis)
df_post <- df %>%
filter(!is.na(meanpain_p)) # Has post-treatment pain data
cat("=== ANALYSIS SAMPLES ===\n\n")
## === ANALYSIS SAMPLES ===
cat("Full ITT sample:", nrow(df_itt), "\n")
## Full ITT sample: 126
cat("Post-treatment sample:", nrow(df_post), "\n")
## Post-treatment sample: 111
cat(" - With WAI data:", sum(!is.na(df_post$wai)), "\n")
## - With WAI data: 110
cat(" - With complete data (pain + WAI):",
sum(!is.na(df_post$wai) & !is.na(df_post$meanpain_p)), "\n\n")
## - With complete data (pain + WAI): 110
# -----------------------------------------------------------------------------
# 1.5 Save Prepared Dataset
# -----------------------------------------------------------------------------
# Save for next steps
saveRDS(df_itt, "data_itt_prepared.rds")
saveRDS(df_post, "data_post_prepared.rds")
cat("✓ Data prepared and saved\n")
## ✓ Data prepared and saved
# Final verification
if(nrow(df_itt) == 126) {
cat("✓ Sample size matches paper (N=126)\n")
} else {
cat("⚠ WARNING: Sample size is", nrow(df_itt), "but paper reports 126\n")
}
## ✓ Sample size matches paper (N=126)
cat("✓ Ready for Step 2: Descriptive Statistics\n")
## ✓ Ready for Step 2: Descriptive Statistics
# =============================================================================
# STEP 2: DESCRIPTIVE STATISTICS & BASELINE EQUIVALENCE
# Focus: WAI as MODERATOR of treatment effects (moderated by race)
# =============================================================================
library(tidyverse)
library(psych)
library(effectsize)
##
## Attaching package: 'effectsize'
## The following object is masked from 'package:psych':
##
## phi
library(table1)
##
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
##
## units, units<-
# Load prepared data
df_itt <- readRDS("data_itt_prepared.rds")
df_post <- readRDS("data_post_prepared.rds")
cat("=== STEP 2: DESCRIPTIVE STATISTICS ===\n\n")
## === STEP 2: DESCRIPTIVE STATISTICS ===
cat("Analysis Samples:\n")
## Analysis Samples:
cat("- ITT sample (all randomized):", nrow(df_itt), "\n")
## - ITT sample (all randomized): 126
cat("- Post-treatment sample:", nrow(df_post), "\n")
## - Post-treatment sample: 111
cat("- With WAI data:", sum(!is.na(df_post$wai)), "\n\n")
## - With WAI data: 110
# =============================================================================
# 2.1 BASELINE CHARACTERISTICS BY TREATMENT GROUP (Table 1)
# =============================================================================
cat("=== TABLE 1: BASELINE CHARACTERISTICS ===\n\n")
## === TABLE 1: BASELINE CHARACTERISTICS ===
# Summary by treatment group
baseline_table <- df_itt %>%
group_by(treatment) %>%
summarise(
N = n(),
# Demographics
`Age M(SD)` = sprintf("%.1f (%.1f)", mean(age, na.rm=T), sd(age, na.rm=T)),
`Female n(%)` = sprintf("%d (%.0f%%)",
sum(female==1, na.rm=T),
mean(female==1, na.rm=T)*100),
`BIPOC n(%)` = sprintf("%d (%.0f%%)",
sum(race_white_vs_bipoc==1, na.rm=T),
mean(race_white_vs_bipoc==1, na.rm=T)*100),
# Pain characteristics
`Back pain n(%)` = sprintf("%d (%.0f%%)",
sum(back_pain==1, na.rm=T),
mean(back_pain==1, na.rm=T)*100),
`Pain duration M(SD)` = sprintf("%.1f (%.1f)",
mean(pain_duration, na.rm=T),
sd(pain_duration, na.rm=T)),
`Opiates n(%)` = sprintf("%d (%.0f%%)",
sum(opiates==1, na.rm=T),
mean(opiates==1, na.rm=T)*100),
# Clinical characteristics
`Psych dx n(%)` = sprintf("%d (%.0f%%)",
sum(psych_dx>0, na.rm=T),
mean(psych_dx>0, na.rm=T)*100),
`PTSD n(%)` = sprintf("%d (%.0f%%)",
sum(ptsd==1, na.rm=T),
mean(ptsd==1, na.rm=T)*100),
`Med conditions M(SD)` = sprintf("%.1f (%.1f)",
mean(med_cond, na.rm=T),
sd(med_cond, na.rm=T)),
`Medications M(SD)` = sprintf("%.1f (%.1f)",
mean(medications, na.rm=T),
sd(medications, na.rm=T)),
`MMSE M(SD)` = sprintf("%.1f (%.1f)",
mean(mmse, na.rm=T),
sd(mmse, na.rm=T)),
# Baseline outcomes
`Pain severity M(SD)` = sprintf("%.2f (%.2f)",
mean(meanpain_b, na.rm=T),
sd(meanpain_b, na.rm=T)),
`Pain interference M(SD)` = sprintf("%.1f (%.1f)",
mean(painint_b, na.rm=T),
sd(painint_b, na.rm=T)),
`Depression M(SD)` = sprintf("%.1f (%.1f)",
mean(depress_b, na.rm=T),
sd(depress_b, na.rm=T)),
`Anxiety M(SD)` = sprintf("%.1f (%.1f)",
mean(anxiety_b, na.rm=T),
sd(anxiety_b, na.rm=T)),
# Attendance
`Attendance M(SD)` = sprintf("%.1f (%.1f)",
mean(attend, na.rm=T),
sd(attend, na.rm=T)),
.groups = "drop"
)
print(baseline_table)
## # A tibble: 2 × 18
## treatment N `Age M(SD)` `Female n(%)` `BIPOC n(%)` `Back pain n(%)`
## <fct> <int> <chr> <chr> <chr> <chr>
## 1 CBT 60 71.2 (6.6) 3 (5%) 43 (72%) 58 (97%)
## 2 EAET 66 72.6 (5.2) 7 (11%) 44 (67%) 63 (95%)
## # ℹ 12 more variables: `Pain duration M(SD)` <chr>, `Opiates n(%)` <chr>,
## # `Psych dx n(%)` <chr>, `PTSD n(%)` <chr>, `Med conditions M(SD)` <chr>,
## # `Medications M(SD)` <chr>, `MMSE M(SD)` <chr>, `Pain severity M(SD)` <chr>,
## # `Pain interference M(SD)` <chr>, `Depression M(SD)` <chr>,
## # `Anxiety M(SD)` <chr>, `Attendance M(SD)` <chr>
# Save for manuscript
write.csv(baseline_table, "table1_baseline_characteristics.csv", row.names = FALSE)
cat("\n✓ Table saved to 'table1_baseline_characteristics.csv'\n\n")
##
## ✓ Table saved to 'table1_baseline_characteristics.csv'
# =============================================================================
# 2.2 TEST BASELINE EQUIVALENCE BETWEEN TREATMENT GROUPS
# =============================================================================
cat("\n=== BASELINE EQUIVALENCE TESTS ===\n\n")
##
## === BASELINE EQUIVALENCE TESTS ===
# Categorical variables - Chi-square tests
cat("Chi-Square Tests:\n")
## Chi-Square Tests:
cat("-----------------\n")
## -----------------
test_sex <- chisq.test(df_itt$treatment, df_itt$female)
## Warning in chisq.test(df_itt$treatment, df_itt$female): Chi-squared
## approximation may be incorrect
cat(sprintf("Sex: χ²(1) = %.2f, p = %.3f\n",
test_sex$statistic, test_sex$p.value))
## Sex: χ²(1) = 0.69, p = 0.405
test_race <- chisq.test(df_itt$treatment, df_itt$race_white_vs_bipoc)
cat(sprintf("Race (White vs BIPOC): χ²(1) = %.2f, p = %.3f\n",
test_race$statistic, test_race$p.value))
## Race (White vs BIPOC): χ²(1) = 0.17, p = 0.679
test_ptsd <- chisq.test(df_itt$treatment, df_itt$ptsd)
cat(sprintf("PTSD: χ²(1) = %.2f, p = %.3f\n\n",
test_ptsd$statistic, test_ptsd$p.value))
## PTSD: χ²(1) = 0.00, p = 1.000
# Continuous variables - t-tests with effect sizes
cat("Independent t-tests:\n")
## Independent t-tests:
cat("--------------------\n")
## --------------------
baseline_tests <- data.frame(
Variable = character(),
t = numeric(),
df = numeric(),
p = numeric(),
d = numeric(),
stringsAsFactors = FALSE
)
# Age
t_age <- t.test(age ~ treatment, data = df_itt)
d_age <- cohens_d(age ~ treatment, data = df_itt)$Cohens_d
cat(sprintf("Age: t(%.1f) = %.2f, p = %.3f, d = %.2f\n",
t_age$parameter, t_age$statistic, t_age$p.value, d_age))
## Age: t(111.8) = -1.24, p = 0.217, d = -0.22
baseline_tests <- rbind(baseline_tests,
data.frame(Variable="Age", t=t_age$statistic,
df=t_age$parameter, p=t_age$p.value, d=d_age))
# Pain severity
t_pain <- t.test(meanpain_b ~ treatment, data = df_itt)
d_pain <- cohens_d(meanpain_b ~ treatment, data = df_itt)$Cohens_d
cat(sprintf("Pain severity: t(%.1f) = %.2f, p = %.3f, d = %.2f\n",
t_pain$parameter, t_pain$statistic, t_pain$p.value, d_pain))
## Pain severity: t(122.8) = 0.81, p = 0.421, d = 0.14
baseline_tests <- rbind(baseline_tests,
data.frame(Variable="Pain severity", t=t_pain$statistic,
df=t_pain$parameter, p=t_pain$p.value, d=d_pain))
# Pain interference
t_painint <- t.test(painint_b ~ treatment, data = df_itt)
d_painint <- cohens_d(painint_b ~ treatment, data = df_itt)$Cohens_d
cat(sprintf("Pain interference: t(%.1f) = %.2f, p = %.3f, d = %.2f\n",
t_painint$parameter, t_painint$statistic, t_painint$p.value, d_painint))
## Pain interference: t(123.7) = 1.46, p = 0.146, d = 0.26
baseline_tests <- rbind(baseline_tests,
data.frame(Variable="Pain interference", t=t_painint$statistic,
df=t_painint$parameter, p=t_painint$p.value, d=d_painint))
# Depression
t_dep <- t.test(depress_b ~ treatment, data = df_itt)
d_dep <- cohens_d(depress_b ~ treatment, data = df_itt)$Cohens_d
cat(sprintf("Depression: t(%.1f) = %.2f, p = %.3f, d = %.2f\n",
t_dep$parameter, t_dep$statistic, t_dep$p.value, d_dep))
## Depression: t(119.4) = -0.20, p = 0.841, d = -0.04
baseline_tests <- rbind(baseline_tests,
data.frame(Variable="Depression", t=t_dep$statistic,
df=t_dep$parameter, p=t_dep$p.value, d=d_dep))
# Anxiety
t_anx <- t.test(anxiety_b ~ treatment, data = df_itt)
d_anx <- cohens_d(anxiety_b ~ treatment, data = df_itt)$Cohens_d
cat(sprintf("Anxiety: t(%.1f) = %.2f, p = %.3f, d = %.2f\n\n",
t_anx$parameter, t_anx$statistic, t_anx$p.value, d_anx))
## Anxiety: t(121.5) = 0.54, p = 0.593, d = 0.10
baseline_tests <- rbind(baseline_tests,
data.frame(Variable="Anxiety", t=t_anx$statistic,
df=t_anx$parameter, p=t_anx$p.value, d=d_anx))
# Summary of baseline equivalence
cat("Baseline Equivalence Summary:\n")
## Baseline Equivalence Summary:
sig_diffs <- sum(baseline_tests$p < 0.05)
if(sig_diffs == 0) {
cat("✓ No significant baseline differences detected (all p > .05)\n")
cat("✓ Groups are well-balanced for causal inference\n\n")
} else {
cat(sprintf("⚠ %d significant difference(s) detected\n", sig_diffs))
cat("Consider controlling for these variables in analyses\n\n")
}
## ✓ No significant baseline differences detected (all p > .05)
## ✓ Groups are well-balanced for causal inference
# =============================================================================
# 2.3 WORKING ALLIANCE INVENTORY (WAI) - MODERATOR
# =============================================================================
cat("\n=== WORKING ALLIANCE INVENTORY (WAI) ANALYSIS ===\n")
##
## === WORKING ALLIANCE INVENTORY (WAI) ANALYSIS ===
cat("(This is the MODERATOR variable)\n\n")
## (This is the MODERATOR variable)
# Overall WAI statistics
cat("Overall WAI Statistics:\n")
## Overall WAI Statistics:
wai_overall <- df_post %>%
summarise(
N = sum(!is.na(wai)),
Mean = mean(wai, na.rm=T),
SD = sd(wai, na.rm=T),
Min = min(wai, na.rm=T),
Max = max(wai, na.rm=T),
Median = median(wai, na.rm=T)
)
print(wai_overall)
## # A tibble: 1 × 6
## N Mean SD Min Max Median
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 110 63.4 13.4 19 84 65
# WAI by treatment group
cat("\nWAI by Treatment Group:\n")
##
## WAI by Treatment Group:
wai_by_tx <- df_post %>%
group_by(treatment) %>%
summarise(
N = sum(!is.na(wai)),
Mean = mean(wai, na.rm=T),
SD = sd(wai, na.rm=T),
.groups = "drop"
)
print(wai_by_tx)
## # A tibble: 2 × 4
## treatment N Mean SD
## <fct> <int> <dbl> <dbl>
## 1 CBT 53 62.8 11.8
## 2 EAET 57 63.9 14.8
# Test if WAI differs by treatment
t_wai_tx <- t.test(wai ~ treatment, data = df_post)
d_wai_tx <- cohens_d(wai ~ treatment, data = df_post)$Cohens_d
## Warning: Missing values detected. NAs dropped.
cat(sprintf("\nDoes WAI differ by treatment?\n"))
##
## Does WAI differ by treatment?
cat(sprintf("t(%.1f) = %.2f, p = %.3f, d = %.2f\n",
t_wai_tx$parameter, t_wai_tx$statistic, t_wai_tx$p.value, d_wai_tx))
## t(105.4) = -0.43, p = 0.671, d = -0.08
if(t_wai_tx$p.value < 0.05) {
cat("✓ Significant difference: Alliance differs between treatments\n")
} else {
cat("○ No significant difference in therapeutic alliance\n")
}
## ○ No significant difference in therapeutic alliance
# WAI by race/ethnicity
cat("\n\nWAI by Race/Ethnicity:\n")
##
##
## WAI by Race/Ethnicity:
wai_by_race <- df_post %>%
filter(!is.na(race_group)) %>%
group_by(race_group) %>%
summarise(
N = sum(!is.na(wai)),
Mean = mean(wai, na.rm=T),
SD = sd(wai, na.rm=T),
.groups = "drop"
)
print(wai_by_race)
## # A tibble: 2 × 4
## race_group N Mean SD
## <fct> <int> <dbl> <dbl>
## 1 White 34 62.9 14.1
## 2 BIPOC 76 63.6 13.2
# Test if WAI differs by race
t_wai_race <- t.test(wai ~ race_white_vs_bipoc, data = df_post)
d_wai_race <- cohens_d(wai ~ race_white_vs_bipoc, data = df_post)$Cohens_d
## Warning: Missing values detected. NAs dropped.
cat(sprintf("\nDoes WAI differ by race?\n"))
##
## Does WAI differ by race?
cat(sprintf("t(%.1f) = %.2f, p = %.3f, d = %.2f\n",
t_wai_race$parameter, t_wai_race$statistic, t_wai_race$p.value, d_wai_race))
## t(59.6) = -0.26, p = 0.793, d = -0.06
if(t_wai_race$p.value < 0.05) {
cat("✓ Significant difference: Alliance differs by race/ethnicity\n")
} else {
cat("○ No significant difference by race/ethnicity\n")
}
## ○ No significant difference by race/ethnicity
# WAI by Treatment × Race (2x2 table)
cat("\n\nWAI by Treatment × Race:\n")
##
##
## WAI by Treatment × Race:
wai_tx_race <- df_post %>%
filter(!is.na(race_group)) %>%
group_by(treatment, race_group) %>%
summarise(
N = sum(!is.na(wai)),
Mean = mean(wai, na.rm=T),
SD = sd(wai, na.rm=T),
.groups = "drop"
)
print(wai_tx_race)
## # A tibble: 4 × 5
## treatment race_group N Mean SD
## <fct> <fct> <int> <dbl> <dbl>
## 1 CBT White 15 61.6 12.0
## 2 CBT BIPOC 38 63.3 11.8
## 3 EAET White 19 63.8 15.9
## 4 EAET BIPOC 38 63.9 14.5
# Test for interaction (preview of moderation analysis)
cat("\nTest for Treatment × Race interaction on WAI:\n")
##
## Test for Treatment × Race interaction on WAI:
wai_interaction_model <- lm(wai ~ treatment * race_white_vs_bipoc,
data = df_post)
wai_int_summary <- summary(wai_interaction_model)
cat(sprintf("Interaction term: β = %.2f, SE = %.2f, t = %.2f, p = %.3f\n",
coef(wai_int_summary)["treatmentEAET:race_white_vs_bipoc", "Estimate"],
coef(wai_int_summary)["treatmentEAET:race_white_vs_bipoc", "Std. Error"],
coef(wai_int_summary)["treatmentEAET:race_white_vs_bipoc", "t value"],
coef(wai_int_summary)["treatmentEAET:race_white_vs_bipoc", "Pr(>|t|)"]))
## Interaction term: β = -1.61, SE = 5.62, t = -0.29, p = 0.775
if(coef(wai_int_summary)["treatmentEAET:race_white_vs_bipoc", "Pr(>|t|)"] < 0.05) {
cat("✓ Significant interaction: Treatment effect on alliance differs by race\n")
} else {
cat("○ No significant interaction\n")
}
## ○ No significant interaction
# =============================================================================
# 2.4 CORRELATIONS FOR MODERATION ANALYSIS
# =============================================================================
cat("\n\n=== CORRELATIONS AMONG KEY VARIABLES ===\n\n")
##
##
## === CORRELATIONS AMONG KEY VARIABLES ===
# Select key variables
cor_vars <- df_post %>%
select(
# Predictor
Treatment = eaet,
# Moderator
WAI = wai,
Race_BIPOC = race_white_vs_bipoc,
Attendance = attend,
# Outcomes at post-treatment
Pain_Change = pain_change_post,
PainInt_Change = painint_change_post,
Depress_Change = depress_change_post,
Anxiety_Change = anxiety_change_post,
# Baseline values
Pain_Baseline = meanpain_b,
Depress_Baseline = depress_b
)
# Correlation matrix
cor_matrix <- cor(cor_vars, use = "pairwise.complete.obs")
print(round(cor_matrix, 2))
## Treatment WAI Race_BIPOC Attendance Pain_Change
## Treatment 1.00 0.04 -0.06 -0.12 -0.37
## WAI 0.04 1.00 0.03 0.45 -0.20
## Race_BIPOC -0.06 0.03 1.00 0.00 -0.12
## Attendance -0.12 0.45 0.00 1.00 -0.18
## Pain_Change -0.37 -0.20 -0.12 -0.18 1.00
## PainInt_Change -0.09 -0.30 0.01 -0.24 0.55
## Depress_Change -0.25 -0.27 -0.19 -0.09 0.46
## Anxiety_Change -0.22 -0.29 -0.14 -0.16 0.48
## Pain_Baseline -0.11 -0.02 0.25 -0.01 -0.40
## Depress_Baseline 0.02 0.13 0.13 -0.01 -0.18
## PainInt_Change Depress_Change Anxiety_Change Pain_Baseline
## Treatment -0.09 -0.25 -0.22 -0.11
## WAI -0.30 -0.27 -0.29 -0.02
## Race_BIPOC 0.01 -0.19 -0.14 0.25
## Attendance -0.24 -0.09 -0.16 -0.01
## Pain_Change 0.55 0.46 0.48 -0.40
## PainInt_Change 1.00 0.50 0.47 -0.12
## Depress_Change 0.50 1.00 0.71 -0.27
## Anxiety_Change 0.47 0.71 1.00 -0.17
## Pain_Baseline -0.12 -0.27 -0.17 1.00
## Depress_Baseline -0.24 -0.50 -0.30 0.43
## Depress_Baseline
## Treatment 0.02
## WAI 0.13
## Race_BIPOC 0.13
## Attendance -0.01
## Pain_Change -0.18
## PainInt_Change -0.24
## Depress_Change -0.50
## Anxiety_Change -0.30
## Pain_Baseline 0.43
## Depress_Baseline 1.00
# Key correlations for moderation
cat("\n\nKey Correlations for WAI Moderation Analysis:\n")
##
##
## Key Correlations for WAI Moderation Analysis:
cat("----------------------------------------------\n")
## ----------------------------------------------
cat(sprintf("Treatment → WAI: r = %.3f (does treatment affect alliance?)\n",
cor_matrix["Treatment", "WAI"]))
## Treatment → WAI: r = 0.041 (does treatment affect alliance?)
cat(sprintf("WAI → Pain Change: r = %.3f (does alliance predict outcomes?)\n",
cor_matrix["WAI", "Pain_Change"]))
## WAI → Pain Change: r = -0.201 (does alliance predict outcomes?)
cat(sprintf("Treatment → Pain Change: r = %.3f (main treatment effect)\n",
cor_matrix["Treatment", "Pain_Change"]))
## Treatment → Pain Change: r = -0.368 (main treatment effect)
cat(sprintf("Race → WAI: r = %.3f (does alliance differ by race?)\n",
cor_matrix["Race_BIPOC", "WAI"]))
## Race → WAI: r = 0.026 (does alliance differ by race?)
cat(sprintf("Attendance → WAI: r = %.3f (does attendance relate to alliance?)\n",
cor_matrix["Attendance", "WAI"]))
## Attendance → WAI: r = 0.449 (does attendance relate to alliance?)
# Interpretation
cat("\nInterpretation for Moderation:\n")
##
## Interpretation for Moderation:
if(abs(cor_matrix["WAI", "Pain_Change"]) > 0.20) {
cat("✓ WAI correlates with pain outcomes - moderation analysis is viable\n")
} else {
cat("○ Weak correlation between WAI and outcomes - may limit moderation findings\n")
}
## ✓ WAI correlates with pain outcomes - moderation analysis is viable
# =============================================================================
# 2.5 VISUALIZATIONS
# =============================================================================
cat("\n\n=== CREATING VISUALIZATIONS ===\n\n")
##
##
## === CREATING VISUALIZATIONS ===
# Plot 1: WAI distribution by treatment
p1 <- ggplot(df_post, aes(x = treatment, y = wai, fill = treatment)) +
geom_boxplot(alpha = 0.7, outlier.shape = NA) +
geom_jitter(width = 0.2, alpha = 0.4, size = 2) +
stat_summary(fun = mean, geom = "point", shape = 23, size = 3,
fill = "red", color = "darkred") +
labs(title = "Working Alliance by Treatment",
subtitle = "Diamond = mean, box = median and IQR",
x = "Treatment",
y = "Working Alliance Inventory (WAI)") +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(face = "bold"))
print(p1)
## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
## Warning: Removed 1 rows containing non-finite values (`stat_summary()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).
# Plot 2: WAI by Treatment and Race
p2 <- ggplot(df_post %>% filter(!is.na(race_group)),
aes(x = treatment, y = wai, fill = race_group)) +
geom_boxplot(alpha = 0.7) +
stat_summary(fun = mean, geom = "point", shape = 23, size = 3,
position = position_dodge(width = 0.75)) +
labs(title = "Working Alliance by Treatment and Race/Ethnicity",
subtitle = "This interaction is key for your moderation analysis",
x = "Treatment",
y = "Working Alliance Inventory (WAI)",
fill = "Race/Ethnicity") +
theme_minimal() +
theme(plot.title = element_text(face = "bold"),
legend.position = "top")
print(p2)
## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
## Warning: Removed 1 rows containing non-finite values (`stat_summary()`).
# Plot 3: Relationship between WAI and Pain Change by Treatment
p3 <- ggplot(df_post, aes(x = wai, y = pain_change_post, color = treatment)) +
geom_point(alpha = 0.6, size = 2.5) +
geom_smooth(method = "lm", se = TRUE, linewidth = 1.2) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
labs(title = "Pain Change by Working Alliance and Treatment",
subtitle = "Testing if alliance moderates treatment effects",
x = "Working Alliance Inventory (WAI)",
y = "Pain Change (negative = improvement)",
color = "Treatment") +
theme_minimal() +
theme(plot.title = element_text(face = "bold"),
legend.position = "top")
print(p3)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).
# Plot 4: Relationship between WAI and Pain Change by Race
p4 <- ggplot(df_post %>% filter(!is.na(race_group)),
aes(x = wai, y = pain_change_post, color = race_group)) +
geom_point(alpha = 0.6, size = 2.5) +
geom_smooth(method = "lm", se = TRUE, linewidth = 1.2) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
labs(title = "Pain Change by Working Alliance and Race/Ethnicity",
subtitle = "Testing if alliance effects differ by race",
x = "Working Alliance Inventory (WAI)",
y = "Pain Change (negative = improvement)",
color = "Race/Ethnicity") +
theme_minimal() +
theme(plot.title = element_text(face = "bold"),
legend.position = "top")
print(p4)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).
## Removed 1 rows containing missing values (`geom_point()`).
# =============================================================================
# 2.6 SUMMARY AND READINESS FOR MODERATION ANALYSIS
# =============================================================================
cat("\n\n=== SUMMARY: READINESS FOR MODERATION ANALYSIS ===\n\n")
##
##
## === SUMMARY: READINESS FOR MODERATION ANALYSIS ===
# Check key assumptions for moderation
checks <- list(
baseline_balanced = sig_diffs == 0,
wai_available = sum(!is.na(df_post$wai)) >= 100,
wai_correlates = abs(cor_matrix["WAI", "Pain_Change"]) > 0.15,
race_data = sum(!is.na(df_post$race_white_vs_bipoc)) >= 100
)
cat("Checklist for Moderation Analysis:\n")
## Checklist for Moderation Analysis:
cat(sprintf("✓ Baseline equivalence: %s\n",
ifelse(checks$baseline_balanced, "PASS", "REVIEW")))
## ✓ Baseline equivalence: PASS
cat(sprintf("✓ Sufficient WAI data (N≥100): %s (N=%d)\n",
ifelse(checks$wai_available, "PASS", "FAIL"),
sum(!is.na(df_post$wai))))
## ✓ Sufficient WAI data (N≥100): PASS (N=110)
cat(sprintf("✓ WAI correlates with outcomes: %s (r=%.2f)\n",
ifelse(checks$wai_correlates, "PASS", "WEAK"),
cor_matrix["WAI", "Pain_Change"]))
## ✓ WAI correlates with outcomes: PASS (r=-0.20)
cat(sprintf("✓ Race/ethnicity data available: %s (N=%d)\n",
ifelse(checks$race_data, "PASS", "FAIL"),
sum(!is.na(df_post$race_white_vs_bipoc))))
## ✓ Race/ethnicity data available: PASS (N=111)
if(all(unlist(checks))) {
cat("\n✓✓✓ ALL CHECKS PASSED ✓✓✓\n")
cat("✓ Ready to proceed with moderation analysis!\n\n")
} else {
cat("\n⚠ Some checks need attention before proceeding\n\n")
}
##
## ✓✓✓ ALL CHECKS PASSED ✓✓✓
## ✓ Ready to proceed with moderation analysis!
# Save correlation matrix
write.csv(round(cor_matrix, 3), "correlation_matrix.csv")
cat("✓ Correlation matrix saved to 'correlation_matrix.csv'\n")
## ✓ Correlation matrix saved to 'correlation_matrix.csv'
cat("\n✓ Step 2 complete!\n")
##
## ✓ Step 2 complete!
cat("✓ Ready for Step 3: Moderation Analysis\n")
## ✓ Ready for Step 3: Moderation Analysis
# =============================================================================
# STEP 3: 3-WAY MODERATION ANALYSIS + SENSITIVITY (Attendance)
# - Treatment (CBT v EAET) x WAI (centered) x Race (White=0 [ref], BIPOC=1)
# - Mixed-effects models with random intercept for group
# - Outcomes: pain_change_post, painint_change_post, depress_change_post
# =============================================================================
library(tidyverse)
library(readxl)
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lmerTest)
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(interactions)
library(emmeans)
## Warning: package 'emmeans' was built under R version 4.3.3
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(broom.mixed)
## Warning: package 'broom.mixed' was built under R version 4.3.3
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools
library(car) # for Anova(type="III")
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:psych':
##
## logit
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
library(glue)
# -----------------------
# 1. Analysis setup
# -----------------------
outcomes <- list(
pain_change_post = "meanpain_b",
painint_change_post = "painint_b",
depress_change_post = "depress_b"
)
covariates <- c("age", "female", "attend_prop")
# -----------------------
# 2. Function to fit model and probe (prints to console / plots)
# -----------------------
fit_and_probe <- function(outcome, baseline_var, data, include_attendance = FALSE) {
dv <- outcome
covs <- c(baseline_var, covariates)
if (include_attendance) covs <- unique(c(covs, "attend"))
fixed_part <- paste0("treatment * wai_c * race_white_vs_bipoc + ", paste(covs, collapse = " + "))
fmla <- as.formula(paste0(dv, " ~ ", fixed_part, " + (1 | group)"))
message("\nFitting model for ", dv, ifelse(include_attendance, " (with attend)", " (main)"))
print(fmla)
mod <- lmer(fmla, data = data, REML = FALSE, control = lmerControl(optimizer="bobyqa"))
message("\n--- Summary ---")
print(summary(mod))
message("\n--- Type III ANOVA ---")
print(Anova(mod, type="III"))
tidy_fix <- broom.mixed::tidy(mod, effects = "fixed", conf.int = TRUE)
print("\n--- Fixed effects table ---")
print(tidy_fix)
# emmeans contrasts
wai_mean <- mean(data$wai_c, na.rm = TRUE)
wai_sd <- sd(data$wai_c, na.rm = TRUE)
probes <- c(wai_mean - wai_sd, wai_mean, wai_mean + wai_sd)
names(probes) <- c("low_wai", "mean_wai", "high_wai")
emm <- emmeans(mod, ~ treatment | race_white_vs_bipoc, at = list(wai_c = probes))
message("\n--- Estimated Marginal Means ---")
print(as.data.frame(emm))
contr <- contrast(emm, method = "pairwise", adjust = "none")
message("\n--- Pairwise contrasts ---")
print(as.data.frame(contr))
# interact_plot
message("\n--- Interaction plot ---")
tryCatch({
interactions::interact_plot(mod,
pred = "wai_c",
modx = "treatment",
mod2 = "race_white_vs_bipoc",
plot.points = TRUE,
interval = TRUE,
int.width = 0.95,
legend.main = "Treatment",
x.label = "WAI (centered)",
y.label = dv,
main.title = paste0("3-way: treatment × WAI × race on ", dv))
}, error = function(e) warning("interact_plot failed: ", e$message))
# simple slopes
for (r in c(0,1)) {
subdat <- data %>% filter(race_white_vs_bipoc == r)
if (nrow(subdat) < 10) {
message("Insufficient sample for race =", r)
} else {
tryCatch({
submod <- lmer(as.formula(paste0(dv, " ~ treatment * wai_c + ", paste(setdiff(covs, "race_white_vs_bipoc"), collapse = " + "), " + (1 | group)")),
data = subdat, REML = FALSE, control = lmerControl(optimizer="bobyqa"))
message("\n--- Simple slopes for race =", r, "---")
print(interactions::sim_slopes(submod, pred = "wai_c", modx = "treatment"))
}, error = function(e) message("sim_slopes failed for race =", r, ":", e$message))
}
}
return(list(model = mod, tidy_fixed = tidy_fix, emmeans = emm, contrasts = contr))
}
# -----------------------
# 3. Run main models
# -----------------------
main_results <- map2(names(outcomes), outcomes, ~ fit_and_probe(.x, .y, df_post, include_attendance = FALSE))
##
## Fitting model for pain_change_post (main)
## pain_change_post ~ treatment * wai_c * race_white_vs_bipoc +
## meanpain_b + age + female + attend_prop + (1 | group)
## <environment: 0x134c71e40>
##
## --- Summary ---
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: fmla
## Data: data
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 380.7 416.4 -176.3 352.7 81
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.21004 -0.61173 0.03432 0.62570 2.28741
##
## Random effects:
## Groups Name Variance Std.Dev.
## group (Intercept) 0.302 0.5496
## Residual 2.201 1.4834
## Number of obs: 95, groups: group, 10
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 4.73822 2.33651 93.37599 2.028
## treatmentEAET -1.30363 0.57870 89.86233 -2.253
## wai_c 0.03432 0.03748 91.99913 0.916
## race_white_vs_bipoc -0.03519 0.54054 88.35007 -0.065
## meanpain_b -0.43944 0.09531 94.07770 -4.611
## age -0.01805 0.03113 92.71804 -0.580
## female -0.96970 0.58504 91.83911 -1.658
## attend_prop -1.60329 0.70048 93.65939 -2.289
## treatmentEAET:wai_c -0.04395 0.04375 91.38556 -1.005
## treatmentEAET:race_white_vs_bipoc -0.59696 0.70963 92.05643 -0.841
## wai_c:race_white_vs_bipoc -0.02250 0.04527 89.80447 -0.497
## treatmentEAET:wai_c:race_white_vs_bipoc -0.01162 0.05313 89.04215 -0.219
## Pr(>|t|)
## (Intercept) 0.0454 *
## treatmentEAET 0.0267 *
## wai_c 0.3623
## race_white_vs_bipoc 0.9482
## meanpain_b 0.0000126 ***
## age 0.5634
## female 0.1008
## attend_prop 0.0243 *
## treatmentEAET:wai_c 0.3178
## treatmentEAET:race_white_vs_bipoc 0.4024
## wai_c:race_white_vs_bipoc 0.6204
## treatmentEAET:wai_c:race_white_vs_bipoc 0.8274
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trEAET wai_c rc_w__ mnpn_b age female attnd_ trEAET:_
## tretmntEAET -0.036
## wai_c 0.045 -0.074
## rc_wht_vs_b -0.251 0.592 -0.139
## meanpain_b -0.260 0.152 0.164 -0.075
## age -0.920 -0.182 -0.017 0.121 0.027
## female -0.142 0.069 -0.266 0.230 -0.208 0.148
## attend_prop -0.015 0.118 -0.171 -0.031 0.001 -0.233 0.029
## trtmnEAET:_ -0.095 0.008 -0.800 0.128 -0.156 0.148 0.198 -0.147
## trtEAET:___ 0.040 -0.806 0.095 -0.742 -0.105 0.113 -0.203 -0.025 -0.041
## w_c:rc_wh__ -0.042 0.061 -0.815 0.135 -0.108 0.019 0.210 0.116 0.659
## tEAET:_:___ 0.147 -0.011 0.655 -0.119 0.067 -0.160 -0.149 0.044 -0.798
## tEAET:__ w_:___
## tretmntEAET
## wai_c
## rc_wht_vs_b
## meanpain_b
## age
## female
## attend_prop
## trtmnEAET:_
## trtEAET:___
## w_c:rc_wh__ -0.099
## tEAET:_:___ 0.044 -0.826
##
## --- Type III ANOVA ---
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: pain_change_post
## Chisq Df Pr(>Chisq)
## (Intercept) 4.1124 1 0.04257 *
## treatment 5.0747 1 0.02428 *
## wai_c 0.8383 1 0.35988
## race_white_vs_bipoc 0.0042 1 0.94809
## meanpain_b 21.2574 1 0.000004016 ***
## age 0.3363 1 0.56199
## female 2.7473 1 0.09742 .
## attend_prop 5.2388 1 0.02209 *
## treatment:wai_c 1.0090 1 0.31513
## treatment:race_white_vs_bipoc 0.7077 1 0.40022
## wai_c:race_white_vs_bipoc 0.2470 1 0.61918
## treatment:wai_c:race_white_vs_bipoc 0.0478 1 0.82692
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "\n--- Fixed effects table ---"
## # A tibble: 12 × 9
## effect term estimate std.error statistic df p.value conf.low conf.high
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed (Interc… 4.74 2.34 2.03 93.4 4.54e-2 0.0986 9.38
## 2 fixed treatme… -1.30 0.579 -2.25 89.9 2.67e-2 -2.45 -0.154
## 3 fixed wai_c 0.0343 0.0375 0.916 92.0 3.62e-1 -0.0401 0.109
## 4 fixed race_wh… -0.0352 0.541 -0.0651 88.4 9.48e-1 -1.11 1.04
## 5 fixed meanpai… -0.439 0.0953 -4.61 94.1 1.26e-5 -0.629 -0.250
## 6 fixed age -0.0181 0.0311 -0.580 92.7 5.63e-1 -0.0799 0.0438
## 7 fixed female -0.970 0.585 -1.66 91.8 1.01e-1 -2.13 0.192
## 8 fixed attend_… -1.60 0.700 -2.29 93.7 2.43e-2 -2.99 -0.212
## 9 fixed treatme… -0.0440 0.0438 -1.00 91.4 3.18e-1 -0.131 0.0430
## 10 fixed treatme… -0.597 0.710 -0.841 92.1 4.02e-1 -2.01 0.812
## 11 fixed wai_c:r… -0.0225 0.0453 -0.497 89.8 6.20e-1 -0.112 0.0674
## 12 fixed treatme… -0.0116 0.0531 -0.219 89.0 8.27e-1 -0.117 0.0940
## NOTE: Results may be misleading due to involvement in interactions
##
## --- Estimated Marginal Means ---
## race_white_vs_bipoc = 0:
## treatment emmean SE df lower.CL upper.CL
## CBT -0.8679795 0.5604971 84.19 -1.982553 0.2465944
## EAET -2.1716139 0.5186525 78.50 -3.204067 -1.1391607
##
## race_white_vs_bipoc = 1:
## treatment emmean SE df lower.CL upper.CL
## CBT -0.9031703 0.5196198 71.20 -1.939213 0.1328728
## EAET -2.8037638 0.3992145 52.46 -3.604678 -2.0028497
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## --- Pairwise contrasts ---
## race_white_vs_bipoc = 0:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 1.303634 0.6221544 101.27 2.095 0.0386
##
## race_white_vs_bipoc = 1:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 1.900594 0.4611582 108.26 4.121 0.0001
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
##
## --- Interaction plot ---
##
## --- Simple slopes for race =0---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =0:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
##
## --- Simple slopes for race =1---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =1:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
##
## Fitting model for painint_change_post (main)
## painint_change_post ~ treatment * wai_c * race_white_vs_bipoc +
## painint_b + age + female + attend_prop + (1 | group)
## <environment: 0x116390828>
##
## --- Summary ---
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: fmla
## Data: data
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 636.6 672.3 -304.3 608.6 81
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.73839 -0.44770 0.03904 0.66257 2.01980
##
## Random effects:
## Groups Name Variance Std.Dev.
## group (Intercept) 3.095 1.759
## Residual 33.206 5.763
## Number of obs: 95, groups: group, 10
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 12.60971 9.73897 93.27345 1.295
## treatmentEAET -4.37889 2.23408 89.87175 -1.960
## wai_c 0.07440 0.14432 92.28792 0.516
## race_white_vs_bipoc -0.91524 2.09207 88.52966 -0.437
## painint_b -0.32852 0.08794 91.84350 -3.736
## age -0.02336 0.12244 93.18699 -0.191
## female -4.53516 2.25546 91.96304 -2.011
## attend_prop -4.56735 2.71029 94.23625 -1.685
## treatmentEAET:wai_c -0.15039 0.17190 91.76546 -0.875
## treatmentEAET:race_white_vs_bipoc 2.40300 2.72870 92.59023 0.881
## wai_c:race_white_vs_bipoc -0.02705 0.17689 90.12834 -0.153
## treatmentEAET:wai_c:race_white_vs_bipoc -0.10430 0.20839 89.32089 -0.500
## Pr(>|t|)
## (Intercept) 0.198596
## treatmentEAET 0.053089 .
## wai_c 0.607417
## race_white_vs_bipoc 0.662827
## painint_b 0.000325 ***
## age 0.849093
## female 0.047281 *
## attend_prop 0.095261 .
## treatmentEAET:wai_c 0.383945
## treatmentEAET:race_white_vs_bipoc 0.380793
## wai_c:race_white_vs_bipoc 0.878799
## treatmentEAET:wai_c:race_white_vs_bipoc 0.617956
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trEAET wai_c rc_w__ pnnt_b age female attnd_ trEAET:_
## tretmntEAET -0.053
## wai_c 0.028 -0.083
## rc_wht_vs_b -0.238 0.603 -0.131
## painint_b -0.448 0.125 0.134 -0.033
## age -0.916 -0.159 -0.003 0.117 0.191
## female -0.100 0.076 -0.258 0.220 -0.191 0.119
## attend_prop -0.048 0.126 -0.160 -0.035 0.077 -0.213 0.006
## trtmnEAET:_ -0.021 0.002 -0.799 0.121 -0.231 0.105 0.210 -0.161
## trtEAET:___ 0.025 -0.803 0.108 -0.756 -0.019 0.107 -0.219 -0.022 -0.053
## w_c:rc_wh__ 0.006 0.056 -0.816 0.130 -0.168 -0.006 0.218 0.102 0.667
## tEAET:_:___ 0.077 -0.001 0.663 -0.118 0.165 -0.127 -0.166 0.056 -0.806
## tEAET:__ w_:___
## tretmntEAET
## wai_c
## rc_wht_vs_b
## painint_b
## age
## female
## attend_prop
## trtmnEAET:_
## trtEAET:___
## w_c:rc_wh__ -0.104
## tEAET:_:___ 0.048 -0.831
##
## --- Type III ANOVA ---
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: painint_change_post
## Chisq Df Pr(>Chisq)
## (Intercept) 1.6764 1 0.1954001
## treatment 3.8418 1 0.0499912 *
## wai_c 0.2658 1 0.6061837
## race_white_vs_bipoc 0.1914 1 0.6617616
## painint_b 13.9569 1 0.0001871 ***
## age 0.0364 1 0.8486782
## female 4.0431 1 0.0443521 *
## attend_prop 2.8399 1 0.0919523 .
## treatment:wai_c 0.7653 1 0.3816606
## treatment:race_white_vs_bipoc 0.7755 1 0.3785129
## wai_c:race_white_vs_bipoc 0.0234 1 0.8784573
## treatment:wai_c:race_white_vs_bipoc 0.2505 1 0.6167252
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "\n--- Fixed effects table ---"
## # A tibble: 12 × 9
## effect term estimate std.error statistic df p.value conf.low conf.high
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed (Interc… 12.6 9.74 1.29 93.3 1.99e-1 -6.73 31.9
## 2 fixed treatme… -4.38 2.23 -1.96 89.9 5.31e-2 -8.82 0.0596
## 3 fixed wai_c 0.0744 0.144 0.516 92.3 6.07e-1 -0.212 0.361
## 4 fixed race_wh… -0.915 2.09 -0.437 88.5 6.63e-1 -5.07 3.24
## 5 fixed painint… -0.329 0.0879 -3.74 91.8 3.25e-4 -0.503 -0.154
## 6 fixed age -0.0234 0.122 -0.191 93.2 8.49e-1 -0.266 0.220
## 7 fixed female -4.54 2.26 -2.01 92.0 4.73e-2 -9.01 -0.0556
## 8 fixed attend_… -4.57 2.71 -1.69 94.2 9.53e-2 -9.95 0.814
## 9 fixed treatme… -0.150 0.172 -0.875 91.8 3.84e-1 -0.492 0.191
## 10 fixed treatme… 2.40 2.73 0.881 92.6 3.81e-1 -3.02 7.82
## 11 fixed wai_c:r… -0.0271 0.177 -0.153 90.1 8.79e-1 -0.378 0.324
## 12 fixed treatme… -0.104 0.208 -0.500 89.3 6.18e-1 -0.518 0.310
## NOTE: Results may be misleading due to involvement in interactions
##
## --- Estimated Marginal Means ---
## race_white_vs_bipoc = 0:
## treatment emmean SE df lower.CL upper.CL
## CBT -4.075853 2.139372 89.21 -8.326601 0.174895
## EAET -8.454744 1.955388 83.59 -12.343526 -4.565962
##
## race_white_vs_bipoc = 1:
## treatment emmean SE df lower.CL upper.CL
## CBT -4.991095 1.963626 76.36 -8.901697 -1.080494
## EAET -6.966988 1.467323 58.79 -9.903320 -4.030656
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## --- Pairwise contrasts ---
## race_white_vs_bipoc = 0:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 4.378891 2.401961 101.48 1.823 0.0712
##
## race_white_vs_bipoc = 1:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 1.975893 1.790221 106.92 1.104 0.2722
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
##
## --- Interaction plot ---
##
## --- Simple slopes for race =0---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =0:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
##
## --- Simple slopes for race =1---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =1:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
##
## Fitting model for depress_change_post (main)
## depress_change_post ~ treatment * wai_c * race_white_vs_bipoc +
## depress_b + age + female + attend_prop + (1 | group)
## <environment: 0x131f54eb8>
##
## --- Summary ---
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: fmla
## Data: data
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 585.2 621.0 -278.6 557.2 81
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.70362 -0.64934 -0.02602 0.66319 2.52348
##
## Random effects:
## Groups Name Variance Std.Dev.
## group (Intercept) 0.02305 0.1518
## Residual 20.62683 4.5417
## Number of obs: 95, groups: group, 10
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 17.95104 6.98393 94.70920 2.570
## treatmentEAET -5.84621 1.72916 91.96052 -3.381
## wai_c -0.07595 0.11120 94.65313 -0.683
## race_white_vs_bipoc -3.97799 1.66155 90.67096 -2.394
## depress_b -0.37738 0.06378 92.51564 -5.917
## age -0.11024 0.09303 94.96453 -1.185
## female -4.76223 1.73485 93.65143 -2.745
## attend_prop -0.74691 2.06600 94.07465 -0.362
## treatmentEAET:wai_c 0.07212 0.13101 94.78887 0.550
## treatmentEAET:race_white_vs_bipoc 3.51578 2.11774 94.17832 1.660
## wai_c:race_white_vs_bipoc 0.06533 0.13863 93.04541 0.471
## treatmentEAET:wai_c:race_white_vs_bipoc -0.22192 0.16272 93.42380 -1.364
## Pr(>|t|)
## (Intercept) 0.01172 *
## treatmentEAET 0.00106 **
## wai_c 0.49626
## race_white_vs_bipoc 0.01872 *
## depress_b 0.0000000548 ***
## age 0.23896
## female 0.00725 **
## attend_prop 0.71852
## treatmentEAET:wai_c 0.58328
## treatmentEAET:race_white_vs_bipoc 0.10021
## wai_c:race_white_vs_bipoc 0.63858
## treatmentEAET:wai_c:race_white_vs_bipoc 0.17590
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trEAET wai_c rc_w__ dprss_ age female attnd_ trEAET:_
## tretmntEAET 0.009
## wai_c 0.091 -0.099
## rc_wht_vs_b -0.226 0.607 -0.135
## depress_b -0.292 -0.030 0.113 -0.175
## age -0.936 -0.186 -0.047 0.112 0.119
## female -0.172 0.097 -0.243 0.230 -0.142 0.163
## attend_prop -0.034 0.112 -0.162 -0.048 0.056 -0.220 -0.024
## trtmnEAET:_ -0.095 0.034 -0.808 0.136 -0.157 0.139 0.202 -0.144
## trtEAET:___ 0.008 -0.819 0.104 -0.766 0.092 0.107 -0.220 0.002 -0.075
## w_c:rc_wh__ -0.035 0.082 -0.812 0.151 -0.204 0.023 0.211 0.107 0.667
## tEAET:_:___ 0.114 -0.026 0.670 -0.139 0.173 -0.142 -0.162 0.045 -0.807
## tEAET:__ w_:___
## tretmntEAET
## wai_c
## rc_wht_vs_b
## depress_b
## age
## female
## attend_prop
## trtmnEAET:_
## trtEAET:___
## w_c:rc_wh__ -0.113
## tEAET:_:___ 0.067 -0.836
##
## --- Type III ANOVA ---
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: depress_change_post
## Chisq Df Pr(>Chisq)
## (Intercept) 6.6066 1 0.0101600 *
## treatment 11.4309 1 0.0007223 ***
## wai_c 0.4665 1 0.4945892
## race_white_vs_bipoc 5.7319 1 0.0166593 *
## depress_b 35.0161 1 0.00000000327 ***
## age 1.4043 1 0.2360055
## female 7.5352 1 0.0060505 **
## attend_prop 0.1307 1 0.7177079
## treatment:wai_c 0.3030 1 0.5819848
## treatment:race_white_vs_bipoc 2.7561 1 0.0968837 .
## wai_c:race_white_vs_bipoc 0.2221 1 0.6374750
## treatment:wai_c:race_white_vs_bipoc 1.8600 1 0.1726251
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "\n--- Fixed effects table ---"
## # A tibble: 12 × 9
## effect term estimate std.error statistic df p.value conf.low conf.high
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed (Interc… 18.0 6.98 2.57 94.7 1.17e-2 4.09 31.8
## 2 fixed treatme… -5.85 1.73 -3.38 92.0 1.06e-3 -9.28 -2.41
## 3 fixed wai_c -0.0760 0.111 -0.683 94.7 4.96e-1 -0.297 0.145
## 4 fixed race_wh… -3.98 1.66 -2.39 90.7 1.87e-2 -7.28 -0.677
## 5 fixed depress… -0.377 0.0638 -5.92 92.5 5.48e-8 -0.504 -0.251
## 6 fixed age -0.110 0.0930 -1.19 95.0 2.39e-1 -0.295 0.0744
## 7 fixed female -4.76 1.73 -2.75 93.7 7.25e-3 -8.21 -1.32
## 8 fixed attend_… -0.747 2.07 -0.362 94.1 7.19e-1 -4.85 3.36
## 9 fixed treatme… 0.0721 0.131 0.550 94.8 5.83e-1 -0.188 0.332
## 10 fixed treatme… 3.52 2.12 1.66 94.2 1.00e-1 -0.689 7.72
## 11 fixed wai_c:r… 0.0653 0.139 0.471 93.0 6.39e-1 -0.210 0.341
## 12 fixed treatme… -0.222 0.163 -1.36 93.4 1.76e-1 -0.545 0.101
## NOTE: Results may be misleading due to involvement in interactions
##
## --- Estimated Marginal Means ---
## race_white_vs_bipoc = 0:
## treatment emmean SE df lower.CL upper.CL
## CBT -0.782627 1.565636 99.99 -3.888807 2.323552
## EAET -6.628839 1.437743 96.83 -9.482424 -3.775253
##
## race_white_vs_bipoc = 1:
## treatment emmean SE df lower.CL upper.CL
## CBT -4.760618 1.422323 89.37 -7.586584 -1.934651
## EAET -7.091050 1.053660 83.29 -9.186631 -4.995470
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## --- Pairwise contrasts ---
## race_white_vs_bipoc = 0:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 5.846211 1.865126 103.87 3.134 0.0022
##
## race_white_vs_bipoc = 1:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 2.330432 1.344195 102.69 1.734 0.0860
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
##
## --- Interaction plot ---
## boundary (singular) fit: see help('isSingular')
##
## --- Simple slopes for race =0---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =0:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
## boundary (singular) fit: see help('isSingular')
##
## --- Simple slopes for race =1---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =1:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
# -----------------------
# 4. Run sensitivity models
# -----------------------
att_results <- map2(names(outcomes), outcomes, ~ fit_and_probe(.x, .y, df_post, include_attendance = TRUE))
##
## Fitting model for pain_change_post (with attend)
## pain_change_post ~ treatment * wai_c * race_white_vs_bipoc +
## meanpain_b + age + female + attend_prop + attend + (1 | group)
## <environment: 0x1147df7e0>
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Summary ---
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: fmla
## Data: data
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 380.7 416.4 -176.3 352.7 81
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.21004 -0.61173 0.03432 0.62570 2.28741
##
## Random effects:
## Groups Name Variance Std.Dev.
## group (Intercept) 0.302 0.5496
## Residual 2.201 1.4834
## Number of obs: 95, groups: group, 10
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 4.73822 2.33651 93.37599 2.028
## treatmentEAET -1.30363 0.57870 89.86233 -2.253
## wai_c 0.03432 0.03748 91.99913 0.916
## race_white_vs_bipoc -0.03519 0.54054 88.35007 -0.065
## meanpain_b -0.43944 0.09531 94.07770 -4.611
## age -0.01805 0.03113 92.71804 -0.580
## female -0.96970 0.58504 91.83911 -1.658
## attend_prop -1.60329 0.70048 93.65939 -2.289
## treatmentEAET:wai_c -0.04395 0.04375 91.38556 -1.005
## treatmentEAET:race_white_vs_bipoc -0.59696 0.70963 92.05643 -0.841
## wai_c:race_white_vs_bipoc -0.02250 0.04527 89.80447 -0.497
## treatmentEAET:wai_c:race_white_vs_bipoc -0.01162 0.05313 89.04215 -0.219
## Pr(>|t|)
## (Intercept) 0.0454 *
## treatmentEAET 0.0267 *
## wai_c 0.3623
## race_white_vs_bipoc 0.9482
## meanpain_b 0.0000126 ***
## age 0.5634
## female 0.1008
## attend_prop 0.0243 *
## treatmentEAET:wai_c 0.3178
## treatmentEAET:race_white_vs_bipoc 0.4024
## wai_c:race_white_vs_bipoc 0.6204
## treatmentEAET:wai_c:race_white_vs_bipoc 0.8274
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trEAET wai_c rc_w__ mnpn_b age female attnd_ trEAET:_
## tretmntEAET -0.036
## wai_c 0.045 -0.074
## rc_wht_vs_b -0.251 0.592 -0.139
## meanpain_b -0.260 0.152 0.164 -0.075
## age -0.920 -0.182 -0.017 0.121 0.027
## female -0.142 0.069 -0.266 0.230 -0.208 0.148
## attend_prop -0.015 0.118 -0.171 -0.031 0.001 -0.233 0.029
## trtmnEAET:_ -0.095 0.008 -0.800 0.128 -0.156 0.148 0.198 -0.147
## trtEAET:___ 0.040 -0.806 0.095 -0.742 -0.105 0.113 -0.203 -0.025 -0.041
## w_c:rc_wh__ -0.042 0.061 -0.815 0.135 -0.108 0.019 0.210 0.116 0.659
## tEAET:_:___ 0.147 -0.011 0.655 -0.119 0.067 -0.160 -0.149 0.044 -0.798
## tEAET:__ w_:___
## tretmntEAET
## wai_c
## rc_wht_vs_b
## meanpain_b
## age
## female
## attend_prop
## trtmnEAET:_
## trtEAET:___
## w_c:rc_wh__ -0.099
## tEAET:_:___ 0.044 -0.826
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Type III ANOVA ---
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: pain_change_post
## Chisq Df Pr(>Chisq)
## (Intercept) 4.1124 1 0.04257 *
## treatment 5.0747 1 0.02428 *
## wai_c 0.8383 1 0.35988
## race_white_vs_bipoc 0.0042 1 0.94809
## meanpain_b 21.2574 1 0.000004016 ***
## age 0.3363 1 0.56199
## female 2.7473 1 0.09742 .
## attend_prop 5.2388 1 0.02209 *
## attend 0
## treatment:wai_c 1.0090 1 0.31513
## treatment:race_white_vs_bipoc 0.7077 1 0.40022
## wai_c:race_white_vs_bipoc 0.2470 1 0.61918
## treatment:wai_c:race_white_vs_bipoc 0.0478 1 0.82692
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "\n--- Fixed effects table ---"
## # A tibble: 12 × 9
## effect term estimate std.error statistic df p.value conf.low conf.high
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed (Interc… 4.74 2.34 2.03 93.4 4.54e-2 0.0986 9.38
## 2 fixed treatme… -1.30 0.579 -2.25 89.9 2.67e-2 -2.45 -0.154
## 3 fixed wai_c 0.0343 0.0375 0.916 92.0 3.62e-1 -0.0401 0.109
## 4 fixed race_wh… -0.0352 0.541 -0.0651 88.4 9.48e-1 -1.11 1.04
## 5 fixed meanpai… -0.439 0.0953 -4.61 94.1 1.26e-5 -0.629 -0.250
## 6 fixed age -0.0181 0.0311 -0.580 92.7 5.63e-1 -0.0799 0.0438
## 7 fixed female -0.970 0.585 -1.66 91.8 1.01e-1 -2.13 0.192
## 8 fixed attend_… -1.60 0.700 -2.29 93.7 2.43e-2 -2.99 -0.212
## 9 fixed treatme… -0.0440 0.0438 -1.00 91.4 3.18e-1 -0.131 0.0430
## 10 fixed treatme… -0.597 0.710 -0.841 92.1 4.02e-1 -2.01 0.812
## 11 fixed wai_c:r… -0.0225 0.0453 -0.497 89.8 6.20e-1 -0.112 0.0674
## 12 fixed treatme… -0.0116 0.0531 -0.219 89.0 8.27e-1 -0.117 0.0940
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## NOTE: Results may be misleading due to involvement in interactions
##
## --- Estimated Marginal Means ---
## race_white_vs_bipoc = 0:
## treatment emmean SE df lower.CL upper.CL
## CBT -0.8679795 0.5604971 84.19 -1.982553 0.2465944
## EAET -2.1716139 0.5186525 78.50 -3.204067 -1.1391607
##
## race_white_vs_bipoc = 1:
## treatment emmean SE df lower.CL upper.CL
## CBT -0.9031703 0.5196198 71.20 -1.939213 0.1328728
## EAET -2.8037638 0.3992145 52.46 -3.604678 -2.0028497
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## --- Pairwise contrasts ---
## race_white_vs_bipoc = 0:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 1.303634 0.6221544 101.27 2.095 0.0386
##
## race_white_vs_bipoc = 1:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 1.900594 0.4611582 108.26 4.121 0.0001
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
##
## --- Interaction plot ---
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Simple slopes for race =0---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =0:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Simple slopes for race =1---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =1:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
##
## Fitting model for painint_change_post (with attend)
## painint_change_post ~ treatment * wai_c * race_white_vs_bipoc +
## painint_b + age + female + attend_prop + attend + (1 | group)
## <environment: 0x132399dd0>
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Summary ---
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: fmla
## Data: data
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 636.6 672.3 -304.3 608.6 81
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.73839 -0.44770 0.03904 0.66257 2.01980
##
## Random effects:
## Groups Name Variance Std.Dev.
## group (Intercept) 3.095 1.759
## Residual 33.206 5.763
## Number of obs: 95, groups: group, 10
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 12.60971 9.73897 93.27345 1.295
## treatmentEAET -4.37889 2.23408 89.87175 -1.960
## wai_c 0.07440 0.14432 92.28792 0.516
## race_white_vs_bipoc -0.91524 2.09207 88.52966 -0.437
## painint_b -0.32852 0.08794 91.84350 -3.736
## age -0.02336 0.12244 93.18699 -0.191
## female -4.53516 2.25546 91.96304 -2.011
## attend_prop -4.56735 2.71029 94.23625 -1.685
## treatmentEAET:wai_c -0.15039 0.17190 91.76546 -0.875
## treatmentEAET:race_white_vs_bipoc 2.40300 2.72870 92.59023 0.881
## wai_c:race_white_vs_bipoc -0.02705 0.17689 90.12834 -0.153
## treatmentEAET:wai_c:race_white_vs_bipoc -0.10430 0.20839 89.32089 -0.500
## Pr(>|t|)
## (Intercept) 0.198596
## treatmentEAET 0.053089 .
## wai_c 0.607417
## race_white_vs_bipoc 0.662827
## painint_b 0.000325 ***
## age 0.849093
## female 0.047281 *
## attend_prop 0.095261 .
## treatmentEAET:wai_c 0.383945
## treatmentEAET:race_white_vs_bipoc 0.380793
## wai_c:race_white_vs_bipoc 0.878799
## treatmentEAET:wai_c:race_white_vs_bipoc 0.617956
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trEAET wai_c rc_w__ pnnt_b age female attnd_ trEAET:_
## tretmntEAET -0.053
## wai_c 0.028 -0.083
## rc_wht_vs_b -0.238 0.603 -0.131
## painint_b -0.448 0.125 0.134 -0.033
## age -0.916 -0.159 -0.003 0.117 0.191
## female -0.100 0.076 -0.258 0.220 -0.191 0.119
## attend_prop -0.048 0.126 -0.160 -0.035 0.077 -0.213 0.006
## trtmnEAET:_ -0.021 0.002 -0.799 0.121 -0.231 0.105 0.210 -0.161
## trtEAET:___ 0.025 -0.803 0.108 -0.756 -0.019 0.107 -0.219 -0.022 -0.053
## w_c:rc_wh__ 0.006 0.056 -0.816 0.130 -0.168 -0.006 0.218 0.102 0.667
## tEAET:_:___ 0.077 -0.001 0.663 -0.118 0.165 -0.127 -0.166 0.056 -0.806
## tEAET:__ w_:___
## tretmntEAET
## wai_c
## rc_wht_vs_b
## painint_b
## age
## female
## attend_prop
## trtmnEAET:_
## trtEAET:___
## w_c:rc_wh__ -0.104
## tEAET:_:___ 0.048 -0.831
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Type III ANOVA ---
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: painint_change_post
## Chisq Df Pr(>Chisq)
## (Intercept) 1.6764 1 0.1954001
## treatment 3.8418 1 0.0499912 *
## wai_c 0.2658 1 0.6061837
## race_white_vs_bipoc 0.1914 1 0.6617616
## painint_b 13.9569 1 0.0001871 ***
## age 0.0364 1 0.8486782
## female 4.0431 1 0.0443521 *
## attend_prop 2.8399 1 0.0919523 .
## attend 0
## treatment:wai_c 0.7653 1 0.3816606
## treatment:race_white_vs_bipoc 0.7755 1 0.3785129
## wai_c:race_white_vs_bipoc 0.0234 1 0.8784573
## treatment:wai_c:race_white_vs_bipoc 0.2505 1 0.6167252
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "\n--- Fixed effects table ---"
## # A tibble: 12 × 9
## effect term estimate std.error statistic df p.value conf.low conf.high
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed (Interc… 12.6 9.74 1.29 93.3 1.99e-1 -6.73 31.9
## 2 fixed treatme… -4.38 2.23 -1.96 89.9 5.31e-2 -8.82 0.0596
## 3 fixed wai_c 0.0744 0.144 0.516 92.3 6.07e-1 -0.212 0.361
## 4 fixed race_wh… -0.915 2.09 -0.437 88.5 6.63e-1 -5.07 3.24
## 5 fixed painint… -0.329 0.0879 -3.74 91.8 3.25e-4 -0.503 -0.154
## 6 fixed age -0.0234 0.122 -0.191 93.2 8.49e-1 -0.266 0.220
## 7 fixed female -4.54 2.26 -2.01 92.0 4.73e-2 -9.01 -0.0556
## 8 fixed attend_… -4.57 2.71 -1.69 94.2 9.53e-2 -9.95 0.814
## 9 fixed treatme… -0.150 0.172 -0.875 91.8 3.84e-1 -0.492 0.191
## 10 fixed treatme… 2.40 2.73 0.881 92.6 3.81e-1 -3.02 7.82
## 11 fixed wai_c:r… -0.0271 0.177 -0.153 90.1 8.79e-1 -0.378 0.324
## 12 fixed treatme… -0.104 0.208 -0.500 89.3 6.18e-1 -0.518 0.310
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## NOTE: Results may be misleading due to involvement in interactions
##
## --- Estimated Marginal Means ---
## race_white_vs_bipoc = 0:
## treatment emmean SE df lower.CL upper.CL
## CBT -4.075853 2.139372 89.21 -8.326601 0.174895
## EAET -8.454744 1.955388 83.59 -12.343526 -4.565962
##
## race_white_vs_bipoc = 1:
## treatment emmean SE df lower.CL upper.CL
## CBT -4.991095 1.963626 76.36 -8.901697 -1.080494
## EAET -6.966988 1.467323 58.79 -9.903320 -4.030656
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## --- Pairwise contrasts ---
## race_white_vs_bipoc = 0:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 4.378891 2.401961 101.48 1.823 0.0712
##
## race_white_vs_bipoc = 1:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 1.975893 1.790221 106.92 1.104 0.2722
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
##
## --- Interaction plot ---
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Simple slopes for race =0---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =0:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Simple slopes for race =1---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =1:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
##
## Fitting model for depress_change_post (with attend)
## depress_change_post ~ treatment * wai_c * race_white_vs_bipoc +
## depress_b + age + female + attend_prop + attend + (1 | group)
## <environment: 0x115fc0d00>
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Summary ---
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: fmla
## Data: data
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 585.2 621.0 -278.6 557.2 81
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.70362 -0.64934 -0.02602 0.66319 2.52348
##
## Random effects:
## Groups Name Variance Std.Dev.
## group (Intercept) 0.02305 0.1518
## Residual 20.62683 4.5417
## Number of obs: 95, groups: group, 10
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 17.95104 6.98393 94.70920 2.570
## treatmentEAET -5.84621 1.72916 91.96052 -3.381
## wai_c -0.07595 0.11120 94.65313 -0.683
## race_white_vs_bipoc -3.97799 1.66155 90.67096 -2.394
## depress_b -0.37738 0.06378 92.51564 -5.917
## age -0.11024 0.09303 94.96453 -1.185
## female -4.76223 1.73485 93.65143 -2.745
## attend_prop -0.74691 2.06600 94.07465 -0.362
## treatmentEAET:wai_c 0.07212 0.13101 94.78887 0.550
## treatmentEAET:race_white_vs_bipoc 3.51578 2.11774 94.17832 1.660
## wai_c:race_white_vs_bipoc 0.06533 0.13863 93.04541 0.471
## treatmentEAET:wai_c:race_white_vs_bipoc -0.22192 0.16272 93.42380 -1.364
## Pr(>|t|)
## (Intercept) 0.01172 *
## treatmentEAET 0.00106 **
## wai_c 0.49626
## race_white_vs_bipoc 0.01872 *
## depress_b 0.0000000548 ***
## age 0.23896
## female 0.00725 **
## attend_prop 0.71852
## treatmentEAET:wai_c 0.58328
## treatmentEAET:race_white_vs_bipoc 0.10021
## wai_c:race_white_vs_bipoc 0.63858
## treatmentEAET:wai_c:race_white_vs_bipoc 0.17590
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trEAET wai_c rc_w__ dprss_ age female attnd_ trEAET:_
## tretmntEAET 0.009
## wai_c 0.091 -0.099
## rc_wht_vs_b -0.226 0.607 -0.135
## depress_b -0.292 -0.030 0.113 -0.175
## age -0.936 -0.186 -0.047 0.112 0.119
## female -0.172 0.097 -0.243 0.230 -0.142 0.163
## attend_prop -0.034 0.112 -0.162 -0.048 0.056 -0.220 -0.024
## trtmnEAET:_ -0.095 0.034 -0.808 0.136 -0.157 0.139 0.202 -0.144
## trtEAET:___ 0.008 -0.819 0.104 -0.766 0.092 0.107 -0.220 0.002 -0.075
## w_c:rc_wh__ -0.035 0.082 -0.812 0.151 -0.204 0.023 0.211 0.107 0.667
## tEAET:_:___ 0.114 -0.026 0.670 -0.139 0.173 -0.142 -0.162 0.045 -0.807
## tEAET:__ w_:___
## tretmntEAET
## wai_c
## rc_wht_vs_b
## depress_b
## age
## female
## attend_prop
## trtmnEAET:_
## trtEAET:___
## w_c:rc_wh__ -0.113
## tEAET:_:___ 0.067 -0.836
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Type III ANOVA ---
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: depress_change_post
## Chisq Df Pr(>Chisq)
## (Intercept) 6.6066 1 0.0101600 *
## treatment 11.4309 1 0.0007223 ***
## wai_c 0.4665 1 0.4945892
## race_white_vs_bipoc 5.7319 1 0.0166593 *
## depress_b 35.0161 1 0.00000000327 ***
## age 1.4043 1 0.2360055
## female 7.5352 1 0.0060505 **
## attend_prop 0.1307 1 0.7177079
## attend 0
## treatment:wai_c 0.3030 1 0.5819848
## treatment:race_white_vs_bipoc 2.7561 1 0.0968837 .
## wai_c:race_white_vs_bipoc 0.2221 1 0.6374750
## treatment:wai_c:race_white_vs_bipoc 1.8600 1 0.1726251
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "\n--- Fixed effects table ---"
## # A tibble: 12 × 9
## effect term estimate std.error statistic df p.value conf.low conf.high
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed (Interc… 18.0 6.98 2.57 94.7 1.17e-2 4.09 31.8
## 2 fixed treatme… -5.85 1.73 -3.38 92.0 1.06e-3 -9.28 -2.41
## 3 fixed wai_c -0.0760 0.111 -0.683 94.7 4.96e-1 -0.297 0.145
## 4 fixed race_wh… -3.98 1.66 -2.39 90.7 1.87e-2 -7.28 -0.677
## 5 fixed depress… -0.377 0.0638 -5.92 92.5 5.48e-8 -0.504 -0.251
## 6 fixed age -0.110 0.0930 -1.19 95.0 2.39e-1 -0.295 0.0744
## 7 fixed female -4.76 1.73 -2.75 93.7 7.25e-3 -8.21 -1.32
## 8 fixed attend_… -0.747 2.07 -0.362 94.1 7.19e-1 -4.85 3.36
## 9 fixed treatme… 0.0721 0.131 0.550 94.8 5.83e-1 -0.188 0.332
## 10 fixed treatme… 3.52 2.12 1.66 94.2 1.00e-1 -0.689 7.72
## 11 fixed wai_c:r… 0.0653 0.139 0.471 93.0 6.39e-1 -0.210 0.341
## 12 fixed treatme… -0.222 0.163 -1.36 93.4 1.76e-1 -0.545 0.101
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## NOTE: Results may be misleading due to involvement in interactions
##
## --- Estimated Marginal Means ---
## race_white_vs_bipoc = 0:
## treatment emmean SE df lower.CL upper.CL
## CBT -0.782627 1.565636 99.99 -3.888807 2.323552
## EAET -6.628839 1.437743 96.83 -9.482424 -3.775253
##
## race_white_vs_bipoc = 1:
## treatment emmean SE df lower.CL upper.CL
## CBT -4.760618 1.422323 89.37 -7.586584 -1.934651
## EAET -7.091050 1.053660 83.29 -9.186631 -4.995470
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## --- Pairwise contrasts ---
## race_white_vs_bipoc = 0:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 5.846211 1.865126 103.87 3.134 0.0022
##
## race_white_vs_bipoc = 1:
## contrast estimate SE df t.ratio p.value
## CBT - EAET 2.330432 1.344195 102.69 1.734 0.0860
##
## Results are averaged over the levels of: wai_c, female
## Degrees-of-freedom method: kenward-roger
##
## --- Interaction plot ---
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## boundary (singular) fit: see help('isSingular')
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Simple slopes for race =0---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =0:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## boundary (singular) fit: see help('isSingular')
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
##
## --- Simple slopes for race =1---
## Warning: Johnson-Neyman intervals are not available for factor moderators.
## sim_slopes failed for race =1:error in evaluating the argument 'x' in selecting a method for function 'print': object 'dv' not found