library(survival)
## Warning: package 'survival' was built under R version 4.2.3
library(survminer)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.2.3
## Loading required package: ggpubr
## Warning: package 'ggpubr' was built under R version 4.2.3
## 
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
## 
##     myeloma
library(knitr)
## Warning: package 'knitr' was built under R version 4.2.3
library(rmarkdown)
library(collett)
## 
## Attaching package: 'collett'
## The following object is masked from 'package:survminer':
## 
##     myeloma
## The following objects are masked from 'package:survival':
## 
##     bladder, kidney, lung, myeloma, ovarian
data(leukaemia)
head(leukaemia)
str(leukaemia)
## 'data.frame':    23 obs. of  8 variables:
##  $ patient  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ time     : int  1199 1111 530 1279 110 243 86 466 262 1850 ...
##  $ status   : int  0 0 0 1 1 1 1 1 1 0 ...
##  $ group    : int  1 1 1 1 1 1 1 1 1 2 ...
##  $ page     : int  24 19 17 17 28 37 17 15 29 37 ...
##  $ dage     : int  40 28 28 20 25 38 26 18 32 36 ...
##  $ precovery: int  1 1 1 1 1 1 0 1 1 1 ...
##  $ ptime    : chr  "29" "22" "34" "22" ...
#changing variables to numeric 
leukaemia$patient <- as.numeric(leukaemia$patient)
leukaemia$time <- as.numeric(leukaemia$time)
leukaemia$status <- as.numeric(leukaemia$status)
leukaemia$group <- as.numeric(leukaemia$group)
leukaemia$page <- as.numeric(leukaemia$page)
leukaemia$dage <- as.numeric(leukaemia$dage)
leukaemia$precovery <- as.numeric(leukaemia$precovery)
leukaemia$ptime <- as.numeric(leukaemia$ptime)
## Warning: NAs introduced by coercion
# Create survival object
surv_object <- Surv(time = leukaemia$time, event = leukaemia$status)

# Fit Kaplan-Meier model
km_fit <- survfit(surv_object ~ group, data = leukaemia)

# Plot Kaplan-Meier curves
plot(km_fit, xlab = "Time (days)", ylab = "Survival Probability", col = c("blue", "red"), main = "KAPLAN MEIER Of Leukaemia Patients")

#The Kaplan-Meier plot shows the survival probability over time for different groups . The plot provides a visual representation of how survival changes over time for each group.

#nelson aalen 
# Fit Nelson-Aalen model
na_fit <- survfit(Surv(time, status) ~ group, data = leukaemia, type = "fleming-harrington")
# Plot Nelson-Aalen estimates
ggsurvplot(na_fit,
           fun = "cumhaz",  # Specify cumulative hazard
           title = "Nelson-Aalen Cumulative Hazard by Treatment",
           xlab = "Time (days)",
           ylab = "Cumulative Hazard",
           
           ggtheme = theme_minimal()
           )

The plot shows the cumulative hazard over time for different groups. The cumulative hazard function increases over time, indicating the accumulation of risk.If the curves for different groups diverge, it suggests that the groups have different hazard rates over time. For example, if one group has a consistently higher cumulative hazard, it indicates a higher risk of the event in that group.

#log rank test

# Create survival object
surv_object <- Surv(time = leukaemia$time, event = leukaemia$status)

# Perform log-rank test
logrank_test <- survdiff(surv_object ~ leukaemia$group)
logrank_test
## Call:
## survdiff(formula = surv_object ~ leukaemia$group)
## 
##                   N Observed Expected (O-E)^2/E (O-E)^2/V
## leukaemia$group=1 9        6     4.34     0.632      1.01
## leukaemia$group=2 7        1     4.94     3.141      5.82
## leukaemia$group=3 7        5     2.72     1.915      2.56
## 
##  Chisq= 6.2  on 2 degrees of freedom, p= 0.04
# Fit Kaplan-Meier model
km_fit <- survfit(surv_object ~ leukaemia$group)

# Plot Kaplan-Meier curves with Log-Rank test p-value
ggsurvplot(km_fit, data = leukaemia, pval = TRUE, 
           )

           ggtheme = theme_minimal()

The plot shows the Kaplan-Meier survival curves for the different groups. The x-axis represents time, and the y-axis represents the survival probability. The p-value displayed on the plot indicates whether there is a statistically significant difference in survival between the groups. A small p-value suggests that the survival curves are significantly different

#cox proportional hazard model
# Create survival object
surv_object <- Surv(time = leukaemia$time, event = leukaemia$status)

# Fit the Cox proportional hazards model
cox_model <- coxph(surv_object ~ group, data = leukaemia)

# View the summary of the Cox model
summary(cox_model)
## Call:
## coxph(formula = surv_object ~ group, data = leukaemia)
## 
##   n= 23, number of events= 12 
## 
##          coef exp(coef) se(coef)     z Pr(>|z|)
## group 0.09108   1.09536  0.38073 0.239    0.811
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## group     1.095     0.9129    0.5194      2.31
## 
## Concordance= 0.528  (se = 0.088 )
## Likelihood ratio test= 0.06  on 1 df,   p=0.8
## Wald test            = 0.06  on 1 df,   p=0.8
## Score (logrank) test = 0.06  on 1 df,   p=0.8

The p-value of 0.811 indicates that this difference is statistically insignificant.