# ============================================================================
#  04-12  Capstone — Business Statistics Project
#  Every chapter, on one data set.
#
#  Run from the project folder:   Rscript r/analysis.R
# ============================================================================

df <- read.csv("data/customers.csv", stringsAsFactors = TRUE)
n  <- nrow(df)
spend <- df$annual_spend

hr <- function(title) {
  cat("\n", strrep("=", 74), "\n", title, "\n", strrep("=", 74), "\n", sep = "")
}

# ════════════════════════════════════════════════════════════════════════════
hr("PART 1 — ORGANIZE & DESCRIBE")
# ════════════════════════════════════════════════════════════════════════════

cat("\n-- data quality --\n")
cat(sprintf("rows                 %d\n", n))
cat(sprintf("duplicate ids        %d\n", sum(duplicated(df$customer_id))))
cat(sprintf("missing values       %d\n", sum(is.na(df))))
cat(sprintf("annual_spend range   %.2f to %.2f\n", min(spend), max(spend)))
cat(sprintf("visits range         %d to %d\n", min(df$visits_per_year),
            max(df$visits_per_year)))
cat("\ncategory levels:\n")
print(sapply(df[sapply(df, is.factor)], function(x) paste(levels(x), collapse = ", ")))

cat("\n-- annual spend --\n")
q <- quantile(spend, c(0.25, 0.5, 0.75))
desc <- c(n = n, mean = mean(spend), median = median(spend), sd = sd(spend),
          cv_pct = 100 * sd(spend) / mean(spend),
          min = min(spend), Q1 = q[[1]], Q3 = q[[3]], max = max(spend),
          IQR = q[[3]] - q[[1]],
          skew_SK = 3 * (mean(spend) - median(spend)) / sd(spend))
print(round(desc, 3))

shape <- if (abs(desc[["skew_SK"]]) < 0.5) "roughly SYMMETRIC" else
         if (desc[["skew_SK"]] > 0) "RIGHT-SKEWED" else "LEFT-SKEWED"
cat(sprintf("\nshape: %s  ->  report the %s\n", shape,
            ifelse(shape == "roughly SYMMETRIC", "MEAN and SD", "MEDIAN and IQR")))

Q1 <- q[[1]]; Q3 <- q[[3]]; iqr <- Q3 - Q1
fences <- c(Q1 - 1.5 * iqr, Q3 + 1.5 * iqr)
out <- spend < fences[1] | spend > fences[2]
cat(sprintf("outlier fences: %.2f to %.2f   -> %d outlier(s)\n",
            fences[1], fences[2], sum(out)))

cat("\n-- by region --\n")
print(round(cbind(n    = table(df$region),
                  mean = tapply(spend, df$region, mean),
                  sd   = tapply(spend, df$region, sd),
                  median = tapply(spend, df$region, median)), 2))

cat("\n-- channel x satisfaction --\n")
print(addmargins(table(df$channel, df$satisfied)))

# ════════════════════════════════════════════════════════════════════════════
hr("PART 2 — ESTIMATE")
# ════════════════════════════════════════════════════════════════════════════

tt <- t.test(spend)
cat(sprintf("\nMean annual spend      $%.2f\n", mean(spend)))
cat(sprintf("95%% CI                 $%.2f to $%.2f\n", tt$conf.int[1], tt$conf.int[2]))
cat(sprintf("margin of error        $%.2f\n", diff(tt$conf.int) / 2))
cat("\n=> \"We are 95% confident that the mean annual spend of all customers\n")
cat(sprintf("   lies between $%.2f and $%.2f.\"\n", tt$conf.int[1], tt$conf.int[2]))

x <- sum(df$satisfied == "Yes"); ph <- x / n
se_p <- sqrt(ph * (1 - ph) / n)
ci_p <- ph + c(-1, 1) * qnorm(0.975) * se_p
cat(sprintf("\nProportion satisfied   %.4f  (%d of %d)\n", ph, x, n))
cat(sprintf("condition check        n*p = %.0f, n*(1-p) = %.0f  (both >= 5)\n",
            n * ph, n * (1 - ph)))
cat(sprintf("95%% CI                 %.4f to %.4f  (%.1f%% to %.1f%%)\n",
            ci_p[1], ci_p[2], 100 * ci_p[1], 100 * ci_p[2]))

need <- ceiling(ph * (1 - ph) * (qnorm(0.975) / 0.02)^2)
cat(sprintf("\nFor a +/-2%% margin we would need n = %d  (currently %d)\n", need, n))

# ════════════════════════════════════════════════════════════════════════════
hr("PART 3 — TEST")
# ════════════════════════════════════════════════════════════════════════════

# ---- TEST 1: one-sample t against the $250 target --------------------------
cat("\n--- TEST 1: does mean spend exceed the $250 target? ---\n")
cat("H0: mu <= 250      H1: mu > 250      right-tailed, alpha = 0.05\n\n")

t1 <- t.test(spend, mu = 250, alternative = "greater")
d1 <- (mean(spend) - 250) / sd(spend)
cat(sprintf("t(%d) = %.4f, p = %.6f, Cohen's d = %.4f\n",
            t1$parameter, t1$statistic, t1$p.value, d1))
cat(sprintf("two-sided 95%% CI: $%.2f to $%.2f\n",
            t.test(spend)$conf.int[1], t.test(spend)$conf.int[2]))
cat(sprintf("=> %s. Mean spend %s the $250 target.\n",
            ifelse(t1$p.value <= 0.05, "REJECT H0", "FAIL TO REJECT H0"),
            ifelse(t1$p.value <= 0.05, "EXCEEDS", "is not shown to exceed")))

# ---- TEST 2: one-way ANOVA across regions ---------------------------------
cat("\n--- TEST 2: does mean spend differ across regions? ---\n")
cat("H0: mu_N = mu_S = mu_E = mu_W      H1: at least one differs\n\n")

vars <- tapply(spend, df$region, var)
cat(sprintf("variance ratio (max/min) = %.3f  -> %s\n", max(vars) / min(vars),
            ifelse(max(vars) / min(vars) < 4, "pooled ANOVA is appropriate",
                   "unequal; prefer Welch")))
cat(sprintf("Bartlett test p = %.4f\n", bartlett.test(annual_spend ~ region, df)$p.value))

aov1 <- aov(annual_spend ~ region, data = df)
print(summary(aov1))
ss <- summary(aov1)[[1]][["Sum Sq"]]
Fp <- summary(aov1)[[1]][["Pr(>F)"]][1]
cat(sprintf("\neta squared = %.4f  -> %.1f%% of spending variation explained by region\n",
            ss[1] / sum(ss), 100 * ss[1] / sum(ss)))
cat(sprintf("Shapiro-Wilk on residuals: p = %.4f\n", shapiro.test(residuals(aov1))$p.value))

if (Fp <= 0.05) {
  cat("\nF IS significant -> post-hoc comparisons:\n")
  print(TukeyHSD(aov1))
} else {
  cat("\nF is NOT significant -> no post-hoc test.\n")
}

# ---- TEST 3: chi-square independence --------------------------------------
cat("\n--- TEST 3: are channel and satisfaction associated? ---\n")
cat("H0: channel and satisfaction are INDEPENDENT      H1: associated\n\n")

tab <- table(df$channel, df$satisfied)
ct <- chisq.test(tab)
cat("observed:\n"); print(addmargins(tab))
cat("\nexpected:\n"); print(round(ct$expected, 2))
cat(sprintf("\nall expected >= 5?  %s\n", all(ct$expected >= 5)))
cat("\ncell contributions:\n"); print(round(ct$residuals^2, 4))
cat(sprintf("\nchi-square = %.4f, df = %d, p = %.4f\n",
            ct$statistic, ct$parameter, ct$p.value))
V <- sqrt(as.numeric(ct$statistic) / (n * min(dim(tab) - 1)))
cat(sprintf("Cramer's V = %.4f\n", V))
cat("\nsatisfaction rate by channel:\n")
print(round(100 * prop.table(tab, 1), 1))
cat(sprintf("=> %s\n", ifelse(ct$p.value <= 0.05,
    "REJECT H0. Channel and satisfaction ARE associated.",
    "FAIL TO REJECT H0. No evidence of an association.")))

# ════════════════════════════════════════════════════════════════════════════
hr("PART 4 — MODEL")
# ════════════════════════════════════════════════════════════════════════════

cat("\ncorrelation of spend with visits:\n")
cat(sprintf("  r = %.4f, r^2 = %.4f, p = %.3e\n",
            cor(df$visits_per_year, spend), cor(df$visits_per_year, spend)^2,
            cor.test(df$visits_per_year, spend)$p.value))

cat("\n--- simple regression: annual_spend ~ visits_per_year ---\n")
s1 <- lm(annual_spend ~ visits_per_year, data = df)
print(round(summary(s1)$coefficients, 4))
cat(sprintf("R^2 = %.4f    s_e = %.2f\n", summary(s1)$r.squared, summary(s1)$sigma))
cat(sprintf("each extra visit is worth $%.2f (95%% CI: $%.2f to $%.2f)\n",
            coef(s1)[2], confint(s1)[2, 1], confint(s1)[2, 2]))

cat("\n--- multiple regression: + channel + region ---\n")
s2 <- lm(annual_spend ~ visits_per_year + channel + region, data = df)
print(round(summary(s2)$coefficients, 4))
cat(sprintf("\nR^2 = %.4f    Adjusted R^2 = %.4f    s_e = %.2f\n",
            summary(s2)$r.squared, summary(s2)$adj.r.squared, summary(s2)$sigma))
cat("reference levels: channel = ", levels(df$channel)[1],
    ", region = ", levels(df$region)[1], "\n", sep = "")

cat("\nmodel comparison (adjusted R^2 and AIC, never raw R^2):\n")
print(data.frame(model = c("visits only", "+ channel + region"),
                 adj_r2 = round(c(summary(s1)$adj.r.squared,
                                  summary(s2)$adj.r.squared), 4),
                 AIC = round(c(AIC(s1), AIC(s2)), 2)), row.names = FALSE)
cat("\nnested F-test:\n"); print(anova(s1, s2))

cat("\ndiagnostics on the full model:\n")
cat(sprintf("  Shapiro-Wilk on residuals   p = %.4f\n",
            shapiro.test(residuals(s2))$p.value))
cat(sprintf("  cor(|resid|, fitted)        %+.4f  (near 0 = equal variance)\n",
            cor(abs(residuals(s2)), fitted(s2))))
cat(sprintf("  influential points          %d\n",
            sum(cooks.distance(s2) > 4 / n)))

nd <- data.frame(visits_per_year = 12, channel = levels(df$channel)[1],
                 region = levels(df$region)[1])
cat("\nprediction for a 12-visit, ", levels(df$channel)[1], ", ",
    levels(df$region)[1], " customer:\n", sep = "")
print(round(predict(s2, nd, interval = "confidence"), 2))
print(round(predict(s2, nd, interval = "prediction"), 2))
cat("Use the PREDICTION interval when the question is about ONE customer.\n")

# ════════════════════════════════════════════════════════════════════════════
hr("PART 5 — EXECUTIVE SUMMARY")
# ════════════════════════════════════════════════════════════════════════════

cat(sprintf("\n1. The average customer spends $%.0f a year (95%% CI: $%.0f-$%.0f),\n",
            mean(spend), tt$conf.int[1], tt$conf.int[2]))
cat(sprintf("   with %.0f%% satisfied (95%% CI: %.0f%%-%.0f%%).\n",
            100 * ph, 100 * ci_p[1], 100 * ci_p[2]))
cat(sprintf("\n2. Mean spend %s the $250 target (one-sample t, p = %.4f).\n",
            ifelse(t1$p.value <= 0.05, "EXCEEDS", "does NOT clearly exceed"),
            t1$p.value))
cat(sprintf("\n3. Regional spending %s (ANOVA p = %.4f, eta-sq = %.3f).\n",
            ifelse(Fp <= 0.05, "DIFFERS significantly", "shows no significant difference"),
            Fp, ss[1] / sum(ss)))
cat(sprintf("\n4. Channel and satisfaction %s (chi-square p = %.4f, V = %.3f).\n",
            ifelse(ct$p.value <= 0.05, "ARE associated", "are NOT associated"),
            ct$p.value, V))
cat(sprintf("\n5. Visits, channel and region explain %.0f%% of spending variation;\n",
            100 * summary(s2)$adj.r.squared))
cat(sprintf("   a typical prediction is within about $%.0f.\n", summary(s2)$sigma))
cat("\nLIMITATION: this is OBSERVATIONAL data. Every finding is an ASSOCIATION.\n")
cat("Only a randomized experiment would license a causal claim.\n")

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

hist(spend, breaks = 25, col = "#8A5FBF", border = "white",
     main = "Annual spend", xlab = "$")
abline(v = mean(spend), col = "#0FA3A3", lwd = 2)
abline(v = median(spend), col = "#0B7A7A", lwd = 2, lty = 2)

boxplot(annual_spend ~ region, data = df, col = "#5B2A86",
        main = "Spend by region", ylab = "$", las = 2)

boxplot(annual_spend ~ channel, data = df, col = "#0FA3A3",
        main = "Spend by channel", ylab = "$", las = 2)

barplot(100 * prop.table(table(df$channel, df$satisfied), 1)[, "Yes"],
        col = "#0B7A7A", border = NA, las = 2, ylim = c(0, 100),
        main = "Satisfaction rate by channel", ylab = "% satisfied")

plot(df$visits_per_year, spend, pch = 19, col = "#5B2A86",
     main = "Spend vs visits", xlab = "visits per year", ylab = "$")
abline(s1, col = "#0FA3A3", lwd = 2)

plot(fitted(s2), residuals(s2), 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 capstone_plots.png\n")