Cohen's d

# Cohen's d: CBT vs. Control for Depression
library(effsize)     # cohen.d()
library(ggplot2)
library(dplyr)

# Simulate data (or load: data <- read.csv("depression_rct.csv"))
set.seed(2025)
data <- data.frame(
  group = rep(c("CBT", "Control"), each=50),
  depression = c(
    rnorm(50, mean=12.3, sd=8.2),   # CBT: M=12.3, SD=8.2
    rnorm(50, mean=23.7, sd=9.5)    # Control: M=23.7, SD=9.5
  )
)

# === STEP 1: Descriptive Statistics ===
data %>%
  group_by(group) %>%
  summarise(
    n = n(),
    M = mean(depression),
    SD = sd(depression),
    Min = min(depression),
    Max = max(depression)
  )

# === STEP 2: Visual Comparison ===
# Overlaid density plots
ggplot(data, aes(x=depression, fill=group)) +
  geom_density(alpha=0.5) +
  geom_vline(data = data %>% group_by(group) %>% 
               summarise(M=mean(depression)),
             aes(xintercept=M, color=group), linetype="dashed", size=1) +
  labs(title="Distribution of Depression Scores by Group",
       subtitle="Dashed lines = group means",
       x="BDI-II Score(0-63)", y="Density") +
  scale_fill_brewer(palette="Set1") +
  scale_color_brewer(palette="Set1") +
  theme_classic()

# === STEP 3: Calculate Cohen's d ===
# Method 1: Using effsize package
cohen_result <- cohen.d(depression ~ group, data=data)
print(cohen_result)
# Output: Cohen's d estimate: -1.29 (large)
# 95% CI: [-1.71, -0.87]

# Method 2: Manual calculation (for understanding)
M_cbt <- mean(data$depression[data$group == "CBT"])
M_control <- mean(data$depression[data$group == "Control"])
SD_cbt <- sd(data$depression[data$group == "CBT"])
SD_control <- sd(data$depression[data$group == "Control"])
n_cbt <- sum(data$group == "CBT")
n_control <- sum(data$group == "Control")

# Pooled standard deviation
SD_pooled <- sqrt(((n_cbt-1)*SD_cbt^2 + (n_control-1)*SD_control^2) / 
                    (n_cbt + n_control - 2))

cohen_d <- (M_cbt - M_control) / SD_pooled

cat("\n=== Manual Calculation ===")
cat("\nMean CBT:", round(M_cbt, 2))
cat("\nMean Control:", round(M_control, 2))
cat("\nMean Difference:", round(M_cbt - M_control, 2))
cat("\nPooled SD:", round(SD_pooled, 2))
cat("\nCohen's d:", round(cohen_d, 2), "\n")

# Interpretation
if (abs(cohen_d) < 0.2) {
  interpretation <- "negligible"
} else if (abs(cohen_d) < 0.5) {
  interpretation <- "small"
} else if (abs(cohen_d) < 0.8) {
  interpretation <- "medium"
} else {
  interpretation <- "large"
}
cat("Effect size:", interpretation, "(Cohen, 1988)\n")

# === STEP 4: Hedges' g (bias-corrected for small samples) ===
# Correction factor J
df <- n_cbt + n_control - 2
J <- 1 - (3 / (4*df - 1))
hedges_g <- cohen_d * J

cat("\nHedges' g(bias-corrected):", round(hedges_g, 2))
cat("\nCorrection factor J:", round(J, 4), "\n")

# === STEP 5: Alternative Effect Sizes ===
# Glass's delta (standardize by control SD only)
glass_delta <- (M_cbt - M_control) / SD_control
cat("\nGlass's Δ (control SD):", round(glass_delta, 2), "\n")

# Unstandardized difference with 95% CI
library(rstatix)
t_result <- t.test(depression ~ group, data=data)
cat("\nUnstandardized difference:", round(M_cbt - M_control, 2), 
    "points on BDI-II")
cat("\n95% CI:", round(t_result$conf.int, 2), "\n")

# === STEP 6: Visualization with Effect Size ===
data_summary <- data %>%
  group_by(group) %>%
  summarise(
    M = mean(depression),
    SD = sd(depression),
    SE = SD/sqrt(n())
  )

ggplot(data_summary, aes(x=group, y=M, fill=group)) +
  geom_bar(stat="identity", width=0.6, alpha=0.8) +
  geom_errorbar(aes(ymin=M-1.96*SE, ymax=M+1.96*SE), width=0.2) +
  annotate("text", x=1.5, y=max(data_summary$M)-5, 
           label=paste0("Cohen's d = ", round(abs(cohen_d), 2), " (large)"),
           size=5, fontface="bold") +
  labs(title="CBT Effect on Depression",
       subtitle=paste0("Mean difference = ", round(M_control - M_cbt, 1), 
                       " points on BDI-II"),
       x="Group", y="Mean BDI-II Score ± 95% CI") +
  scale_fill_brewer(palette="Set1") +
  theme_classic() +
  theme(legend.position="none")

# === APA-Style Reporting ===
cat("\n=== APA-Style Report ===\n")
cat("Cohen's d was calculated to quantify the magnitude of CBT effect on\n")
cat("depression. The CBT group(M =", round(M_cbt, 1), ", SD =", round(SD_cbt, 1),
    ") showed\n")
cat("substantially lower depression than the control group(M =", round(M_control, 1),
    ", SD =", round(SD_control, 1), "),\n")
cat("d =", round(cohen_d, 2), "(95% CI [", round(cohen_result$conf.int[1], 2),
    ",", round(cohen_result$conf.int[2], 2), "]).\n")
cat("This represents a large effect(Cohen, 1988), indicating CBT produced\n")
cat("a clinically meaningful reduction in depression symptoms.\n")