survival_script.R

library(survival)
library(coxme)
library(dplyr)
library(survminer)

# Read the data
data <- read.csv('out/ACCplatoon.csv')


# Drop rows with missing values
dfdata <- na.omit(data)

# Define the model formula
model <- coxme(Surv(duration, died) ~ mean_speed + length + platoon_compo + ACC_in_platoon + (1 | platoon_compo), data = dfdata)


# Print the model summary
summary(model)
# Fit the null model (model with no covariates)
null_formula <- Surv(duration, died) ~ 1
null_model <- coxph(null_formula, data = data)

# Extract the log-likelihoods
null_log_likelihood <- logLik(null_model)
final_log_likelihood <- logLik(model)

# Print the log-likelihoods
cat("Null Model Log-Likelihood:", null_log_likelihood, "\n")
cat("Final Model Log-Likelihood:", final_log_likelihood, "\n")

# Calculate rho-squared
rho_squared <- 1 - exp((final_log_likelihood - null_log_likelihood) / -length(data$duration))
cat("Rho-squared (ρ²):", rho_squared, "\n")
# Extract the frailty terms
frailty_terms <- ranef(model)
frailty_terms_numeric <- as.numeric(unlist(frailty_terms))


# Perform Kolmogorov-Smirnov test for log-normal distribution
ks_test <- ks.test(frailty_terms_numeric, "pnorm", mean = mean(frailty_terms_numeric), sd = sd(frailty_terms_numeric)/2)
print(ks_test)


data$ACC_in_platoon <- as.factor(data$ACC_in_platoon)

# Estimation des courbes de Kaplan-Meier par groupe
km_fit <- survfit(Surv(duration, died) ~ ACC_in_platoon, data = data)

# Visualisation avec survminer
ggsurvplot(km_fit,
           data = data,
           conf.int = TRUE,
           risk.table = TRUE,
           pval = TRUE,
           title = "Survival comparison between platoons with and without ACCs",
           xlab = "Time",
           ylab = "Survival probability",
           palette = c( "#2E9FDF","#E7B800"),
           legend.title = "ACC in platoon",
           legend.labs = c("No", "Yes"))