# =============================================================================
# STEP 2: DESCRIPTIVE STATISTICS FOR MEDIATION VARIABLES
# =============================================================================

library(dplyr)

cat("\n=== CORRELATIONS FOR MEDIATION ===\n\n")
## 
## === CORRELATIONS FOR MEDIATION ===
cor_vars <- df_6m %>%
  dplyr::select(dplyr::any_of(c("eaet", "attend", "wai_c", "meanpain_6m", "painint_6m", "depress_6m")))

missing_vars <- setdiff(c("eaet", "attend", "wai_c", "meanpain_6m", "painint_6m", "depress_6m"),
                        names(cor_vars))
if (length(missing_vars) > 0) {
  warning("Missing variables from df_6m: ", paste(missing_vars, collapse = ", "))
}

if (ncol(cor_vars) > 1) {
  cor_vars <- cor(cor_vars, use = "pairwise.complete.obs") %>%
    round(3)
  print(cor_vars)
} else {
  stop("Not enough valid variables found to compute correlations.")
}
##               eaet attend  wai_c meanpain_6m painint_6m depress_6m
## eaet         1.000 -0.118  0.045      -0.330     -0.109     -0.144
## attend      -0.118  1.000  0.469      -0.033     -0.032     -0.126
## wai_c        0.045  0.469  1.000      -0.090      0.075     -0.071
## meanpain_6m -0.330 -0.033 -0.090       1.000      0.664      0.445
## painint_6m  -0.109 -0.032  0.075       0.664      1.000      0.488
## depress_6m  -0.144 -0.126 -0.071       0.445      0.488      1.000
# =============================================================================
# STEP 3: SIMPLE MEDIATION #1 - WAI AS MEDIATOR
# Treatment → WAI → 6-month Outcomes
# =============================================================================

cat("\n\n")
cat("=============================================================================\n")
## =============================================================================
cat("SIMPLE MEDIATION #1: WAI AS MEDIATOR\n")
## SIMPLE MEDIATION #1: WAI AS MEDIATOR
cat("Path: Treatment → WAI (post-treatment) → 6-month Outcomes\n")
## Path: Treatment → WAI (post-treatment) → 6-month Outcomes
cat("=============================================================================\n\n")
## =============================================================================
# -----------------------------------------------------------------------------
# Outcome 1: Pain Severity at 6 months
# -----------------------------------------------------------------------------

cat("--- OUTCOME: Pain Severity at 6 months ---\n\n")
## --- OUTCOME: Pain Severity at 6 months ---
# Step 1: Treatment → WAI (Path a)
cat("Step 1: Treatment → WAI (Path a)\n")
## Step 1: Treatment → WAI (Path a)
model_tx_wai_pain <- lm(wai ~ eaet + meanpain_b + age + female + race_white_vs_bipoc,
                        data = df_6m)
summary(model_tx_wai_pain)
## 
## Call:
## lm(formula = wai ~ eaet + meanpain_b + age + female + race_white_vs_bipoc, 
##     data = df_6m)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.328  -8.315   1.575  10.660  21.811 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         38.222651  20.315996   1.881   0.0629 .
## eaet                 0.731239   2.755252   0.265   0.7913  
## meanpain_b          -0.002985   0.796545  -0.004   0.9970  
## age                  0.318431   0.258928   1.230   0.2217  
## female               4.771969   4.905129   0.973   0.3330  
## race_white_vs_bipoc  1.916581   3.199584   0.599   0.5506  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.78 on 97 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.02424,    Adjusted R-squared:  -0.02606 
## F-statistic: 0.4819 on 5 and 97 DF,  p-value: 0.789
# Step 2: Treatment + WAI → Outcome (Paths b and c')
cat("\nStep 2: Treatment + WAI → Pain at 6m (Paths b and c')\n")
## 
## Step 2: Treatment + WAI → Pain at 6m (Paths b and c')
model_full_pain <- lm(meanpain_6m ~ eaet + wai + meanpain_b + age + female + race_white_vs_bipoc,
                      data = df_6m)
summary(model_full_pain)
## 
## Call:
## lm(formula = meanpain_6m ~ eaet + wai + meanpain_b + age + female + 
##     race_white_vs_bipoc, data = df_6m)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7796 -0.9936  0.2147  1.0934  3.9205 
## 
## Coefficients:
##                      Estimate Std. Error t value   Pr(>|t|)    
## (Intercept)          4.922512   2.561420   1.922     0.0576 .  
## eaet                -1.095775   0.341333  -3.210     0.0018 ** 
## wai                 -0.009828   0.012574  -0.782     0.4364    
## meanpain_b           0.545966   0.098644   5.535 0.00000027 ***
## age                 -0.023245   0.032315  -0.719     0.4737    
## female              -0.441859   0.610406  -0.724     0.4709    
## race_white_vs_bipoc -0.040899   0.396968  -0.103     0.9182    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.706 on 96 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.347,  Adjusted R-squared:  0.3062 
## F-statistic: 8.503 on 6 and 96 DF,  p-value: 0.0000002076
# Bootstrap mediation test
cat("\nBootstrap Mediation Analysis (5000 iterations):\n")
## 
## Bootstrap Mediation Analysis (5000 iterations):
med_wai_pain <- mediate(model_tx_wai_pain, model_full_pain,
                        treat = "eaet", mediator = "wai",
                        boot = TRUE, sims = 5000, boot.ci.type = "perc")
## Running nonparametric bootstrap
summary(med_wai_pain)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                  Estimate 95% CI Lower 95% CI Upper p-value   
## ACME           -0.0071863   -0.1302531    0.0764553  0.8612   
## ADE            -1.0957754   -1.7813494   -0.4443440  0.0016 **
## Total Effect   -1.1029618   -1.8196045   -0.4471005  0.0016 **
## Prop. Mediated  0.0065155   -0.0807171    0.1229927  0.8612   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 103 
## 
## 
## Simulations: 5000
# Save key results
wai_pain_results <- data.frame(
  Outcome = "Pain Severity (6m)",
  Mediator = "WAI",
  Path_a = coef(model_tx_wai_pain)["eaet"],
  Path_b = coef(model_full_pain)["wai"],
  Direct_Effect = med_wai_pain$d0,
  Indirect_Effect = med_wai_pain$z0,
  Total_Effect = med_wai_pain$tau.coef,
  Prop_Mediated = med_wai_pain$n0,
  ACME_p = med_wai_pain$z0.p,
  Direct_p = med_wai_pain$d0.p
)

# -----------------------------------------------------------------------------
# Outcome 2: Pain Interference at 6 months
# -----------------------------------------------------------------------------

cat("\n\n--- OUTCOME: Pain Interference at 6 months ---\n\n")
## 
## 
## --- OUTCOME: Pain Interference at 6 months ---
model_tx_wai_painint <- lm(wai ~ eaet + painint_b + age + female + race_white_vs_bipoc,
                           data = df_6m)
model_full_painint <- lm(painint_6m ~ eaet + wai + painint_b + age + female + race_white_vs_bipoc,
                         data = df_6m)

cat("Bootstrap Mediation Analysis:\n")
## Bootstrap Mediation Analysis:
med_wai_painint <- mediate(model_tx_wai_painint, model_full_painint,
                           treat = "eaet", mediator = "wai",
                           boot = TRUE, sims = 5000, boot.ci.type = "perc")
## Running nonparametric bootstrap
summary(med_wai_painint)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                 Estimate 95% CI Lower 95% CI Upper p-value
## ACME           -0.036847    -0.586868     0.301267  0.8192
## ADE            -0.074847    -2.524368     2.290552  0.9192
## Total Effect   -0.111694    -2.660433     2.278803  0.8796
## Prop. Mediated  0.329891    -1.330739     1.728477  0.8628
## 
## Sample Size Used: 103 
## 
## 
## Simulations: 5000
wai_painint_results <- data.frame(
  Outcome = "Pain Interference (6m)",
  Mediator = "WAI",
  Path_a = coef(model_tx_wai_painint)["eaet"],
  Path_b = coef(model_full_painint)["wai"],
  Direct_Effect = med_wai_painint$d0,
  Indirect_Effect = med_wai_painint$z0,
  Total_Effect = med_wai_painint$tau.coef,
  Prop_Mediated = med_wai_painint$n0,
  ACME_p = med_wai_painint$z0.p,
  Direct_p = med_wai_painint$d0.p
)

# -----------------------------------------------------------------------------
# Outcome 3: Depression at 6 months
# -----------------------------------------------------------------------------

cat("\n\n--- OUTCOME: Depression at 6 months ---\n\n")
## 
## 
## --- OUTCOME: Depression at 6 months ---
model_tx_wai_depress <- lm(wai ~ eaet + depress_b + age + female + race_white_vs_bipoc,
                           data = df_6m)
model_full_depress <- lm(depress_6m ~ eaet + wai + depress_b + age + female + race_white_vs_bipoc,
                         data = df_6m)

cat("Bootstrap Mediation Analysis:\n")
## Bootstrap Mediation Analysis:
med_wai_depress <- mediate(model_tx_wai_depress, model_full_depress,
                           treat = "eaet", mediator = "wai",
                           boot = TRUE, sims = 5000, boot.ci.type = "perc")
## Running nonparametric bootstrap
summary(med_wai_depress)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                 Estimate 95% CI Lower 95% CI Upper p-value  
## ACME           -0.063045    -0.661909     0.416600  0.7728  
## ADE            -2.045571    -3.988597    -0.012495  0.0480 *
## Total Effect   -2.108616    -4.128859    -0.016855  0.0492 *
## Prop. Mediated  0.029899    -0.399793     0.454050  0.7604  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 103 
## 
## 
## Simulations: 5000
wai_depress_results <- data.frame(
  Outcome = "Depression (6m)",
  Mediator = "WAI",
  Path_a = coef(model_tx_wai_depress)["eaet"],
  Path_b = coef(model_full_depress)["wai"],
  Direct_Effect = med_wai_depress$d0,
  Indirect_Effect = med_wai_depress$z0,
  Total_Effect = med_wai_depress$tau.coef,
  Prop_Mediated = med_wai_depress$n0,
  ACME_p = med_wai_depress$z0.p,
  Direct_p = med_wai_depress$d0.p
)

# Combine results
wai_mediation_summary <- bind_rows(wai_pain_results, wai_painint_results, wai_depress_results)
view(wai_mediation_summary)
# =============================================================================
# STEP 4: SIMPLE MEDIATION #2 - ATTENDANCE AS MEDIATOR
# Treatment → Attendance → 6-month Outcomes
# =============================================================================

cat("\n\n")
cat("=============================================================================\n")
## =============================================================================
cat("SIMPLE MEDIATION #2: ATTENDANCE AS MEDIATOR\n")
## SIMPLE MEDIATION #2: ATTENDANCE AS MEDIATOR
cat("Path: Treatment → Attendance → 6-month Outcomes\n")
## Path: Treatment → Attendance → 6-month Outcomes
cat("=============================================================================\n\n")
## =============================================================================
# -----------------------------------------------------------------------------
# Outcome 1: Pain Severity at 6 months
# -----------------------------------------------------------------------------

cat("--- OUTCOME: Pain Severity at 6 months ---\n\n")
## --- OUTCOME: Pain Severity at 6 months ---
# Step 1: Treatment → Attendance (Path a)
cat("Step 1: Treatment → Attendance (Path a)\n")
## Step 1: Treatment → Attendance (Path a)
model_tx_attend_pain <- lm(attend ~ eaet + meanpain_b + age + female + race_white_vs_bipoc,
                           data = df_6m)
summary(model_tx_attend_pain)
## 
## Call:
## lm(formula = attend ~ eaet + meanpain_b + age + female + race_white_vs_bipoc, 
##     data = df_6m)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1511 -1.0564  0.5295  1.6830  2.6105 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         -0.97059    3.40682  -0.285    0.776  
## eaet                -0.67806    0.47686  -1.422    0.159  
## meanpain_b           0.02461    0.13201   0.186    0.853  
## age                  0.09673    0.04381   2.208    0.030 *
## female               1.08832    0.83763   1.299    0.197  
## race_white_vs_bipoc  0.50147    0.53831   0.932    0.354  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.183 on 84 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared:  0.07752,    Adjusted R-squared:  0.02261 
## F-statistic: 1.412 on 5 and 84 DF,  p-value: 0.2284
# Step 2: Treatment + Attendance → Outcome (Paths b and c')
cat("\nStep 2: Treatment + Attendance → Pain at 6m (Paths b and c')\n")
## 
## Step 2: Treatment + Attendance → Pain at 6m (Paths b and c')
model_full_attend_pain <- lm(meanpain_6m ~ eaet + attend + meanpain_b + age + female + race_white_vs_bipoc,
                             data = df_6m)
summary(model_full_attend_pain)
## 
## Call:
## lm(formula = meanpain_6m ~ eaet + attend + meanpain_b + age + 
##     female + race_white_vs_bipoc, data = df_6m)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.745 -1.117  0.263  1.117  3.987 
## 
## Coefficients:
##                     Estimate Std. Error t value   Pr(>|t|)    
## (Intercept)          4.90700    2.80123   1.752    0.08351 .  
## eaet                -1.15470    0.39659  -2.912    0.00462 ** 
## attend              -0.06260    0.08967  -0.698    0.48704    
## meanpain_b           0.51404    0.10851   4.737 0.00000884 ***
## age                 -0.02430    0.03704  -0.656    0.51359    
## female              -0.34317    0.69528  -0.494    0.62291    
## race_white_vs_bipoc  0.06136    0.44469   0.138    0.89058    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.794 on 83 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared:  0.3152, Adjusted R-squared:  0.2657 
## F-statistic: 6.367 on 6 and 83 DF,  p-value: 0.00001526
# Bootstrap mediation test
cat("\nBootstrap Mediation Analysis (5000 iterations):\n")
## 
## Bootstrap Mediation Analysis (5000 iterations):
med_attend_pain <- mediate(model_tx_attend_pain, model_full_attend_pain,
                           treat = "eaet", mediator = "attend",
                           boot = TRUE, sims = 5000, boot.ci.type = "perc")
## Running nonparametric bootstrap
summary(med_attend_pain)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                 Estimate 95% CI Lower 95% CI Upper p-value   
## ACME            0.042449    -0.094530     0.235394  0.5664   
## ADE            -1.154702    -1.956904    -0.394502  0.0036 **
## Total Effect   -1.112253    -1.869287    -0.413990  0.0016 **
## Prop. Mediated -0.038165    -0.224431     0.132717  0.5672   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 90 
## 
## 
## Simulations: 5000
attend_pain_results <- data.frame(
  Outcome = "Pain Severity (6m)",
  Mediator = "Attendance",
  Path_a = coef(model_tx_attend_pain)["eaet"],
  Path_b = coef(model_full_attend_pain)["attend"],
  Direct_Effect = med_attend_pain$d0,
  Indirect_Effect = med_attend_pain$z0,
  Total_Effect = med_attend_pain$tau.coef,
  Prop_Mediated = med_attend_pain$n0,
  ACME_p = med_attend_pain$z0.p,
  Direct_p = med_attend_pain$d0.p
)

# -----------------------------------------------------------------------------
# Outcome 2: Pain Interference at 6 months
# -----------------------------------------------------------------------------

cat("\n\n--- OUTCOME: Pain Interference at 6 months ---\n\n")
## 
## 
## --- OUTCOME: Pain Interference at 6 months ---
model_tx_attend_painint <- lm(attend ~ eaet + painint_b + age + female + race_white_vs_bipoc,
                              data = df_6m)
model_full_attend_painint <- lm(painint_6m ~ eaet + attend + painint_b + age + female + race_white_vs_bipoc,
                                data = df_6m)

cat("Bootstrap Mediation Analysis:\n")
## Bootstrap Mediation Analysis:
med_attend_painint <- mediate(model_tx_attend_painint, model_full_attend_painint,
                              treat = "eaet", mediator = "attend",
                              boot = TRUE, sims = 5000, boot.ci.type = "perc")
## Running nonparametric bootstrap
## Warning in predict.lm(new.fit.M, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.M, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
summary(med_attend_painint)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME            0.17250     -0.27349      0.90238  0.4640
## ADE            -0.91761     -3.90879      1.78416  0.4936
## Total Effect   -0.74510     -3.57892      1.76307  0.5580
## Prop. Mediated -0.23152     -2.71786      1.96726  0.6420
## 
## Sample Size Used: 90 
## 
## 
## Simulations: 5000
attend_painint_results <- data.frame(
  Outcome = "Pain Interference (6m)",
  Mediator = "Attendance",
  Path_a = coef(model_tx_attend_painint)["eaet"],
  Path_b = coef(model_full_attend_painint)["attend"],
  Direct_Effect = med_attend_painint$d0,
  Indirect_Effect = med_attend_painint$z0,
  Total_Effect = med_attend_painint$tau.coef,
  Prop_Mediated = med_attend_painint$n0,
  ACME_p = med_attend_painint$z0.p,
  Direct_p = med_attend_painint$d0.p
)

# -----------------------------------------------------------------------------
# Outcome 3: Depression at 6 months
# -----------------------------------------------------------------------------

cat("\n\n--- OUTCOME: Depression at 6 months ---\n\n")
## 
## 
## --- OUTCOME: Depression at 6 months ---
model_tx_attend_depress <- lm(attend ~ eaet + depress_b + age + female + race_white_vs_bipoc,
                              data = df_6m)
model_full_attend_depress <- lm(depress_6m ~ eaet + attend + depress_b + age + female + race_white_vs_bipoc,
                                data = df_6m)

cat("Bootstrap Mediation Analysis:\n")
## Bootstrap Mediation Analysis:
med_attend_depress <- mediate(model_tx_attend_depress, model_full_attend_depress,
                              treat = "eaet", mediator = "attend",
                              boot = TRUE, sims = 5000, boot.ci.type = "perc")
## Running nonparametric bootstrap
## Warning in predict.lm(new.fit.M, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.M, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.M, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.M, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.t):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
## Warning in predict.lm(new.fit.Y, type = "response", newdata = pred.data.c):
## prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases
summary(med_attend_depress)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                 Estimate 95% CI Lower 95% CI Upper p-value   
## ACME            0.320239    -0.086266     0.971576  0.1324   
## ADE            -3.460175    -5.760767    -1.183761  0.0024 **
## Total Effect   -3.139936    -5.414516    -0.926277  0.0056 **
## Prop. Mediated -0.101989    -0.474113     0.031205  0.1380   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 90 
## 
## 
## Simulations: 5000
attend_depress_results <- data.frame(
  Outcome = "Depression (6m)",
  Mediator = "Attendance",
  Path_a = coef(model_tx_attend_depress)["eaet"],
  Path_b = coef(model_full_attend_depress)["attend"],
  Direct_Effect = med_attend_depress$d0,
  Indirect_Effect = med_attend_depress$z0,
  Total_Effect = med_attend_depress$tau.coef,
  Prop_Mediated = med_attend_depress$n0,
  ACME_p = med_attend_depress$z0.p,
  Direct_p = med_attend_depress$d0.p
)

# Combine results
attend_mediation_summary <- bind_rows(attend_pain_results, attend_painint_results, attend_depress_results)