################################
### HW2 Q4
### Cox model on pharmacoSmoking: ttr ~ age + grp + employment + race + levelSmoking.
### The question supplies its own output table; check each statement's arithmetic
### from THAT table. A refit on the packaged asaur data confirms the reference
### levels (its coefficients differ slightly - different data version).
### ! Rscript "/mnt/c/Users/kirch/OneDrive/Documents/dev/ut-11-hdi/docs/part1/hw2_q4.R"
################################

suppressMessages({library(survival); library(asaur)})

## beta values exactly as printed in the question
b = c(age = -0.03134, grppatchOnly = 0.65163, employmentother = 0.81910,
      employmentpt = 0.59316, racehispanic = -0.43145, raceother = -0.39418,
      racewhite = -0.31526, levelSmokinglight = 0.15678)
p = c(0.000151, 0.000150, 0.009226, 0.030827, 0.450410, 0.652757, 0.177213, 0.552406)

cat("==== Table from the question: beta -> exp(beta) ====\n")
print(data.frame(beta = b, exp.beta = round(exp(b), 5), p = p))

## reference levels: R uses the first factor level
m = coxph(Surv(ttr, relapse) ~ age + grp + employment + race + levelSmoking,
          data = pharmacoSmoking)
cat("\nreference levels (from refit on asaur::pharmacoSmoking):\n")
for (v in c("grp", "employment", "race", "levelSmoking"))
  cat(sprintf("  %-13s ref = %s\n", v, levels(pharmacoSmoking[[v]])[1]))
cat("(refit coefficient names match the table:",
    all(names(coef(m)) == names(b)), ")\n")

cat("\n==== Statement checks ====\n")
cat("[A] patch-only vs combination: HR = exp(0.65163) =", round(exp(b["grppatchOnly"]), 3),
    "-> hazard is 1.92x HIGHER, not 'decreased'; and Cox output is about hazard, not probability\n")
cat("[B] signs of the significant betas (p<0.05): age", sign(b["age"]),
    "| patchOnly", sign(b["grppatchOnly"]), "| empl other", sign(b["employmentother"]),
    "| empl pt", sign(b["employmentpt"]), "-> only age DEcreases the hazard\n")
cat("[C] light vs heavy smoking (ref = heavy): HR = exp(0.15678) =",
    round(exp(b["levelSmokinglight"]), 3), "\n")
cat("[D] 'other' vs 'pt' employment: HR = exp(0.81910 - 0.59316) =",
    round(exp(b["employmentother"] - b["employmentpt"]), 3),
    " (2.27 is 'other' vs FULL TIME, the reference)\n")

cat("\nANSWER: [C] holding the rest constant, light smokers' hazard of relapse is",
    round(exp(b["levelSmokinglight"]), 2), "times that of heavy smokers\n")