# ============================================================================
#  04-03  Descriptive Statistics Dashboard
#  Chapters 03-01, 03-02, 04-01, 04-02
#
#  Run from the project folder:   Rscript r/analysis.R
# ============================================================================

df <- read.csv("data/employees.csv", stringsAsFactors = TRUE)
sal <- df$salary
n <- length(sal)

stat_mode <- function(x) {
  tb <- table(x)
  as.numeric(names(tb)[tb == max(tb)])
}

describe_one <- function(x) {
  q <- quantile(x, c(0.25, 0.75))
  c(n        = length(x),
    mean     = mean(x),
    median   = median(x),
    sd       = sd(x),
    var      = var(x),
    cv_pct   = 100 * sd(x) / mean(x),
    min      = min(x),
    Q1       = q[[1]],
    Q3       = q[[2]],
    max      = max(x),
    IQR      = q[[2]] - q[[1]],
    skew_SK  = 3 * (mean(x) - median(x)) / sd(x))
}

# ── 1. OVERALL SUMMARY ──────────────────────────────────────────────────────
cat("\n===== 1. SALARY — OVERALL =====\n")
print(round(describe_one(sal), 3))
cat("mode(s):", paste(stat_mode(sal), collapse = ", "), "\n")
cat("trimmed mean (10% each end):", round(mean(sal, trim = 0.1), 2), "\n")

cat("\nPopulation vs sample spread (for contrast):\n")
cat(sprintf("  sample sd     %.2f   (divides by n-1)\n", sd(sal)))
cat(sprintf("  population sd %.2f   (divides by N)\n", sqrt(var(sal) * (n - 1) / n)))

# ── 2. BY DEPARTMENT ────────────────────────────────────────────────────────
cat("\n===== 2. SALARY BY DEPARTMENT =====\n")
by_dept <- do.call(rbind, lapply(split(sal, df$department), describe_one))
print(round(by_dept, 2))

cat("\nMost / least variable department (by CV, not by sd):\n")
cv <- by_dept[, "cv_pct"]
cat("  most  ", names(which.max(cv)), sprintf("(%.1f%%)\n", max(cv)))
cat("  least ", names(which.min(cv)), sprintf("(%.1f%%)\n", min(cv)))

# ── 3. THE WEIGHTED-MEAN CHECK ──────────────────────────────────────────────
cat("\n===== 3. WEIGHTED MEAN CHECK =====\n")
means  <- by_dept[, "mean"]
counts <- by_dept[, "n"]

cat(sprintf("overall mean          %.2f\n", mean(sal)))
cat(sprintf("weighted mean of means %.2f   <- must match\n",
            weighted.mean(means, counts)))
cat(sprintf("PLAIN average of means %.2f   <- only equal when every n is equal\n",
            mean(means)))

# ── 4. OUTLIERS, BOTH RULES ─────────────────────────────────────────────────
cat("\n===== 4. OUTLIERS =====\n")

Q1 <- quantile(sal, 0.25); Q3 <- quantile(sal, 0.75); iqr <- Q3 - Q1
lower <- Q1 - 1.5 * iqr;   upper <- Q3 + 1.5 * iqr
extreme_lo <- Q1 - 3 * iqr; extreme_hi <- Q3 + 3 * iqr

df$z <- as.numeric(scale(sal))
df$iqr_flag <- sal < lower | sal > upper
df$z_flag   <- abs(df$z) > 3

cat(sprintf("IQR fences:    %.0f  to  %.0f\n", lower, upper))
cat(sprintf("Extreme fences: %.0f  to  %.0f\n", extreme_lo, extreme_hi))
cat("IQR-rule outliers:", sum(df$iqr_flag), "\n")
cat("z-rule  outliers:", sum(df$z_flag),  "\n")

flagged <- df[df$iqr_flag | df$z_flag,
              c("employee_id", "department", "salary", "z", "iqr_flag", "z_flag")]
if (nrow(flagged)) {
  print(flagged[order(-flagged$salary), ], row.names = FALSE, digits = 4)
} else {
  cat("(none)\n")
}

if (sum(df$iqr_flag) != sum(df$z_flag)) {
  cat("\nThe two rules DISAGREE. The extreme value inflates s, shrinking its own\n",
      "z-score -- the masking effect. The resistant IQR rule is the one to trust.\n")
}

cat("\nEffect of the extreme value on each statistic:\n")
keep <- !df$iqr_flag
comp <- rbind(with_outlier    = describe_one(sal),
              without_outlier = describe_one(sal[keep]))
print(round(comp[, c("mean", "median", "sd", "IQR")], 2))

# ── 5. MOST UNUSUAL EMPLOYEES ───────────────────────────────────────────────
cat("\n===== 5. FIVE MOST UNUSUAL SALARIES (by |z|) =====\n")
top <- df[order(-abs(df$z))[1:5], c("employee_id", "department", "salary", "z")]
print(top, row.names = FALSE, digits = 4)

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

hist(sal, breaks = 15, col = "#8A5FBF", border = "white",
     main = "Salary distribution", xlab = "Salary")
abline(v = mean(sal),   col = "#0FA3A3", lwd = 2)
abline(v = median(sal), col = "#0B7A7A", lwd = 2, lty = 2)
legend("topright", c("mean", "median"), col = c("#0FA3A3", "#0B7A7A"),
       lwd = 2, lty = c(1, 2), bty = "n")

boxplot(salary ~ department, data = df, col = "#0FA3A3",
        main = "Salary by department", xlab = "", ylab = "Salary", las = 2)

barplot(sort(by_dept[, "mean"]), col = "#5B2A86", border = NA, las = 2,
        main = "Mean salary by department", ylab = "Salary")

plot(df$years_service, sal, pch = 19, col = "#5B2A86",
     main = "Salary vs years of service",
     xlab = "Years of service", ylab = "Salary")
abline(lm(sal ~ df$years_service), col = "#0FA3A3", lwd = 2)

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