################################
### HW3 Q12-Q15: hw3_heart.csv (297 patients)
### Outcome hd = 1 if class in {1,2,3,4}, 0 if class == 0.
### Q12: AUC of trestbps alone for predicting hd (no CV).
### Needs: install.packages("pROC")   # for the cross-check only
### ! Rscript "/mnt/c/Users/kirch/OneDrive/Documents/dev/ut-11-hdi/docs/part2/hw3_q12_15.R"
################################

setwd("/mnt/c/Users/kirch/OneDrive/Documents/dev/ut-11-hdi/docs/part2")
heart = read.csv("hw3_heart.csv")
heart$hd = 1*(heart$class > 0)
cat("n =", nrow(heart), " | class table:"); print(table(heart$class))
cat("hd = 1 (disease):", sum(heart$hd), "  hd = 0:", sum(heart$hd == 0), "\n")

## nonparametric AUC = P(X > Y), ties count 1/2  (notes p. 41 Mann-Whitney form)
auc.hand = function(score, y) {
  x = score[y == 1]; yy = score[y == 0]
  mean(outer(x, yy, function(a, b) (a > b) + 0.5*(a == b)))
}

cat("\n==== Q12: AUC of trestbps (higher BP -> more likely disease) ====\n")
cat("cases mean trestbps   :", mean(heart$trestbps[heart$hd == 1]), "\n")
cat("controls mean trestbps:", mean(heart$trestbps[heart$hd == 0]), "\n")
auc12 = auc.hand(heart$trestbps, heart$hd)
cat("AUC by hand (Mann-Whitney, ties = 1/2):", auc12, "\n")
if (requireNamespace("pROC", quietly = TRUE)) {
  r = pROC::roc(heart$hd, heart$trestbps, direction = "<", quiet = TRUE)
  cat("pROC::roc(hd, trestbps, direction = '<')$auc:", as.numeric(r$auc), "\n")
}
cat("Q12 ANSWER (2 dp)                     :", round(auc12, 2), "\n")

cat("\n==== Q13: incremental value of trestbps, logistic models, training AUC ====\n")
m1 = glm(hd ~ age + thalach + chol,            family = "binomial", data = heart)
m2 = glm(hd ~ age + trestbps + chol + thalach, family = "binomial", data = heart)
cat("Model 2 coefficients:\n"); print(round(summary(m2)$coefficients, 4))
p1 = predict(m1, type = "response")
p2 = predict(m2, type = "response")
auc1 = auc.hand(p1, heart$hd)
auc2 = auc.hand(p2, heart$hd)
cat("AUC model 1 (age + thalach + chol)            :", auc1, "\n")
cat("AUC model 2 (age + trestbps + chol + thalach) :", auc2, "\n")
if (requireNamespace("pROC", quietly = TRUE))
  cat("pROC cross-check:", as.numeric(pROC::roc(heart$hd, p1, direction = "<", quiet = TRUE)$auc),
      as.numeric(pROC::roc(heart$hd, p2, direction = "<", quiet = TRUE)$auc), "\n")
iv = auc2 - auc1
cat("incremental value = AUC2 - AUC1               :", iv, "\n")
cat("Q13 ANSWER (2 dp)                             :", round(iv, 2), "\n")

cat("\n==== Q14: AIC(Model 1) - AIC(Model 2) ====\n")
cat("logLik m1:", as.numeric(logLik(m1)), " df =", attr(logLik(m1), "df"), "\n")
cat("logLik m2:", as.numeric(logLik(m2)), " df =", attr(logLik(m2), "df"), "\n")
cat("AIC by hand m1 = -2*ll + 2*4 :", -2*as.numeric(logLik(m1)) + 2*4, "\n")
cat("AIC by hand m2 = -2*ll + 2*5 :", -2*as.numeric(logLik(m2)) + 2*5, "\n")
cat("AIC(m1)                      :", AIC(m1), "\n")
cat("AIC(m2)                      :", AIC(m2), "\n")
d.aic = AIC(m1) - AIC(m2)
cat("AIC(m1) - AIC(m2)            :", d.aic, "\n")
cat("Q14 ANSWER (2 dp)            :", round(d.aic, 2), "\n")

cat("\n==== Q15: LOOCV Brier score, logistic model with all 13 predictors ====\n")
preds13 = c("age","sex","cp","trestbps","chol","fbs","restecg","thalach","exang",
            "oldpeak","slope","ca","thal")
f = as.formula(paste("hd ~", paste(preds13, collapse = " + ")))
cat("formula:", deparse(f), "\n")
n = nrow(heart)
brier = numeric(n)
for (i in 1:n) {                                   # Video 24 LOOCV loop
  m = suppressWarnings(glm(f, family = "binomial", data = heart[-i, ]))
  p = predict(m, newdata = heart[i, ], type = "response")
  brier[i] = (p - heart$hd[i])^2
}
cat("n fits                       :", n, "\n")
cat("LOOCV Brier = mean (p_i-A_i)^2:", mean(brier), "\n")
cat("training Brier, for contrast :",
    mean((predict(glm(f, family = "binomial", data = heart), type = "response") - heart$hd)^2), "\n")
cat("Q15 ANSWER (2 dp)            :", round(mean(brier), 2), "\n")

## sensitivity: coded categoricals (cp, restecg, slope, thal) as factors
h2 = heart; for (v in c("cp","restecg","slope","thal")) h2[[v]] = factor(h2[[v]])
b2 = numeric(n)
for (i in 1:n) {
  m = suppressWarnings(glm(f, family = "binomial", data = h2[-i, ]))
  b2[i] = (predict(m, newdata = h2[i, ], type = "response") - h2$hd[i])^2
}
cat("(sensitivity: cp/restecg/slope/thal as factors -> LOOCV Brier =", round(mean(b2), 4), ")\n")