Partial Correlation

# Partial Correlation: Meditation-Stress controlling for Sleep
# Research: Does meditation predict stress after accounting for sleep quality?

# Install/load packages
library(ppcor)        # For partial correlation
library(psych)        # For pairs.panels
library(car)          # For vif
library(MVN)          # For multivariate normality tests
library(tidyverse)

# Simulate realistic data (or load: data <- read.csv("meditation_stress.csv"))
set.seed(2025)
library(MASS)

# Create correlated variables with known structure:
# Sleep (Z) causes both Meditation (X) and Stress (Y)
# Meditation (X) causes Stress (Y)
mu <- c(150, 20, 10)  # Means: meditation=150min, stress=20, sleep=10
Sigma <- matrix(c(
  2500,  -30,  -50,   # Meditation variance & covariances
   -30,   64,   15,   # Stress variance & covariances  
   -50,   15,   25    # Sleep variance & covariances
), nrow=3, byrow=TRUE)

data_matrix <- mvrnorm(n=120, mu=mu, Sigma=Sigma)
data <- data.frame(
  meditation_min = pmax(0, pmin(300, data_matrix[,1])),  # Bounded 0-300
  stress_score = pmax(0, pmin(40, data_matrix[,2])),     # Bounded 0-40
  sleep_quality = pmax(0, pmin(21, data_matrix[,3]))     # Bounded 0-21
)

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

# 1. Linearity: Scatterplot matrix
pairs.panels(data, 
             method = "pearson",
             hist.col = "steelblue",
             density = TRUE,
             ellipses = TRUE,
             main = "Scatterplot Matrix: Linearity Check")
# Look for: linear patterns (not curved), correlation ellipses

# 2. Multivariate normality
cat("\n=== Multivariate Normality Tests ===\n")
mvn_test <- mvn(data, mvnTest = "mardia")
print(mvn_test$multivariateNormality)
# Mardia skewness and kurtosis p > .05 → multivariate normal

# Univariate normality for each variable
cat("\n=== Univariate Normality(Shapiro-Wilk) ===\n")
shapiro.test(data$meditation_min)
shapiro.test(data$stress_score) 
shapiro.test(data$sleep_quality)

# Q-Q plots
par(mfrow=c(1,3))
qqnorm(data$meditation_min, main="Meditation Q-Q Plot")
qqline(data$meditation_min)
qqnorm(data$stress_score, main="Stress Q-Q Plot")
qqline(data$stress_score)
qqnorm(data$sleep_quality, main="Sleep Q-Q Plot")
qqline(data$sleep_quality)

# 3. Multivariate outliers (Mahalanobis distance)
mahal <- mahalanobis(data, colMeans(data), cov(data))
mahal_cutoff <- qchisq(0.999, df=3)  # χ² critical value, df = # variables
outliers <- which(mahal > mahal_cutoff)
cat("\n=== Multivariate Outliers ===\n")
cat("Number of outliers(p < .001):", length(outliers), "\n")
if(length(outliers) > 0) {
  cat("Outlier cases:", outliers, "\n")
}

# 4. Homoscedasticity: residual plots from Y~Z and X~Z
par(mfrow=c(1,2))
fit_yz <- lm(stress_score ~ sleep_quality, data=data)
plot(fit_yz, which=1, main="Y~Z Residuals(Stress~Sleep)")

fit_xz <- lm(meditation_min ~ sleep_quality, data=data)
plot(fit_xz, which=1, main="X~Z Residuals(Meditation~Sleep)")
# Look for: constant spread (no funnel shape)

# === STEP 2: Compute Zero-Order (Bivariate) Correlation ===
cat("\n=== Zero-Order Correlation(before controlling) ===\n")
zero_order <- cor.test(data$meditation_min, data$stress_score)
print(zero_order)
# r_XY = correlation ignoring sleep

# === STEP 3: Compute Partial Correlation ===
cat("\n=== Partial Correlation(controlling for Sleep) ===\n")

# Method 1: Using ppcor package (easiest)
partial_result <- pcor.test(data$meditation_min, 
                             data$stress_score, 
                             data$sleep_quality)
print(partial_result)
# r_XY·Z = partial correlation after removing sleep effects

# Method 2: Manual calculation (for understanding)
# Step 1: Regress Y on Z, get residuals (Y with Z removed)
resid_y <- residuals(lm(stress_score ~ sleep_quality, data=data))

# Step 2: Regress X on Z, get residuals (X with Z removed)  
resid_x <- residuals(lm(meditation_min ~ sleep_quality, data=data))

# Step 3: Correlate the residuals
partial_manual <- cor.test(resid_x, resid_y)
cat("\nManual partial correlation(residual method):\n")
print(partial_manual)

# Visualize residual correlation
ggplot(data.frame(resid_x, resid_y), aes(x=resid_x, y=resid_y)) +
  geom_point(alpha=0.5) +
  geom_smooth(method="lm", se=TRUE, color="red") +
  labs(title="Partial Correlation Visualization",
       subtitle="Meditation-Stress residuals after controlling Sleep",
       x="Meditation residuals(sleep removed)",
       y="Stress residuals(sleep removed)") +
  theme_classic()

# === STEP 4: Compare Zero-Order vs Partial ===
cat("\n=== Comparison ===\n")
cat(sprintf("Zero-order r_XY:        %.3f (p = %.4f)\n", 
            zero_order$estimate, zero_order$p.value))
cat(sprintf("Partial r_XY·Z:         %.3f (p = %.4f)\n", 
            partial_result$estimate, partial_result$p.value))
cat(sprintf("Change in r:            %.3f\n", 
            partial_result$estimate - zero_order$estimate))
cat(sprintf("Proportion explained by sleep: %.1f%%\n",
            100 * (zero_order$estimate - partial_result$estimate) / zero_order$estimate))

# === STEP 5: Effect Size ===
cat("\n=== Effect Sizes ===\n")
cat(sprintf("Partial r²:             %.3f (%.1f%% unique variance)\n",
            partial_result$estimate^2, 
            100*partial_result$estimate^2))

# Interpret partial r (Cohen's guidelines)
if(abs(partial_result$estimate) < 0.1) {
  interpretation <- "negligible"
} else if(abs(partial_result$estimate) < 0.3) {
  interpretation <- "small"  
} else if(abs(partial_result$estimate) < 0.5) {
  interpretation <- "medium"
} else {
  interpretation <- "large"
}
cat(sprintf("Interpretation:         %s effect\n", interpretation))

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

boot_partial <- function(data, indices) {
  d <- data[indices,]
  pcor.test(d$meditation_min, d$stress_score, d$sleep_quality)$estimate
}

set.seed(2025)
boot_results <- boot(data, boot_partial, 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 partial correlation was computed to assess the relationship between\nmeditation practice and perceived stress while controlling for sleep quality.\nAssumptions were checked: multivariate normality was satisfied(Mardia test\np > .05), linearity was confirmed via scatterplot matrix, and no extreme\nmultivariate outliers were detected(Mahalanobis D² < χ²_critical).\n\nThe zero-order correlation between meditation and stress was r = %.3f, p = %.3f.\nAfter controlling for sleep quality, the partial correlation was r_partial = %.3f,\np = %.3f, 95%% CI [%.3f, %.3f]. Sleep quality accounted for %.1f%% of the\noriginal correlation. The partial correlation r² = %.3f indicates that meditation\nexplains %.1f%% of unique variance in stress, beyond sleep quality. This represents\na %s effect size(Cohen, 1988).\n",
  zero_order$estimate, zero_order$p.value,
  partial_result$estimate, partial_result$p.value,
  boot_ci$percent[4], boot_ci$percent[5],
  100 * (zero_order$estimate - partial_result$estimate) / zero_order$estimate,
  partial_result$estimate^2,
  100 * partial_result$estimate^2,
  interpretation
))