################################
### HW1 Q12-Q14
### x = (190,47,5,84,89,103,311), delta = (1,1,0,1,0,1,0)
### Q12: Kaplan-Meier S(100).  Q13: Nelson-Aalen Lambda(85).
### Q14: Fleming-Harrington S(190) = exp(-Lambda_NA(190)).
### ! Rscript "/mnt/c/Users/kirch/OneDrive/Documents/dev/ut-11-hdi/docs/part1/q12_14.R"
################################

suppressMessages(library(survival))

x     = c(190, 47, 5, 84, 89, 103, 311)
delta = c(1,   1,  0, 1,  0,  1,   0)

## Hand computation: event times with risk sets
## t_k :  47   84   103   190
## n_k :   6    5     3     2
## d_k :   1    1     1     1

cat("==== Q12: Kaplan-Meier S(100) ====\n")
S.km.100 = (1 - 1/6) * (1 - 1/5)
cat("hand: (5/6)*(4/5)                :", S.km.100, "\n")
km = survfit(Surv(x, delta) ~ 1)
cat("survfit KM at t=100              :",
    approx(km$time, km$surv, xout = 100, method = "constant")$y, "\n")
cat("Q12 ANSWER (2 dp)                :", round(S.km.100, 2), "\n\n")

cat("full KM table from survfit:\n")
print(summary(km))

cat("\n==== Q13: Nelson-Aalen Lambda(85) ====\n")
Lambda.na.85 = 1/6 + 1/5
cat("hand: 1/6 + 1/5                  :", Lambda.na.85, "\n")
na = survfit(Surv(x, delta) ~ 1, stype = 2, ctype = 1)
cat("survfit cumhaz at t=85           :",
    approx(na$time, na$cumhaz, xout = 85, method = "constant")$y, "\n")
cat("Q13 ANSWER (2 dp)                :", round(Lambda.na.85, 2), "\n")

cat("\nNelson-Aalen cumulative hazard by time:\n")
print(data.frame(time = na$time, cumhaz = na$cumhaz))

cat("\n==== Q14: Fleming-Harrington S(190) ====\n")
Lambda.na.190 = 1/6 + 1/5 + 1/3 + 1/2
cat("hand: Lambda_NA(190) = 1/6+1/5+1/3+1/2 :", Lambda.na.190, "\n")
S.fh.190 = exp(-Lambda.na.190)
cat("hand: S_FH(190) = exp(-1.2)            :", S.fh.190, "\n")
fh = survfit(Surv(x, delta) ~ 1, type = "fh")
cat("survfit type='fh' at t=190             :",
    approx(fh$time, fh$surv, xout = 190, method = "constant")$y, "\n")
cat("Q14 ANSWER (2 dp)                      :", round(S.fh.190, 2), "\n")