Pearson Correlation (r)

# Pearson Correlation: Height and Weight
# Classic example of linear association between continuous variables

library(tidyverse)
library(ggplot2)
library(car)         # For Q-Q plots and outlier diagnostics
library(psych)       # For describe and corr.test
library(DescTools)   # For CIs

# Simulate realistic data (or load: data <- read.csv("height_weight.csv"))
set.seed(2025)
n <- 100

# Height: normal distribution, mean=170cm (mix male/female), SD=10cm
height <- rnorm(n, mean=170, sd=10)

# Weight: linearly related to height with noise
# True relationship: weight ≈ -100 + 1.0*height + error
weight <- -100 + 1.0 * height + rnorm(n, mean=0, sd=5)

data <- data.frame(
  subject_id = 1:n,
  height = height,
  weight = weight
)

head(data, 10)
psych::describe(data[, c("height", "weight")])

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

cat("\n=== Assumption Checks ===\n")

# 1. Linearity: Scatterplot with regression line
ggplot(data, aes(x=height, y=weight)) +
  geom_point(alpha=0.6, size=2.5) +
  geom_smooth(method="lm", color="blue", se=TRUE, linewidth=1) +
  geom_smooth(method="loess", color="red", se=FALSE, linetype="dashed") +
  labs(title="Height vs Weight: Linearity Check",
       subtitle="Blue=linear fit, Red=loess smooth(should overlap if linear)",
       x="Height(cm)",
       y="Weight(kg)") +
  theme_classic(base_size=12)

cat("Check: Linear(blue) and loess(red) lines should be similar.\n")

# 2. Normality: Q-Q plots
par(mfrow=c(1, 2))
qqPlot(data$height, main="Q-Q Plot: Height", ylab="Height(cm)")
qqPlot(data$weight, main="Q-Q Plot: Weight", ylab="Weight(kg)")
par(mfrow=c(1, 1))

# Shapiro-Wilk test (formal test, but use with caution - sensitive to large n)
shapiro_height <- shapiro.test(data$height)
shapiro_weight <- shapiro.test(data$weight)
cat(sprintf("\nShapiro-Wilk normality tests:\n"))
cat(sprintf("  Height: W=%.3f, p=%.3f %s\n", 
            shapiro_height$statistic, shapiro_height$p.value,
            ifelse(shapiro_height$p.value > 0.05, "(normal)", "(non-normal)")))
cat(sprintf("  Weight: W=%.3f, p=%.3f %s\n", 
            shapiro_weight$statistic, shapiro_weight$p.value,
            ifelse(shapiro_weight$p.value > 0.05, "(normal)", "(non-normal)")))

# 3. Outliers: Cook's distance and influence
model <- lm(weight ~ height, data=data)
cooks_d <- cooks.distance(model)
outlier_threshold <- 4 / n
outliers <- which(cooks_d > outlier_threshold)

cat(sprintf("\nOutliers(Cook's D > %.3f): %d cases\n", 
            outlier_threshold, length(outliers)))
if (length(outliers) > 0) {
  cat("Outlier IDs:", outliers, "\n")
  print(data[outliers, ])
}

# Scatterplot with outliers highlighted
data$outlier <- cooks_d > outlier_threshold
ggplot(data, aes(x=height, y=weight, color=outlier)) +
  geom_point(size=3, alpha=0.7) +
  scale_color_manual(values=c("FALSE"="black", "TRUE"="red")) +
  labs(title="Outlier Detection: Cook's Distance",
       subtitle=sprintf("%d outliers(red points)", length(outliers))) +
  theme_classic()

# 4. Homoscedasticity: Residual plot
data$residuals <- residuals(model)
data$fitted <- fitted(model)

ggplot(data, aes(x=fitted, y=residuals)) +
  geom_point(alpha=0.6) +
  geom_hline(yintercept=0, color="red", linetype="dashed") +
  geom_smooth(se=FALSE, color="blue") +
  labs(title="Residuals vs Fitted: Homoscedasticity Check",
       subtitle="Should show random scatter around zero",
       x="Fitted Values",
       y="Residuals") +
  theme_classic()

# === STEP 2: Compute Pearson Correlation ===

cat("\n=== Pearson Correlation Results ===\n")

# Method 1: Base R cor.test
result <- cor.test(data$height, data$weight, 
                   method="pearson", 
                   alternative="two.sided")
print(result)

cat(sprintf("\nr = %.3f\n", result$estimate))
cat(sprintf("95%% CI: [%.3f, %.3f]\n", result$conf.int[1], result$conf.int[2]))
cat(sprintf("t(%d) = %.2f, p %s\n", 
            result$parameter,
            result$statistic, 
            ifelse(result$p.value < 0.001, "< .001", 
                   sprintf("= %.4f", result$p.value))))

# r-squared (proportion of variance explained)
r_squared <- result$estimate^2
cat(sprintf("\nr² = %.3f (%.1f%% of variance in weight explained by height)\n",
            r_squared, r_squared * 100))

# === STEP 3: Effect Size Interpretation ===

cat("\n=== Effect Size Guidelines(Cohen, 1988) ===\n")
cat("r = 0.10: small, r = 0.30: medium, r = 0.50: large\n\n")

r_val <- result$estimate
if (abs(r_val) < 0.10) {
  strength <- "negligible"
} else if (abs(r_val) < 0.30) {
  strength <- "small"
} else if (abs(r_val) < 0.50) {
  strength <- "medium"
} else {
  strength <- "large"
}

direction <- ifelse(r_val > 0, "positive", "negative")
cat(sprintf("Observed effect: %s(%s association)\n", strength, direction))

# === STEP 4: Bootstrap Confidence Interval (alternative to parametric CI) ===

cat("\n=== Bootstrap 95% CI(n=1000 resamples) ===\n")
set.seed(2025)
boot_r <- replicate(1000, {
  indices <- sample(1:n, n, replace=TRUE)
  cor(data$height[indices], data$weight[indices])
})

boot_ci <- quantile(boot_r, c(0.025, 0.975))
cat(sprintf("Bootstrap CI: [%.3f, %.3f]\n", boot_ci[1], boot_ci[2]))
cat("Compare with parametric CI - should be similar if normality holds.\n")

# === STEP 5: Compare with Spearman (robustness check) ===

spearman_result <- cor.test(data$height, data$weight, method="spearman")
cat(sprintf("\nSpearman's ρ = %.3f (p %s)\n",
            spearman_result$estimate,
            ifelse(spearman_result$p.value < 0.001, "< .001", 
                   sprintf("= %.4f", spearman_result$p.value))))
cat("Note: Pearson and Spearman should be similar if linearity holds.\n")
cat(sprintf("Difference: %.3f (large differences suggest non-linearity/outliers)\n",
            abs(r_val - spearman_result$estimate)))

# === STEP 6: Sensitivity Analysis - Remove Outliers ===

if (length(outliers) > 0) {
  cat("\n=== Sensitivity Analysis: Excluding Outliers ===\n")
  data_clean <- data[!data$outlier, ]
  r_clean <- cor.test(data_clean$height, data_clean$weight, method="pearson")$estimate
  cat(sprintf("r without outliers = %.3f (original r = %.3f)\n", r_clean, r_val))
  cat(sprintf("Difference = %.3f (>0.10 suggests outliers influential)\n", 
              abs(r_clean - r_val)))
} else {
  cat("\nNo outliers to remove - sensitivity analysis not needed.\n")
}

# === STEP 7: Visualize Correlation Strength ===

ggplot(data[!data$outlier, ], aes(x=height, y=weight)) +
  geom_point(alpha=0.5, size=2) +
  geom_smooth(method="lm", color="blue", fill="lightblue") +
  labs(title=sprintf("Height-Weight Correlation: r = %.2f", r_val),
       subtitle=sprintf("r² = %.2f: %.0f%% variance explained", 
                        r_squared, r_squared*100),
       x="Height(cm)",
       y="Weight(kg)") +
  annotate("text", x=min(data$height)+5, y=max(data$weight)-5,
           label=sprintf("r = %.2f***\n95%% CI [%.2f, %.2f]",
                        r_val, result$conf.int[1], result$conf.int[2]),
           hjust=0, size=5) +
  theme_classic(base_size=12)

# === APA-Style Reporting ===

cat("\n=== APA-Style Report ===\n")
cat(sprintf(
"A Pearson product-moment correlation was computed to assess the linear 
relationship between height and weight in a sample of %d adults. Preliminary 
analyses showed the relationship was linear with both variables approximately 
normally distributed(Shapiro-Wilk p > .05), and no influential outliers were 
detected(Cook's D < %.2f). There was a %s, significant positive correlation 
between height and weight, r(%d) = %.2f, 95%% CI [%.2f, %.2f], p < .001. 
Height explained %.0f%% of the variance in weight(r² = %.2f). These findings 
indicate that taller individuals tend to weigh more, consistent with 
anthropometric research. The large effect size suggests height is a strong 
predictor of weight in this population.\n",
n, 4/n, strength, result$parameter, r_val, 
result$conf.int[1], result$conf.int[2],
r_squared * 100, r_squared
))