# Install required packages
install.packages(c("readxl", "survival", "survminer", "cmprsk"))

# Load packages
library(readxl)
library(survival)
library(survminer)
library(cmprsk)

# Import the dataset from Excel
dati <- read_xlsx("survival_exercise_data.xlsx", sheet = "survival_data")

# Inspect the first rows
head(dati)

# Check variable structure
str(dati)

# Convert categorical variables to factors
dati$sex <- factor(dati$sex)
dati$treatment <- factor(dati$treatment)
dati$age_group <- factor(dati$age_group)

############################################################
# PART 1: CLASSICAL SURVIVAL ANALYSIS
############################################################

# Create the survival object
# days = follow-up time
# status = 1 indicates the event occurred
# status = 0 indicates censoring
surv_obj <- with(dati, Surv(days, status == 1))

# Kaplan-Meier curve by treatment group
km_treatment <- survfit(surv_obj ~ treatment, data = dati)

ggsurvplot(
  km_treatment,
  data = dati,
  risk.table = TRUE,
  pval = TRUE,
  conf.int = TRUE,
  xlab = "Time in days",
  ylab = "Survival probability",
  legend.title = "Treatment group"
)

# Kaplan-Meier curve by age group
km_age <- survfit(surv_obj ~ age_group, data = dati)


ggsurvplot(
  km_age,
  data = dati,
  risk.table = TRUE,
  pval = TRUE,
  conf.int = TRUE,
  xlab = "Time in days",
  ylab = "Survival probability",
  legend.title = "Age group"
)

# Log-rank test by treatment group
logrank_treatment <- survdiff(surv_obj ~ treatment, data = dati)
logrank_treatment

# Log-rank test by age group
logrank_age <- survdiff(surv_obj ~ age_group, data = dati)
logrank_age

############################################################
# COX PROPORTIONAL HAZARDS MODEL
############################################################

# Fit a Cox proportional hazards model
# The model includes treatment, age, and sex
cox_model <- coxph(
  Surv(days, status == 1) ~ treatment + age + sex,
  data = dati
)

# Display model results
summary(cox_model)

# Extract hazard ratios
exp(coef(cox_model))

# Extract 95% confidence intervals for hazard ratios
exp(confint(cox_model))

############################################################
# CHECK PROPORTIONAL HAZARDS ASSUMPTION
############################################################

# Test the proportional hazards assumption
# Based on Schoenfeld residuals
ph_test <- cox.zph(cox_model)
ph_test

# Plot Schoenfeld residuals
plot(ph_test)

############################################################
# PART 2: COMPETING RISKS ANALYSIS
############################################################

# In the competing risks setting:
# cr_status = 0 means censored
# cr_status = 1 means event of interest
# cr_status = 2 means competing event

# Check distribution of competing risk status
table(dati$cr_status)

# Cumulative incidence function by treatment group
cif_treatment <- cuminc(
  ftime = dati$cr_days,
  fstatus = dati$cr_status,
  group = dati$treatment
)

# Display cumulative incidence results
cif_treatment

# Plot cumulative incidence by treatment group
plot(
  cif_treatment,
  xlab = "Time in days",
  ylab = "Cumulative incidence",
  lty = 1,
  col = 1:4
)

# Cumulative incidence function by age group
cif_age <- cuminc(
  ftime = dati$cr_days,
  fstatus = dati$cr_status,
  group = dati$age_group
)

# Display cumulative incidence results
cif_age

# Plot cumulative incidence by age group
plot(
  cif_age,
  xlab = "Time in days",
  ylab = "Cumulative incidence",
  lty = 1,
  col = 1:4
)

############################################################
# OPTIONAL: FINE-GRAY MODEL
############################################################

# Create numeric covariates for the Fine-Gray model
dati$treat_num <- ifelse(dati$treatment == "New treatment", 1, 0)
dati$sex_num <- ifelse(dati$sex == "Male", 1, 0)

# Fit Fine-Gray model for the event of interest
# failcode = 1 indicates the event of interest
# cencode = 0 indicates censoring
fg_model <- crr(
  ftime = dati$cr_days,
  fstatus = dati$cr_status,
  cov1 = dati[, c("treat_num", "age", "sex_num")],
  failcode = 1,
  cencode = 0
)

# Display Fine-Gray model results
summary(fg_model)