Harrell’s C-index

### Load library, functions 
source("C:/local_work_ming/workspaces/r_workspace/R_customized_global_functions/my_functional_scripts/my_global_functions.R")
source("C:/local_work_ming/workspaces/r_workspace/R_customized_global_functions/my_functional_scripts/my_libraries.R")
source("C:/local_work_ming/workspaces/r_workspace/R_customized_global_functions/my_functional_scripts/R_intergroup_diff_test_functions.R")
# Load the survival package
library(survival)
library(rms)

# Datensatz laden
data.path <- "C:/local_work_ming/workspaces/r_workspace/florian_scurt_data/data"

——- Old Risk Score + Univariate ———————————————–

old.uni.data <- read.xlsx(paste(data.path, "1 Masterdatei_Vasculitis Modell Alt Univariat.xlsx", sep = "/"), sheetIndex = 1)
names(old.uni.data) <- c("ID", "TimeToDialysis", "DialysisStatus", "OldRiskGroup")
old.uni.data$OldRiskGroup <- factor(old.uni.data$OldRiskGroup, levels = seq(1, 3, 1))
table(old.uni.data$OldRiskGroup)

 1  2  3 
50 70 44 

# Fit a Cox proportional hazards model
# Here, we use age and sex as predictors for the survival time of lung cancer patients
old_uni_cox_model <- coxph(Surv(TimeToDialysis, DialysisStatus) ~ OldRiskGroup, data = old.uni.data, id=ID)
print("---- Cox proportional hazards model -----------")

[1] "---- Cox proportional hazards model -----------"

print(old_uni_cox_model)

Call:
coxph(formula = Surv(TimeToDialysis, DialysisStatus) ~ OldRiskGroup, 
    data = old.uni.data, id = ID)

                 coef exp(coef) se(coef)     z        p
OldRiskGroup2  1.9683    7.1584   0.4754 4.141 3.47e-05
OldRiskGroup3  2.7448   15.5615   0.4900 5.602 2.12e-08

Likelihood ratio test=51.36  on 2 df, p=7.044e-12
n= 164, number of events= 75 

# Calculate Harrell's C-index using the concordance index function (concordance)
# The concordance function is part of the survival package and directly calculates the C-index
old_uni_c_index <- concordance(old_uni_cox_model)$concordance

# Print the C-index
print(paste("Old Risk Score + Univariate Harrell's C-index:", old_uni_c_index))

[1] "Old Risk Score + Univariate Harrell's C-index: 0.738438017581858"

——- New Risk Score + Univariate ———————————————–

new.uni.data <- read.xlsx(paste(data.path, "2 Masterdatei_Vasculitis Modell Neu Univariat.xlsx", sep = "/"), sheetIndex = 1)
names(new.uni.data) <- c("ID", "TimeToDialysis", "DialysisStatus", "NewRiskGroup")
new.uni.data$NewRiskGroup <- factor(new.uni.data$NewRiskGroup, levels = seq(1, 4, 1))
table(new.uni.data$NewRiskGroup)

 1  2  3  4 
51 55 46 12 

#head(new.uni.data)

# Fit a Cox proportional hazards model
# Here, we use age and sex as predictors for the survival time of lung cancer patients
new_uni_cox_model <- coxph(Surv(TimeToDialysis, DialysisStatus) ~ NewRiskGroup, data = new.uni.data, id=ID)
print("---- Cox proportional hazards model -----------")

[1] "---- Cox proportional hazards model -----------"

print(new_uni_cox_model)

Call:
coxph(formula = Surv(TimeToDialysis, DialysisStatus) ~ NewRiskGroup, 
    data = new.uni.data, id = ID)

                 coef exp(coef) se(coef)     z        p
NewRiskGroup2  1.0245    2.7858   0.4112 2.492   0.0127
NewRiskGroup3  2.1035    8.1951   0.3963 5.308 1.11e-07
NewRiskGroup4  2.4922   12.0882   0.4778 5.216 1.82e-07

Likelihood ratio test=49.25  on 3 df, p=1.153e-10
n= 164, number of events= 75 

# Calculate Harrell's C-index using the concordance index function (concordance)
# The concordance function is part of the survival package and directly calculates the C-index
new_uni_c_index <- concordance(new_uni_cox_model)$concordance

# Print the C-index
print(paste("New Risk Score + Univariate Harrell's C-index:", new_uni_c_index))

[1] "New Risk Score + Univariate Harrell's C-index: 0.753471779844566"

——- Old Risk Score + Multi-Variable ———————————————–

old.multi.data <- read.xlsx(paste(data.path, "3 Masterdatei_Vasculitis Modell Alt Multivariat.xlsx", sep = "/"), sheetIndex = 1)
names(old.multi.data) <- c("ID", "TimeToDialysis", "DialysisStatus", "OldRiskGroup", "CVD", "Anemia", "NephroticSyndrome", "Immunology_C3", "Blood_CrP", "IT_Plasmapheresis")
old.multi.data$OldRiskGroup <- factor(old.multi.data$OldRiskGroup, levels = seq(1, 3, 1))
table(old.multi.data$OldRiskGroup)

 1  2  3 
50 70 44 

old.multi.data$CVD <- factor(old.multi.data$CVD)
table(old.multi.data$CVD)

  0   1 
110  47 

old.multi.data$Anemia <- factor(old.multi.data$Anemia)
table(old.multi.data$Anemia)

 0  1 
80 76 

old.multi.data$NephroticSyndrome <- factor(old.multi.data$NephroticSyndrome)
table(old.multi.data$NephroticSyndrome)

  0   1 
120  44 

old.multi.data$IT_Plasmapheresis <- factor(old.multi.data$IT_Plasmapheresis)
table(old.multi.data$IT_Plasmapheresis)

  0   1 
121  43 

# Fit a Cox proportional hazards model
# Here, we use age and sex as predictors for the survival time of lung cancer patients
old_multi_cox_model <- coxph(Surv(TimeToDialysis, DialysisStatus) ~ OldRiskGroup+CVD+Anemia+NephroticSyndrome+Immunology_C3+Blood_CrP+IT_Plasmapheresis, data = old.multi.data, id=ID)
print("---- Cox proportional hazards model -----------")

[1] "---- Cox proportional hazards model -----------"

print(old_multi_cox_model)

Call:
coxph(formula = Surv(TimeToDialysis, DialysisStatus) ~ OldRiskGroup + 
    CVD + Anemia + NephroticSyndrome + Immunology_C3 + Blood_CrP + 
    IT_Plasmapheresis, data = old.multi.data, id = ID)

                        coef exp(coef)  se(coef)      z        p
OldRiskGroup2       2.061751  7.859720  0.540262  3.816 0.000136
OldRiskGroup3       2.489304 12.052885  0.558080  4.460 8.18e-06
CVD1                0.354267  1.425136  0.280692  1.262 0.206905
Anemia1             0.288055  1.333830  0.279548  1.030 0.302808
NephroticSyndrome1  0.939968  2.559900  0.267119  3.519 0.000433
Immunology_C3      -1.382202  0.251025  0.469297 -2.945 0.003227
Blood_CrP           0.004115  1.004123  0.001652  2.490 0.012768
IT_Plasmapheresis1  0.797814  2.220682  0.259226  3.078 0.002086

Likelihood ratio test=94.33  on 8 df, p=< 2.2e-16
n= 155, number of events= 71 
   (9 Beobachtungen als fehlend gelöscht)

# Calculate Harrell's C-index using the concordance index function (concordance)
# The concordance function is part of the survival package and directly calculates the C-index
old_multi_c_index <- concordance(old_multi_cox_model)$concordance

# Print the C-index
print(paste("Old Risk Score + Multi-Variable Harrell's C-index:", old_multi_c_index))

[1] "Old Risk Score + Multi-Variable Harrell's C-index: 0.838056680161943"

——- New Risk Score + Multi-Variable ———————————————–

new.multi.data <- read.xlsx(paste(data.path, "4 Masterdatei_Vasculitis Modell Neu Multivariat.xlsx", sep = "/"), sheetIndex = 1)
names(new.multi.data) <- c("ID", "TimeToDialysis", "DialysisStatus", "NewRiskGroup", "CVD", "Anemia", "NephroticSyndrome", "Immunology_C3", "Blood_CrP", "IT_Plasmapheresis")
new.multi.data$NewRiskGroup <- factor(new.multi.data$NewRiskGroup, levels = seq(1, 4, 1))
table(new.multi.data$NewRiskGroup)

 1  2  3  4 
51 55 46 12 

new.multi.data$CVD <- factor(new.multi.data$CVD)
table(new.multi.data$CVD)

  0   1 
110  47 

new.multi.data$Anemia <- factor(new.multi.data$Anemia)
table(new.multi.data$Anemia)

 0  1 
80 76 

new.multi.data$NephroticSyndrome <- factor(new.multi.data$NephroticSyndrome)
table(new.multi.data$NephroticSyndrome)

  0   1 
120  44 

new.multi.data$IT_Plasmapheresis <- factor(new.multi.data$IT_Plasmapheresis)
table(new.multi.data$IT_Plasmapheresis)

  0   1 
121  43 

#head(new.multi.data)

# Fit a Cox proportional hazards model
# Here, we use age and sex as predictors for the survival time of lung cancer patients
new_multi_cox_model <- coxph(Surv(TimeToDialysis, DialysisStatus) ~ NewRiskGroup+CVD+Anemia+NephroticSyndrome+Immunology_C3+Blood_CrP+IT_Plasmapheresis, data = new.multi.data, id=ID)
print("---- Cox proportional hazards model -----------")

[1] "---- Cox proportional hazards model -----------"

print(new_multi_cox_model)

Call:
coxph(formula = Surv(TimeToDialysis, DialysisStatus) ~ NewRiskGroup + 
    CVD + Anemia + NephroticSyndrome + Immunology_C3 + Blood_CrP + 
    IT_Plasmapheresis, data = new.multi.data, id = ID)

                        coef exp(coef)  se(coef)      z        p
NewRiskGroup2       0.674518  1.963087  0.418840  1.610 0.107301
NewRiskGroup3       1.530595  4.620924  0.417158  3.669 0.000243
NewRiskGroup4       2.150905  8.592629  0.528344  4.071 4.68e-05
CVD1                0.471371  1.602189  0.277868  1.696 0.089813
Anemia1             0.363665  1.438592  0.259338  1.402 0.160831
NephroticSyndrome1  0.888129  2.430578  0.269117  3.300 0.000966
Immunology_C3      -1.528116  0.216944  0.482032 -3.170 0.001524
Blood_CrP           0.002902  1.002906  0.001601  1.813 0.069888
IT_Plasmapheresis1  0.671105  1.956398  0.261650  2.565 0.010321

Likelihood ratio test=86.47  on 9 df, p=8.311e-15
n= 155, number of events= 71 
   (9 Beobachtungen als fehlend gelöscht)

# Calculate Harrell's C-index using the concordance index function (concordance)
# The concordance function is part of the survival package and directly calculates the C-index
new_multi_c_index <- concordance(new_multi_cox_model)$concordance

# Print the C-index
print(paste("New Risk Score + Multi-Variable Harrell's C-index:", new_multi_c_index))

[1] "New Risk Score + Multi-Variable Harrell's C-index: 0.833728884545581"

Comparison

print(paste("Old Risk Score + Univariate Harrell's C-index:", old_uni_c_index))

[1] "Old Risk Score + Univariate Harrell's C-index: 0.738438017581858"

print(paste("New Risk Score + Univariate Harrell's C-index:", new_uni_c_index))

[1] "New Risk Score + Univariate Harrell's C-index: 0.753471779844566"

print(paste("Old Risk Score + Multi-Variable Harrell's C-index:", old_multi_c_index))

[1] "Old Risk Score + Multi-Variable Harrell's C-index: 0.838056680161943"

print(paste("New Risk Score + Multi-Variable Harrell's C-index:", new_multi_c_index))

[1] "New Risk Score + Multi-Variable Harrell's C-index: 0.833728884545581"