# ============================================================================
#  04-11  Correlation & Regression Analyzer
#  Chapters 13-01, 13-02, 13-03
#
#  Run from the project folder:   Rscript r/analysis.R
# ============================================================================

h <- read.csv("data/houses.csv", stringsAsFactors = TRUE)
n <- nrow(h)

# ── 1. CORRELATION ──────────────────────────────────────────────────────────
cat("\n===== 1. CORRELATION MATRIX =====\n")

num <- h[, c("sqft", "bedrooms", "age_years", "price")]
cat("\nPearson r:\n");  print(round(cor(num), 4))
cat("\nSpearman rho (rank-based, robust to curvature and outliers):\n")
print(round(cor(num, method = "spearman"), 4))

cat("\nSignificance of each correlation with price:\n")
for (v in c("sqft", "bedrooms", "age_years")) {
  ct <- cor.test(num[[v]], num$price)
  cat(sprintf("  %-10s r = %+.4f   r^2 = %.4f   t = %7.3f   p = %.3e\n",
              v, ct$estimate, ct$estimate^2, ct$statistic, ct$p.value))
}

cat("\nWATCH THE PREDICTOR-PREDICTOR CORRELATIONS:\n")
cat(sprintf("  sqft vs bedrooms  r = %.4f", cor(h$sqft, h$bedrooms)))
cat(if (abs(cor(h$sqft, h$bedrooms)) > 0.7)
      "   <- high; expect multicollinearity in step 7\n" else "\n")

# ── 2. SIMPLE REGRESSION ────────────────────────────────────────────────────
cat("\n===== 2. SIMPLE REGRESSION — price on sqft =====\n")

m1 <- lm(price ~ sqft, data = h)
s1 <- summary(m1)
b0 <- coef(m1)[1]; b1 <- coef(m1)[2]

cat(sprintf("\nFITTED LINE:  price-hat = %.2f + %.4f * sqft\n", b0, b1))
cat(sprintf("\nSLOPE      each additional square foot is associated with a\n"))
cat(sprintf("           $%.2f increase in predicted price ($%.0f per 100 sqft)\n",
            b1, 100 * b1))
cat(sprintf("INTERCEPT  $%.2f at sqft = 0 -- far outside the observed range\n", b0))
cat(sprintf("           (%.0f to %.0f sqft), so it is a mathematical anchor only\n",
            min(h$sqft), max(h$sqft)))
cat(sprintf("r^2        %.4f  -> %.1f%% of price variation explained by floor area\n",
            s1$r.squared, 100 * s1$r.squared))
cat(sprintf("s_e        $%.2f  -> a typical prediction misses by about this much\n",
            s1$sigma))

# ── 3. INFERENCE ON THE SLOPE ───────────────────────────────────────────────
cat("\n===== 3. INFERENCE ON THE SLOPE =====\n")
cat("H0: beta1 = 0   (floor area has no linear predictive value)\n")
cat("H1: beta1 != 0\n\n")
print(round(s1$coefficients, 6))
cat(sprintf("\ncritical t(0.025, %d) = %.4f\n", m1$df.residual,
            qt(0.975, m1$df.residual)))
cat("\n95% confidence intervals for the coefficients:\n")
print(round(confint(m1), 4))
cat(sprintf("\n=> each extra square foot is worth between $%.2f and $%.2f,\n",
            confint(m1)[2, 1], confint(m1)[2, 2]))
cat("   with 95% confidence. That range is what a decision-maker needs.\n")

cat("\nThe three equivalent tests (simple regression only):\n")
cat(sprintf("  slope t          %.4f\n", s1$coefficients[2, 3]))
cat(sprintf("  correlation t    %.4f\n", cor.test(h$sqft, h$price)$statistic))
cat(sprintf("  sqrt(model F)    %.4f\n", sqrt(s1$fstatistic[1])))

cat("\nANOVA decomposition:\n"); print(anova(m1))

# ── 4. RESIDUAL DIAGNOSTICS (the LINE assumptions) ──────────────────────────
cat("\n===== 4. RESIDUAL DIAGNOSTICS =====\n")

r <- residuals(m1)
cat(sprintf("sum of residuals        %.10f  (must be ~0)\n", sum(r)))
cat(sprintf("Shapiro-Wilk on resid   p = %.4f  -> %s\n",
            shapiro.test(r)$p.value,
            ifelse(shapiro.test(r)$p.value > 0.05,
                   "NORMALITY looks fine", "normality is questionable")))
cor_rf <- cor(abs(r), fitted(m1))
cat(sprintf("cor(|resid|, fitted)    %+.4f  -> %s\n", cor_rf,
            ifelse(abs(cor_rf) < 0.2, "EQUAL VARIANCE looks fine",
                   "possible heteroscedasticity -- inspect the plot")))

infl <- which(cooks.distance(m1) > 4 / n)
cat(sprintf("influential points (Cook's D > 4/n): %d\n", length(infl)))
if (length(infl)) print(h[infl, c("house_id", "sqft", "price")], row.names = FALSE)

# ── 5. CONFIDENCE vs PREDICTION INTERVAL ────────────────────────────────────
cat("\n===== 5. CONFIDENCE vs PREDICTION INTERVAL at sqft = 2000 =====\n")

nd <- data.frame(sqft = 2000)
ci <- predict(m1, nd, interval = "confidence")
pi <- predict(m1, nd, interval = "prediction")

cat(sprintf("point estimate            $%.2f\n", ci[1, "fit"]))
cat(sprintf("95%% CONFIDENCE interval   $%.2f to $%.2f   (width $%.2f)\n",
            ci[1, "lwr"], ci[1, "upr"], ci[1, "upr"] - ci[1, "lwr"]))
cat(sprintf("95%% PREDICTION interval   $%.2f to $%.2f   (width $%.2f)\n",
            pi[1, "lwr"], pi[1, "upr"], pi[1, "upr"] - pi[1, "lwr"]))
cat(sprintf("\nThe prediction interval is %.1fx wider, because it carries BOTH the\n",
            (pi[1,"upr"] - pi[1,"lwr"]) / (ci[1,"upr"] - ci[1,"lwr"])))
cat("uncertainty in the line AND the house-to-house scatter.\n")
cat("Quote the CONFIDENCE interval for 'the average 2000 sqft house';\n")
cat("quote the PREDICTION interval for 'this particular house'.\n")

# ── 6. MULTIPLE REGRESSION ──────────────────────────────────────────────────
cat("\n===== 6. MULTIPLE REGRESSION =====\n")

m2 <- lm(price ~ sqft + bedrooms, data = h)
m3 <- lm(price ~ sqft + bedrooms + age_years, data = h)
m4 <- lm(price ~ sqft + bedrooms + age_years + garage, data = h)

cat("\nFull model (m4):\n"); print(round(summary(m4)$coefficients, 4))
cat(sprintf("\nR^2 = %.4f    Adjusted R^2 = %.4f    s_e = %.2f\n",
            summary(m4)$r.squared, summary(m4)$adj.r.squared, summary(m4)$sigma))
cat(sprintf("F(%d, %d) = %.3f, p = %.3e\n",
            summary(m4)$fstatistic[2], summary(m4)$fstatistic[3],
            summary(m4)$fstatistic[1],
            pf(summary(m4)$fstatistic[1], summary(m4)$fstatistic[2],
               summary(m4)$fstatistic[3], lower.tail = FALSE)))

cat("\nHOW TO READ EACH COEFFICIENT: the change in predicted price per\n")
cat("one-unit rise in that predictor, HOLDING THE OTHERS CONSTANT.\n")
cat("'garageYes' is the difference between a house with a garage and one\n")
cat("without, at the same sqft, bedrooms and age. (No = the reference level.)\n")

# ── 7. MODEL COMPARISON ─────────────────────────────────────────────────────
cat("\n===== 7. MODEL COMPARISON =====\n")

mods <- list(m1 = m1, m2 = m2, m3 = m3, m4 = m4)
print(data.frame(
  model  = names(mods),
  k      = sapply(mods, function(m) length(coef(m)) - 1),
  r2     = round(sapply(mods, function(m) summary(m)$r.squared), 4),
  adj_r2 = round(sapply(mods, function(m) summary(m)$adj.r.squared), 4),
  s_e    = round(sapply(mods, function(m) summary(m)$sigma), 2),
  AIC    = round(sapply(mods, AIC), 2)
), row.names = FALSE)

cat("\nCompare with ADJUSTED R^2 or AIC, never with raw R^2 --\n")
cat("raw R^2 can only rise as predictors are added.\n")
cat("\nNested F-tests:\n"); print(anova(m1, m2, m3, m4))

# ── 8. MULTICOLLINEARITY ────────────────────────────────────────────────────
cat("\n===== 8. MULTICOLLINEARITY (VIF) =====\n")

vif_manual <- function(model) {
  preds <- attr(terms(model), "term.labels")
  sapply(preds, function(p) {
    others <- setdiff(preds, p)
    if (!length(others)) return(1)
    f <- as.formula(paste(p, "~", paste(others, collapse = " + ")))
    1 / (1 - summary(lm(f, data = model$model))$r.squared)
  })
}
v <- vif_manual(m4)
print(round(v, 3))
cat("\nVIF < 5 fine   5-10 investigate   > 10 serious\n")
if (any(v > 5)) {
  cat("=> ", paste(names(v)[v > 5], collapse = ", "),
      " show inflated variance. A non-significant t here may reflect\n", sep = "")
  cat("   collinearity, NOT irrelevance -- check before dropping anything.\n")
} else {
  cat("=> No serious multicollinearity.\n")
}

# ── PLOTS ───────────────────────────────────────────────────────────────────
png("regression_plots.png", width = 1200, height = 900)
par(mfrow = c(2, 3), mar = c(4.5, 4.5, 3, 1))

plot(h$sqft, h$price, pch = 19, col = "#5B2A86",
     main = sprintf("price vs sqft  (r = %.3f)", cor(h$sqft, h$price)),
     xlab = "sqft", ylab = "price")
abline(m1, col = "#0FA3A3", lwd = 2)
newx <- data.frame(sqft = seq(min(h$sqft), max(h$sqft), length.out = 100))
pc <- predict(m1, newx, interval = "confidence")
pp <- predict(m1, newx, interval = "prediction")
lines(newx$sqft, pc[, 2], col = "#0FA3A3", lty = 2)
lines(newx$sqft, pc[, 3], col = "#0FA3A3", lty = 2)
lines(newx$sqft, pp[, 2], col = "#8A5FBF", lty = 3)
lines(newx$sqft, pp[, 3], col = "#8A5FBF", lty = 3)
legend("topleft", c("fit", "95% CI", "95% PI"), lty = 1:3,
       col = c("#0FA3A3", "#0FA3A3", "#8A5FBF"), bty = "n")

plot(fitted(m1), r, pch = 19, col = "#5B2A86",
     main = "Residuals vs fitted", xlab = "fitted", ylab = "residual")
abline(h = 0, lty = 2, col = "grey40")

qqnorm(r, pch = 19, col = "#5B2A86", main = "Normal Q-Q of residuals")
qqline(r, col = "#0FA3A3", lwd = 2)

hist(r, breaks = 20, col = "#0B7A7A", border = "white",
     main = "Residual histogram", xlab = "residual")

boxplot(price ~ garage, data = h, col = c("#5B2A86", "#0FA3A3"),
        main = "Price by garage", ylab = "price")

plot(fitted(m4), residuals(m4), pch = 19, col = "#0B7A7A",
     main = "Full model: residuals vs fitted", xlab = "fitted", ylab = "residual")
abline(h = 0, lty = 2, col = "grey40")

par(mfrow = c(1, 1))
dev.off()
cat("\nWrote regression_plots.png\n")