Point-Biserial Correlation

# Point-Biserial Correlation: Meditation (binary) and Stress (continuous)
# Research: Are meditators less stressed than non-meditators?

# Install/load packages
library(tidyverse)
library(psych)       # For describe.by
library(car)         # For Levene's test
library(effsize)     # For Cohen's d

# Simulate realistic data (or load: data <- read.csv("meditation_stress.csv"))
set.seed(2025)
data <- data.frame(
  meditator = rep(c(0, 1), each=75),  # Binary: 0=no, 1=yes
  stress_score = c(
    rnorm(75, mean=26.5, sd=7.2),      # Non-meditators: higher stress
    rnorm(75, mean=22.1, sd=6.8)       # Meditators: lower stress
  )
)

# Convert to factor for clarity (optional)
data$meditator_factor <- factor(data$meditator, 
                                levels=c(0,1), 
                                labels=c("Non-meditator", "Meditator"))

# === STEP 1: Check Assumptions ===

# 1. Confirm binary variable has exactly 2 levels
cat("=== Binary Variable Check ===\n")
cat("Unique values:", unique(data$meditator), "\n")
table(data$meditator_factor)

# 2. Normality within each group (Shapiro-Wilk)
cat("\n=== Normality Tests(Shapiro-Wilk) ===\n")
by(data$stress_score, data$meditator_factor, shapiro.test)
# Look for p > .05 in both groups

# Q-Q plots by group
par(mfrow=c(1,2))
qqnorm(data$stress_score[data$meditator==0], main="Non-meditators Q-Q Plot")
qqline(data$stress_score[data$meditator==0])
qqnorm(data$stress_score[data$meditator==1], main="Meditators Q-Q Plot")
qqline(data$stress_score[data$meditator==1])

# 3. Homogeneity of variance (Levene's test)
cat("\n=== Levene's Test ===\n")
levene_result <- leveneTest(stress_score ~ meditator_factor, data=data)
print(levene_result)
# p > .05 indicates equal variances

# 4. Outliers (Boxplots)
ggplot(data, aes(x=meditator_factor, y=stress_score, fill=meditator_factor)) +
  geom_boxplot(alpha=0.6) +
  geom_jitter(width=0.1, alpha=0.3) +
  labs(title="Stress Levels by Meditation Practice",
       x="Group", y="Perceived Stress Score(0-40)") +
  scale_fill_brewer(palette="Set2") +
  theme_classic() +
  theme(legend.position="none")
# Look for outliers (points >1.5 IQR from box)

# === STEP 2: Descriptive Statistics ===
cat("\n=== Descriptive Statistics ===\n")
describe.by(data$stress_score, group=data$meditator_factor)

# === STEP 3: Compute Point-Biserial Correlation ===
cat("\n=== Point-Biserial Correlation ===\n")

# Method 1: Using cor.test (treats as Pearson r)
rpb_result <- cor.test(data$meditator, data$stress_score, method="pearson")
print(rpb_result)

cat(sprintf("\nr_pb = %.3f, 95%% CI [%.3f, %.3f], p = %.4f\n",
            rpb_result$estimate,
            rpb_result$conf.int[1],
            rpb_result$conf.int[2],
            rpb_result$p.value))

# Method 2: Manual calculation (for understanding)
# r_pb = (M₁ - M₀) / SD_total * sqrt(p * q)
# where p = proportion in group 1, q = 1-p
M0 <- mean(data$stress_score[data$meditator==0])
M1 <- mean(data$stress_score[data$meditator==1])
SD_total <- sd(data$stress_score)
p <- mean(data$meditator)  # Proportion in group 1
q <- 1 - p

rpb_manual <- ((M1 - M0) / SD_total) * sqrt(p * q)
cat(sprintf("\nManual r_pb calculation: %.3f\n", rpb_manual))

# === STEP 4: Relationship to Independent t-test ===
cat("\n=== Independent t-test(equivalent test) ===\n")
ttest_result <- t.test(stress_score ~ meditator_factor, data=data, var.equal=TRUE)
print(ttest_result)

# Convert t to r_pb: r_pb = sqrt(t² / (t² + df))
t_value <- ttest_result$statistic
df <- ttest_result$parameter
rpb_from_t <- sqrt(t_value^2 / (t_value^2 + df))
cat(sprintf("\nr_pb from t-test: %.3f (matches cor.test)\n", rpb_from_t))

# === STEP 5: Effect Sizes ===
cat("\n=== Effect Sizes ===\n")

# r_pb and r_pb²
cat(sprintf("Point-biserial r:       %.3f\n", rpb_result$estimate))
cat(sprintf("r_pb²:                  %.3f (%.1f%% variance explained)\n",
            rpb_result$estimate^2, 100*rpb_result$estimate^2))

# Cohen's d (standardized mean difference)
cohen_d <- cohen.d(stress_score ~ meditator_factor, data=data)
print(cohen_d)

# Interpret effect sizes
if(abs(rpb_result$estimate) < 0.1) {
  r_interpretation <- "negligible"
} else if(abs(rpb_result$estimate) < 0.3) {
  r_interpretation <- "small"
} else if(abs(rpb_result$estimate) < 0.5) {
  r_interpretation <- "medium"
} else {
  r_interpretation <- "large"
}

cat(sprintf("\nInterpretation: %s effect(Cohen, 1988)\n", r_interpretation))

# Relationship between r_pb and Cohen's d
# When groups equal size: d ≈ 2r / sqrt(1 - r²)
d_from_r <- 2 * rpb_result$estimate / sqrt(1 - rpb_result$estimate^2)
cat(sprintf("Cohen's d from r_pb:    %.3f\n", d_from_r))

# === STEP 6: Visualize Results ===

# Violin plot with individual points
ggplot(data, aes(x=meditator_factor, y=stress_score, fill=meditator_factor)) +
  geom_violin(alpha=0.5, trim=FALSE) +
  geom_boxplot(width=0.15, alpha=0.7) +
  geom_jitter(width=0.05, alpha=0.3) +
  stat_summary(fun=mean, geom="point", size=3, color="red", shape=18) +
  stat_summary(fun=mean, geom="text", aes(label=sprintf("M=%.1f", ..y..)), 
               vjust=-1, color="red", size=4) +
  labs(title=sprintf("Meditation and Stress: r_pb = %.3f, p = %.3f", 
                     rpb_result$estimate, rpb_result$p.value),
       x="Meditation Practice", y="Perceived Stress Score") +
  scale_fill_brewer(palette="Set2") +
  theme_classic() +
  theme(legend.position="none")

# === STEP 7: Bootstrap 95% CI (sensitivity check) ===
library(boot)

boot_rpb <- function(data, indices) {
  d <- data[indices,]
  cor(d$meditator, d$stress_score)
}

set.seed(2025)
boot_results <- boot(data, boot_rpb, R=1000)
boot_ci <- boot.ci(boot_results, type="perc")
cat("\n=== Bootstrap 95% CI ===\n")
print(boot_ci)

# === APA-Style Reporting ===
cat("\n=== APA-Style Report ===\n")
cat(sprintf(
  "A point-biserial correlation was computed to assess the relationship between\nmeditation practice(yes/no) and perceived stress. Assumptions were checked:\nnormality was confirmed within each group(Shapiro-Wilk p > .05), homogeneity\nof variance was satisfied(Levene's test, p = %.3f), and no extreme outliers\nwere detected.\n\nThere was a significant negative association between meditation practice and\nstress, r_pb = %.3f, p = %.3f, 95%% CI [%.3f, %.3f]. Meditation practice\nexplained %.1f%% of the variance in stress scores(r²_pb = %.3f). Meditators\n(M = %.2f, SD = %.2f) reported significantly lower stress than non-meditators\n(M = %.2f, SD = %.2f), t(%d) = %.2f, p = %.3f, Cohen's d = %.2f. This represents\na %s effect size(Cohen, 1988), suggesting meditation is associated with\nmeaningfully lower stress levels.\n",
  levene_result$`Pr(>F)`[1],
  rpb_result$estimate, rpb_result$p.value,
  rpb_result$conf.int[1], rpb_result$conf.int[2],
  100*rpb_result$estimate^2, rpb_result$estimate^2,
  M1, sd(data$stress_score[data$meditator==1]),
  M0, sd(data$stress_score[data$meditator==0]),
  ttest_result$parameter, ttest_result$statistic, ttest_result$p.value,
  cohen_d$estimate,
  r_interpretation
))