################################
### HW2 Q5-Q8: hw2_braincancer2.csv (88 patients, time in months, status 1=death)
### Q5: log-rank test p-value, SRS vs SRT (stereo).
### ! Rscript "/mnt/c/Users/kirch/OneDrive/Documents/dev/ut-11-hdi/docs/part1/hw2_q5_8.R"
################################

suppressMessages(library(survival))
setwd("/mnt/c/Users/kirch/OneDrive/Documents/dev/ut-11-hdi/docs/part1")
bc = read.csv("hw2_braincancer2.csv")
cat("rows:", nrow(bc), " events:", sum(bc$status), "\n")
print(table(bc$stereo, bc$status, dnn = c("stereo", "status")))

cat("\n==== Q5: log-rank test, stereo (SRS vs SRT) ====\n")
lr = survdiff(Surv(time, status) ~ stereo, data = bc)
print(lr)
cat("chi-square statistic            :", lr$chisq, "\n")
cat("p-value = pchisq(chisq, df=1, lower=F):", pchisq(lr$chisq, df = 1, lower.tail = FALSE), "\n")
cat("Q5 ANSWER (2 dp)                :", round(lr$pvalue, 2), "\n")

cat("\n==== Q6: Kaplan-Meier S(24) in the SRT group ====\n")
km.srt = survfit(Surv(time, status) ~ 1, data = bc, subset = (stereo == "SRT"))
cat("SRT event times nearest to 24 months (summary rows around t=24):\n")
s = summary(km.srt)
idx = which(s$time > 18 & s$time < 30)
print(data.frame(time = round(s$time[idx], 3), n.risk = s$n.risk[idx],
                 n.event = s$n.event[idx], surv = round(s$surv[idx], 4)))
S24 = approx(km.srt$time, km.srt$surv, xout = 24, method = "constant")$y
cat("S_KM(24), SRT  [step function value at last event time <= 24]:", S24, "\n")
cat("summary(km.srt, times = 24) cross-check                      :",
    summary(km.srt, times = 24)$surv, "\n")
cat("Q6 ANSWER (2 dp)                                             :", round(S24, 2), "\n")

cat("\n==== Q7: Cox model ki + gtv + stereo; HR for SRT vs SRS ====\n")
bc$stereo = factor(bc$stereo)          # levels: SRS (reference), SRT
cat("stereo levels (first = reference):", levels(bc$stereo), "\n")
cox7 = coxph(Surv(time, status) ~ ki + gtv + stereo, data = bc)
print(summary(cox7)$coefficients)
b.srt = coef(cox7)["stereoSRT"]
cat("beta(stereoSRT)               :", b.srt, "\n")
cat("HR = exp(beta), SRT vs SRS    :", exp(b.srt), "\n")
cat("95% CI for HR                 :", exp(confint(cox7))["stereoSRT", ], "\n")
cat("Q7 ANSWER (2 dp)              :", round(exp(b.srt), 2), "\n")

cat("\n==== Q8: test H0: S_SRS(30) = S_SRT(30) using KM + Greenwood ====\n")
## same recipe as the Lecture Video 9 exercise: difference of KM estimates,
## variance = sum of the two Greenwood variances, z-test, two-sided p-value
km.g = survfit(Surv(time, status) ~ stereo, data = bc)
s30 = summary(km.g, times = 30)
print(data.frame(group = s30$strata, time = s30$time, n.risk = s30$n.risk,
                 surv = round(s30$surv, 4), std.err = round(s30$std.err, 4)))
S.srs = s30$surv[1];  se.srs = s30$std.err[1]
S.srt = s30$surv[2];  se.srt = s30$std.err[2]
diff = S.srt - S.srs
se.diff = sqrt(se.srs^2 + se.srt^2)
z = diff / se.diff
p = 2 * (1 - pnorm(abs(z)))
cat("S_SRT(30) - S_SRS(30)                :", diff, "\n")
cat("se of difference = sqrt(se1^2+se2^2) :", se.diff, "\n")
cat("95% CI for difference                :", diff - 1.96*se.diff, diff + 1.96*se.diff, "\n")
cat("z = diff / se                        :", z, "\n")
cat("two-sided p = 2*(1-pnorm(|z|))       :", p, "\n")
cat("Q8 ANSWER (2 dp)                     :", round(p, 2), "\n")