Supplementary code to paper entatiled "The prognostic value of PI-RADS score in CyberKnife ultra-hypofractionated radiotherapy for localized prostate cancer." by Miszyk et. al. Code by Konrad Stawiski MD, PhD (konrad@konsta.com.pl, https://konsta.com.pl).

Analysis.R

# -*- coding: utf-8 -*-
options(reticulate.conda_binary = "/opt/conda/bin/conda") # conda path

library(survival)
library(ranger)
library(ggplot2)
library(dplyr)
library(readr)

completed = read_csv("pirads_v0_complete.csv")

completed$N_ognisk = as.factor(completed$N_ognisk)
completed$Gleason.Grade.Group = as.factor(completed$Gleason.Grade.Group)
completed$Risk.group = as.factor(completed$Risk.group)
completed$ht.przed.rt = as.factor(completed$ht.przed.rt)

# dokładane zmienne
completed$PSAdensity = completed$PSA.max / completed$objetosc
completed$suma_ratio = completed$V_suma / completed$objetosc
completed$Product = completed$wymiar_a * completed$wymiar_b

completed_before = completed
completed$dce = as.numeric(as.factor(completed$dce))
completed$lokalizacja_PIRADS = as.numeric(as.factor(completed$lokalizacja_PIRADS))
completed$ograniczenie = as.numeric(as.factor(completed$ograniczenie))
completed$strona = as.numeric(as.factor(completed$strona))
completed$TNM = as.numeric(as.factor(completed$TNM))

library(finalfit)
dependent_mfs  <- "Surv(MFStime, MFS)"
#dependent_dss <- "Surv(time, status_dss)"
#dependent_crr <- "Surv(time, status_crr)"
explanatory   <- colnames(completed)[5:ncol(completed)]



s = completed %>%
  finalfit(dependent_mfs, c(explanatory), digits = c(4,4,4,4))
s
data.table::fwrite(s, "coxph_unbalanced.csv")

#completed %>%
#  surv_plot(dependent_mfs, "wymiar_a", 
#            xlab="Time", pval=TRUE)
hist(completed$V_PiRads)
summary(completed$V_PiRads)
temp = dplyr::select(completed, -DFS, -DFStime)

library(party)
mod = ctree(Surv(MFStime, MFS) ~ V_PiRads,data = completed)
plot(mod)
table(completed$V_PiRads <= 352)
predict(mod, newdata = completed[1:10,], type="node")

library(readr)
pirads_v0_comp_balanced <- read_csv("pirads_v0_comp+balanced.csv")
mod = ctree(Surv(MFStime, MFS) ~ V_PiRads,data = pirads_v0_comp_balanced)
plot(mod)
table(completed$V_PiRads <= 352)
table(predict(mod, newdata = completed, type="node"))
predict(mod, newdata = completed, type="prob")
library(rms)
pred = as.factor(predict(mod, newdata = completed, type="node"))
model = cph(Surv(completed$MFStime, completed$MFS) ~ pred, x=T, y=T)
v = rms::validate(model, method = "boot", B=1000, dxy=T)
c_index = 0.5 * (v["Dxy","index.orig"] + 1)
c_index_val = 0.5 * (v["Dxy","index.corrected"] + 1)
reps = 500
dxy = numeric()  
for(i in 1 : reps) {
  try({ v = rms::validate(model, method = "boot", dxy=T)
  c_index_val = 0.5 * (v["Dxy","index.corrected"] + 1)
  dxy <- c(dxy,c_index_val) }) }
std_mean <- function(x) sd(x)/sqrt(length(x))
std_mean(dxy)
quantile(dxy, c(.025, .975))


# KOMBINACJE
pirads_v0_comp_balanced$PSAdensity = pirads_v0_comp_balanced$PSA.max / pirads_v0_comp_balanced$objetosc
pirads_v0_comp_balanced$suma_ratio = pirads_v0_comp_balanced$V_suma / pirads_v0_comp_balanced$objetosc
pirads_v0_comp_balanced$Product = pirads_v0_comp_balanced$wymiar_a * pirads_v0_comp_balanced$wymiar_b

explanatory   <- colnames(completed)[5:ncol(completed)]
explanatory = explanatory[-which(explanatory == "wymiar_x")]
explanatory = explanatory[-which(explanatory == "wymiar_y")]
explanatory = explanatory[-which(explanatory == "wymiar_z")]
explanatory = explanatory[-which(explanatory == "wymiar_b")]
explanatory = explanatory[-which(explanatory == "strona")]
explanatory = explanatory[-which(explanatory == "Risk.group")]
explanatory = explanatory[-which(explanatory == "ht.przed.rt")]
explanatory

kombinacje = c(combn(explanatory, 2, simplify = F), combn(explanatory, 3, simplify = F))
wyniki = list()

for(i in 1:length(kombinacje)) {
  print(i)
  try({ 
  temp = dplyr::select(pirads_v0_comp_balanced, MFS, MFStime, all_of(kombinacje[[i]]))
  mod = ctree(Surv(MFStime, MFS) ~ .,data = temp, controls = ctree_control(maxdepth = 2, testtype = "MonteCarlo")) 
  # plot(mod)
  predict(mod, newdata = completed, type="prob")
  
  pred = as.factor(predict(mod, newdata = completed, type="node"))
  model = cph(Surv(completed$MFStime, completed$MFS) ~ pred, x=T, y=T)
  v = rms::validate(model, method = "boot", B=100, dxy=T)
  c_index = 0.5 * (v["Dxy","index.orig"] + 1)
  c_index_val = 0.5 * (v["Dxy","index.corrected"] + 1)
  reps = 100
  dxy = numeric()  
  for(i in 1 : reps) {
    try({ v = rms::validate(model, method = "boot", B=100, dxy=T)
    c_index_val = 0.5 * (v["Dxy","index.corrected"] + 1)
    dxy <- c(dxy,c_index_val) }) }
  std_mean <- function(x) sd(x)/sqrt(length(x))
  std_mean(dxy)
  quantile(dxy, c(.025, .975))
  
  wyniki = c(wyniki, list(list(kombinacje[[i]], mod, pred, model, v, c_index, c_index_val, quantile(dxy, c(.025, .975)), dxy)))
  saveRDS(wyniki, "wyniki_finalne.RDS")
  })
}