# ====================================================================
# Week 8: Regression Analysis
# Kmex Consult - R Data Science Course
# Open in RStudio and run each block with Ctrl+Enter (Cmd+Enter on Mac).
# ====================================================================

# Understanding and modeling relationships between variables

# ====================================================================
# WEEK 8 OVERVIEW
# ====================================================================
# Regression is how we model relationships between variables.
# KEY CONCEPTS:
# - Correlation: How variables move together (-1 to +1)
# - Regression: Predicting one variable from another
# - Causality: X causes Y (much harder to establish)
# THE EQUATION:
# y = β₀ + β₁*x + ε
# Where:
# - y = outcome (what we predict)
# - x = predictor (what we use)
# - β₀ = intercept (y when x=0)
# - β₁ = slope (change in y per unit x)
# - ε = error (random variation)
# THIS WEEK YOU'LL LEARN:
# ✓ Correlation coefficients
# ✓ Testing correlation significance
# ✓ Simple linear regression
# ✓ Multiple regression
# ✓ Interpreting regression output
# ✓ Model diagnostics
# ✓ Logistic regression (binary outcome)
# ✓ Assumptions and violations
# ✓ Causality vs correlation

# ====================================================================
# SESSION 1: CORRELATION AND SIMPLE REGRESSION
# ====================================================================

# --- Duration: 2 hours ---

# --- 1.1 Correlation Analysis ---
# CORRELATION - ASSOCIATION BETWEEN VARIABLES:
# Pearson Correlation (for continuous variables):


x <- c(1, 2, 3, 4, 5, 6)
y <- c(2, 4, 5, 4, 5, 7)

# Calculate correlation
r <- cor(x, y)
print(r) # [1] 0.769

# Interpretation:
# r = 1: Perfect positive (as x increases, y increases)
# r = 0: No correlation
# r = -1: Perfect negative (as x increases, y decreases)

# -1 to -0.7: Strong negative
# -0.7 to -0.3: Moderate negative
# -0.3 to 0: Weak negative
# 0 to 0.3: Weak positive
# 0.3 to 0.7: Moderate positive
# 0.7 to 1: Strong positive

# TEST SIGNIFICANCE OF CORRELATION:


# H₀: ρ = 0 (no correlation in population)
# H₁: ρ ≠ 0 (correlation exists)

cor.test(x, y)

# Output:
# Pearson's product-moment correlation
# t = 2.3381, df = 4, p-value = 0.1096
# cor = 0.769

# Interpretation: p = 0.11 > 0.05
# Weak evidence against H₀
# Could be due to small sample (n=6)

# With larger sample:
set.seed(42)
x_large <- rnorm(100, mean=50, sd=10)
y_large <- 0.5*x_large + rnorm(100, sd=10)
cor.test(x_large, y_large)
# p < 0.001 (significant!)

# CORRELATION MATRIX:


df <- tibble(
age = c(25, 30, 28, 35, 32),
salary = c(50000, 65000, 60000, 75000, 70000),
years_exp = c(2, 8, 6, 12, 10)
)

# Correlation between all numeric variables
cor(df)
# age salary years_exp
# age 1.0 0.9925408 0.9884393
# salary 0.9 1.0000000 0.9966555
# years_exp 0.9 0.9966555 1.0000000

# Visualize with heatmap
library(corrplot)
corrplot(cor(df), method="circle", type="upper")

# IMPORTANT: CORRELATION ≠ CAUSATION!


# Example: Ice cream sales vs drowning deaths
# Both increase in summer
# Correlation is high
# But ice cream doesn't cause drowning
# (Both related to warm temperature)

# Simpson's Paradox: Trend reverses when stratified

# Example: Admission rates by gender
# Overall: More men admitted
# Within departments: More women admitted!
# (Women applied to harder departments)

# KEY: Always investigate confounders

# SPEARMAN CORRELATION (For ranks/non-normal):


# When data violates normality assumption

x <- c(1, 2, 3, 4, 5)
y <- c(2, 4, 5, 4, 7) # Not perfectly linear

# Spearman (rank-based)
cor(x, y, method="spearman")

# More robust to outliers
x_outlier <- c(1, 2, 3, 4, 100)
y_outlier <- c(2, 4, 5, 4, 7)

cor(x_outlier, y_outlier, method="pearson") # Affected
# [1] 0.6149

cor(x_outlier, y_outlier, method="spearman") # Robust
# [1] 0.7

# SIMPLE LINEAR REGRESSION:


# Predict salary from years of experience

df <- tibble(
years_exp = c(1, 2, 3, 5, 8, 10),
salary = c(30000, 35000, 40000, 50000, 65000, 80000)
)

# Fit model
model <- lm(salary ~ years_exp, data = df)

# Get results
summary(model)

# Output shows:
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 14857 3447 4.312 0.0107 *
# years_exp 6642 558 11.899 0.0003 ***

# Interpretation:
# Intercept = 14857 (salary at 0 years)
# Slope = 6642 (salary increase per year)
# Model: salary = 14857 + 6642*years_exp

# R² = 0.9727 (explains 97% of variation)
# This is excellent fit!

# Predict new values
predict(model, newdata = tibble(years_exp = 7))
# [1] 61451 (predicted salary for 7 years experience)

# With confidence interval
predict(model, newdata = tibble(years_exp = 7),
interval = "confidence", level = 0.95)
# fit lwr upr
# 1 61451 58420 64482

# Visualize
ggplot(df, aes(x = years_exp, y = salary)) +
geom_point(size = 3) +
geom_smooth(method = "lm", se = TRUE) +
labs(title = "Salary vs Experience",
x = "Years Experience", y = "Salary ($)")

# MODEL DIAGNOSTICS:


# Check if regression assumptions hold

# 1. Residuals normally distributed
shapiro.test(resid(model)) # p > 0.05 → OK

# 2. Constant variance (homoscedasticity)
plot(model) # Shows residual plot

# 3. Independence (relevant for time series)
# Check design, not data

# 4. Linearity
# Visual inspection via scatter plot

# Identify problematic observations
influential_obs <- cooks.distance(model)
plot(influential_obs, type = "h")
# Points beyond Cook's distance threshold need investigation

# MULTIPLE REGRESSION:


# Predict salary from multiple variables

df <- tibble(
salary = c(50000, 65000, 60000, 75000, 70000),
years_exp = c(2, 8, 6, 12, 10),
education = c(16, 18, 16, 20, 18), # Years schooling
age = c(25, 45, 40, 50, 48)
)

# Fit model
model_multi <- lm(salary ~ years_exp + education + age, data = df)

summary(model_multi)

# Interpretation:
# Each coefficient = effect controlling for others
# years_exp coefficient = salary change per year (holding education/age constant)