library(survival)

source('../../Data/get_centenarian_data.R')
DF = get_centenarian_data()
	
## Construct a survival curve "by hand":
##
DF$Est_Pr_Event = DF$Deaths/DF$number_at_risk
DF$Est_Pr_No_Event = 1 - DF$Est_Pr_Event

KM_estimate = c(1, cumprod(DF$Est_Pr_No_Event))
DF$Pr_Event_Not_Happened_by_T = KM_estimate[-length(KM_estimate)]
DF$Hazard = DF$Est_Pr_Event / DF$Pr_Event_Not_Happened_by_T ## estimate the hazard function
print(DF)

par(mfrow=c(2, 1))

ages = DF$Age_ALB
##KM_estimate = log(KM_estimate+1.e-6) ## take the logarithm
ages = ages[1:(length(KM_estimate)-1)] + 1 - 100
KM_estimate = KM_estimate

sf = stepfun(ages, KM_estimate, right=FALSE)
plot(sf, pch=19, cex=0.75, xlab='additional years from 100', , ylab='survival fraction', main='survival curve')
grid()

Hazard = c(0, DF$Hazard)
sf = stepfun(ages, Hazard, right=FALSE)
plot(sf, ylim=c(0, 7000), pch=19, cex=0.75, xlab='additional years from 100', ylab='hazard function', main='hazard function')
grid()

par(mfrow=c(1, 1))