****if deces: 0=death; 1 alive)*************
#Prepare: packages, data
library(foreign, lib.loc = "C:/Program Files/R/R-3.6.1/library")
library(survival)
#read data and group for variable
data=read.csv("C:/Users/BINH THANG-TRAN/Dropbox/PhD/Data/BS TRi Ca/Data/data.csv")
head(data)
## ID dossier chirurgien age sexe tabagisme Apparition ASIA_preop Frankel_preop
## 1 1 928667 ZW 55 1 1 2 D D
## 2 2 5473345 ZW 71 1 0 2 E E
## 3 3 5466545 GB 68 2 1 2 D D
## 4 4 553346 DS 71 1 1 1 D D
## 5 5 734329 DS 63 1 1 1 E E
## 6 6 390270 GB 65 2 1 1 D D
## ambulation_preop fctSphincter approcheChx dateChx ambulationPO
## 1 2 2 1 10/6/2018 2
## 2 1 2 1 1/7/2018 1
## 3 2 1 1 1/3/2018 2
## 4 2 2 1 16/12/2017 2
## 5 1 2 1 16/12/2017 1
## 6 2 2 2 28/07/2016 1
## amelioration_ambulance ASIA_PO FrankelPO duree_sejour ATCD_RoRx RoRx_PO
## 1 0 D D 13 0 1
## 2 0 E E 4 1 0
## 3 0 D D 12 0 1
## 4 0 D D 25 1 1
## 5 0 E E 5 1 1
## 6 1 E D 16 1 0
## TypeSurgery chimiotxPO Tokuhashi Tokuhashi_cat SINS SINS_cat
## 1 CO 1 5 1 12 2
## 2 CO 0 9 2 8 2
## 3 CO 1 6 1 10 2
## 4 LA 0 6 1 6 1
## 5 CO 0 6 1 11 2
## 6 CO 1 5 1 10 2
## ImprovementPainPO deces date_deces SurviePosOpMOIS LastFollowUp BloodLoss
## 1 1 0 31/10/2018 4.77 31/10/2018 1500
## 2 1 0 31/10/2018 4.07 15/08/2018 500
## 3 1 0 31/10/2018 8.13 1/4/2018 300
## 4 1 1 12/2/2018 1.93 9/2/2018 200
## 5 1 1 1/2/2018 1.57 30/01/2018 2000
## 6 1 1 12/8/2017 12.67 18/07/2017 200
## DurationSurgery Histology
## 1 299 1
## 2 124 1
## 3 108 1
## 4 110 1
## 5 126 1
## 6 120 1
names(data) #name of variables
## [1] "ID" "dossier" "chirurgien"
## [4] "age" "sexe" "tabagisme"
## [7] "Apparition" "ASIA_preop" "Frankel_preop"
## [10] "ambulation_preop" "fctSphincter" "approcheChx"
## [13] "dateChx" "ambulationPO" "amelioration_ambulance"
## [16] "ASIA_PO" "FrankelPO" "duree_sejour"
## [19] "ATCD_RoRx" "RoRx_PO" "TypeSurgery"
## [22] "chimiotxPO" "Tokuhashi" "Tokuhashi_cat"
## [25] "SINS" "SINS_cat" "ImprovementPainPO"
## [28] "deces" "date_deces" "SurviePosOpMOIS"
## [31] "LastFollowUp" "BloodLoss" "DurationSurgery"
## [34] "Histology"
##explaination
###subset dataset -Independent vars:
data1 <- subset(data, select=c(ID, sexe, age, tabagisme,Histology, ASIA_preop, ASIA_PO, Tokuhashi_cat,ambulation_preop, ambulationPO, RoRx_PO,chimiotxPO, SurviePosOpMOIS, deces))
Independent variables
data1$age_gr <- ifelse(data1$age >= 60,c("0"), c("1")) #Age: 60+; <60; (age) -need to group
data1$sexe = factor(data1$sexe) # Sex: male - female (sexe)
data1$tabagisme = factor(data1$tabagisme) #Tobacco use: Yes - No (tabagisme)
data1$Histology = factor(data1$Histology) #Histologic type: ADC; Non-ADC; SCLC (Histology)
data1$ASIA_preop = factor(data1$ASIA_preop) #Pre. ASIA score: A, C, D, E (ASIA_preop)
data1$ASIA_PO = factor(data1$ASIA_PO) #Post. ASIA score: C, D, E (ASIA_PO)
data1$Tokuhashi_cat = factor(data1$Tokuhashi_cat) #Revised Tokuhashi score: 0-8; 9-11 (Tokuhashi_cat)
data1$ambulation_preop = factor(data1$ambulation_preop) #Pre. ambulatory status: No, With help; Independent (ambulation_preop)
data1$ambulationPO = factor(data1$ambulationPO) #Post. ambulatory status: No, With help; Independent (ambulationPO)
data1$RoRx_PO = factor(data1$RoRx_PO) #Post. Radiotherapy: No- Yes (RoRx_PO)
data1$chimiotxPO = factor(data1$chimiotxPO) #Post. Chemotherapy: No- Yes (chimiotxPO)
####View new dataset again
head(data1)
## ID sexe age tabagisme Histology ASIA_preop ASIA_PO Tokuhashi_cat
## 1 1 1 55 1 1 D D 1
## 2 2 1 71 0 1 E E 2
## 3 3 2 68 1 1 D D 1
## 4 4 1 71 1 1 D D 1
## 5 5 1 63 1 1 E E 1
## 6 6 2 65 1 1 D E 1
## ambulation_preop ambulationPO RoRx_PO chimiotxPO SurviePosOpMOIS deces age_gr
## 1 2 2 1 1 4.77 0 1
## 2 1 1 0 0 4.07 0 0
## 3 2 2 1 1 8.13 0 0
## 4 2 2 1 0 1.93 1 0
## 5 1 1 1 0 1.57 1 0
## 6 2 1 0 1 12.67 1 0
View(data1)
###Outcome: SurviePosOpMOIS: time-to-event deces: 1: censorted (alive); 0 event (die)
baseline = Surv(data1$SurviePosOpMOIS, data1$deces==0)
km = survfit(baseline ~ 1)
summary(km)
## Call: survfit(formula = baseline ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4.07 44 1 0.977 0.0225 0.934 1
## 4.77 36 1 0.950 0.0345 0.885 1
## 8.13 22 1 0.907 0.0536 0.808 1
## 28.67 5 1 0.726 0.1678 0.461 1
Kaplan – Meier plot:
plot(km, xlab="Time to death", ylab="Prob of survival")
#table 1: general chracteristics of participants
#install package "moonBook", ztable
library(moonBook)
require(ztable)
## Loading required package: ztable
## Welcome to package ztable ver 0.2.0
require(magrittr)
## Loading required package: magrittr
options(ztable.type="html")
#table 1: general chracteristics of participants
mytable(data1)
##
## Descriptive Statistics
## -----------------------------------
## N Total
## -----------------------------------
## ID 87 44.0 ± 25.3
## sexe 87
## - 1 45 (51.7%)
## - 2 42 (48.3%)
## age 87 61.3 ± 8.8
## tabagisme 87
## - 0 35 (40.2%)
## - 1 52 (59.8%)
## Histology 87
## - 0 22 (25.3%)
## - 1 58 (66.7%)
## - 2 7 (8.0%)
## ASIA_preop 87
## - A 1 (1.1%)
## - C 5 (5.7%)
## - D 39 (44.8%)
## - E 42 (48.3%)
## ASIA_PO 87
## - C 2 (2.3%)
## - D 27 (31.0%)
## - E 58 (66.7%)
## Tokuhashi_cat 87
## - 1 72 (82.8%)
## - 2 15 (17.2%)
## ambulation_preop 87
## - 0 22 (25.3%)
## - 1 41 (47.1%)
## - 2 24 (27.6%)
## ambulationPO 87
## - 0 4 (4.6%)
## - 1 48 (55.2%)
## - 2 35 (40.2%)
## RoRx_PO 87
## - 0 28 (32.2%)
## - 1 59 (67.8%)
## chimiotxPO 87
## - 0 54 (62.1%)
## - 1 33 (37.9%)
## SurviePosOpMOIS 87 7.5 ± 12.0
## deces 87
## - 0 4 (4.6%)
## - 1 83 (95.4%)
## age_gr 87
## - 0 51 (58.6%)
## - 1 36 (41.4%)
## -----------------------------------
#Univariate cox-model
library(survival)
#model 1: sex
cox1 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$sexe)
summary(cox1)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$sexe)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$sexe2 -0.5513 0.5762 1.0309 -0.535 0.593
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$sexe2 0.5762 1.736 0.0764 4.346
##
## Concordance= 0.665 (se = 0.097 )
## Likelihood ratio test= 0.28 on 1 df, p=0.6
## Wald test = 0.29 on 1 df, p=0.6
## Score (logrank) test = 0.29 on 1 df, p=0.6
#model 2: age group
cox2 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$age_gr)
summary(cox2)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$age_gr)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$age_gr1 -0.5595 0.5715 1.2260 -0.456 0.648
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$age_gr1 0.5715 1.75 0.05169 6.319
##
## Concordance= 0.558 (se = 0.148 )
## Likelihood ratio test= 0.22 on 1 df, p=0.6
## Wald test = 0.21 on 1 df, p=0.6
## Score (logrank) test = 0.21 on 1 df, p=0.6
#model 3: tobacco use
cox3 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$tabagisme)
summary(cox3)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$tabagisme)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$tabagisme1 1.635 5.129 1.194 1.369 0.171
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$tabagisme1 5.129 0.195 0.4937 53.29
##
## Concordance= 0.568 (se = 0.162 )
## Likelihood ratio test= 2.22 on 1 df, p=0.1
## Wald test = 1.87 on 1 df, p=0.2
## Score (logrank) test = 2.23 on 1 df, p=0.1
#Model 4: Histology
cox4 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$Histology)
## Warning in fitter(X, Y, strats, offset, init, control, weights = weights, :
## Loglik converged before variable 1 ; coefficient may be infinite.
summary(cox4)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$Histology)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$Histology1 1.965e+01 3.411e+08 1.548e+04 0.001 0.999
## data1$Histology2 1.120e-01 1.119e+00 3.365e+04 0.000 1.000
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$Histology1 3.411e+08 2.932e-09 0 Inf
## data1$Histology2 1.118e+00 8.941e-01 0 Inf
##
## Concordance= 0.646 (se = 0.038 )
## Likelihood ratio test= 2.89 on 2 df, p=0.2
## Wald test = 0 on 2 df, p=1
## Score (logrank) test = 1.75 on 2 df, p=0.4
#model 5: Pre. ASIA score
cox5 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$ASIA_preop)
## Warning in fitter(X, Y, strats, offset, init, control, weights = weights, :
## Loglik converged before variable 1 ; coefficient may be infinite.
summary(cox5)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$ASIA_preop)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$ASIA_preopC -1.562e+01 1.652e-07 1.408e+04 -0.001 0.999
## data1$ASIA_preopD 9.303e-01 2.535e+00 1.227e+00 0.758 0.448
## data1$ASIA_preopE NA NA 0.000e+00 NA NA
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$ASIA_preopC 1.652e-07 6.052e+06 0.0000 Inf
## data1$ASIA_preopD 2.535e+00 3.944e-01 0.2288 28.1
## data1$ASIA_preopE NA NA NA NA
##
## Concordance= 0.568 (se = 0.157 )
## Likelihood ratio test= 0.72 on 2 df, p=0.7
## Wald test = 0.57 on 2 df, p=0.8
## Score (logrank) test = 0.68 on 2 df, p=0.7
#model 6: Post. ASIA score
cox6 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$ASIA_PO)
## Warning in fitter(X, Y, strats, offset, init, control, weights = weights, :
## Loglik converged before variable 1,2 ; coefficient may be infinite.
summary(cox6)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$ASIA_PO)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$ASIA_POD 1.733e+01 3.368e+07 1.429e+04 0.001 0.999
## data1$ASIA_POE 1.519e+01 3.950e+06 1.429e+04 0.001 0.999
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$ASIA_POD 33676279 2.969e-08 0 Inf
## data1$ASIA_POE 3950104 2.532e-07 0 Inf
##
## Concordance= 0.684 (se = 0.155 )
## Likelihood ratio test= 3.28 on 2 df, p=0.2
## Wald test = 2.99 on 2 df, p=0.2
## Score (logrank) test = 4.37 on 2 df, p=0.1
#model 7: Revised Tokuhashi score
cox7 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$Tokuhashi_cat)
summary(cox7)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$Tokuhashi_cat)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$Tokuhashi_cat2 0.4573 1.5798 1.1112 0.412 0.681
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$Tokuhashi_cat2 1.58 0.633 0.179 13.95
##
## Concordance= 0.587 (se = 0.161 )
## Likelihood ratio test= 0.17 on 1 df, p=0.7
## Wald test = 0.17 on 1 df, p=0.7
## Score (logrank) test = 0.17 on 1 df, p=0.7
#model 8: Pre. ambulatory status
cox8 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$ambulation_preop)
## Warning in fitter(X, Y, strats, offset, init, control, weights = weights, :
## Loglik converged before variable 1,2 ; coefficient may be infinite.
summary(cox8)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$ambulation_preop)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$ambulation_preop1 1.777e+01 5.239e+07 1.731e+04 0.001 0.999
## data1$ambulation_preop2 1.969e+01 3.546e+08 1.731e+04 0.001 0.999
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$ambulation_preop1 52385120 1.909e-08 0 Inf
## data1$ambulation_preop2 354636157 2.820e-09 0 Inf
##
## Concordance= 0.694 (se = 0.127 )
## Likelihood ratio test= 4.2 on 2 df, p=0.1
## Wald test = 2.67 on 2 df, p=0.3
## Score (logrank) test = 4.52 on 2 df, p=0.1
#model 9: Post. ambulatory status
cox9 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$ambulationPO)
summary(cox9)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$ambulationPO)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$ambulationPO1 -1.7728 0.1699 1.2388 -1.431 0.152
## data1$ambulationPO2 NA NA 0.0000 NA NA
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$ambulationPO1 0.1699 5.887 0.01498 1.925
## data1$ambulationPO2 NA NA NA NA
##
## Concordance= 0.646 (se = 0.153 )
## Likelihood ratio test= 2.21 on 1 df, p=0.1
## Wald test = 2.05 on 1 df, p=0.2
## Score (logrank) test = 2.6 on 1 df, p=0.1
#model 10: Post. Radiotherapy
cox10 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$RoRx_PO)
summary(cox10)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$RoRx_PO)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$RoRx_PO1 -0.7463 0.4741 1.2372 -0.603 0.546
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$RoRx_PO1 0.4741 2.109 0.04196 5.357
##
## Concordance= 0.612 (se = 0.14 )
## Likelihood ratio test= 0.33 on 1 df, p=0.6
## Wald test = 0.36 on 1 df, p=0.5
## Score (logrank) test = 0.38 on 1 df, p=0.5
#model 11: Post. Chemotherapy
cox11 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$chimiotxPO)
summary(cox11)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$chimiotxPO)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$chimiotxPO1 -0.05409 0.94735 1.23130 -0.044 0.965
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$chimiotxPO1 0.9474 1.056 0.08481 10.58
##
## Concordance= 0.549 (se = 0.138 )
## Likelihood ratio test= 0 on 1 df, p=1
## Wald test = 0 on 1 df, p=1
## Score (logrank) test = 0 on 1 df, p=1
#multivariate model
#model 12: multivariate model
cox12 = coxph(Surv(data1$SurviePosOpMOIS, data1$deces==0) ~ data1$sexe+data1$age_gr+data1$RoRx_PO +data1$tabagisme+ data1$Histology+data1$ASIA_PO+data1$ASIA_preop +data1$Tokuhashi_cat+ data1$ambulation_preop+ data1$ambulationPO+data1$RoRx_PO+data1$chimiotxPO)
## Warning in fitter(X, Y, strats, offset, init, control, weights = weights, : Ran
## out of iterations and did not converge
summary(cox12)
## Call:
## coxph(formula = Surv(data1$SurviePosOpMOIS, data1$deces == 0) ~
## data1$sexe + data1$age_gr + data1$RoRx_PO + data1$tabagisme +
## data1$Histology + data1$ASIA_PO + data1$ASIA_preop +
## data1$Tokuhashi_cat + data1$ambulation_preop + data1$ambulationPO +
## data1$RoRx_PO + data1$chimiotxPO)
##
## n= 87, number of events= 4
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data1$sexe2 -4.301e+01 2.085e-19 2.737e+03 -0.016 0.987
## data1$age_gr1 -2.617e+01 4.324e-12 2.684e+03 -0.010 0.992
## data1$RoRx_PO1 -4.354e+01 1.238e-19 4.676e+03 -0.009 0.993
## data1$tabagisme1 3.550e+01 2.603e+15 2.956e+03 0.012 0.990
## data1$Histology1 -3.404e+01 1.651e-15 7.646e+03 -0.004 0.996
## data1$Histology2 -1.116e+02 3.350e-49 2.618e+04 -0.004 0.997
## data1$ASIA_POD -1.255e+02 3.030e-55 3.176e+03 -0.040 0.968
## data1$ASIA_POE -1.827e+02 4.454e-80 3.177e+03 -0.058 0.954
## data1$ASIA_preopC -1.810e+01 1.373e-08 1.788e+05 0.000 1.000
## data1$ASIA_preopD -2.928e+01 1.926e-13 3.301e+03 -0.009 0.993
## data1$ASIA_preopE -1.334e+01 1.601e-06 3.302e+03 -0.004 0.997
## data1$Tokuhashi_cat2 5.517e+01 9.158e+23 3.467e+03 0.016 0.987
## data1$ambulation_preop1 2.257e+01 6.349e+09 2.956e+03 0.008 0.994
## data1$ambulation_preop2 4.683e+01 2.173e+20 2.921e+03 0.016 0.987
## data1$ambulationPO1 -2.708e+00 6.666e-02 3.193e+03 -0.001 0.999
## data1$ambulationPO2 2.373e+00 1.073e+01 3.193e+03 0.001 0.999
## data1$chimiotxPO1 1.513e+00 4.541e+00 4.296e+03 0.000 1.000
##
## exp(coef) exp(-coef) lower .95 upper .95
## data1$sexe2 2.085e-19 4.797e+18 0 Inf
## data1$age_gr1 4.324e-12 2.313e+11 0 Inf
## data1$RoRx_PO1 1.238e-19 8.079e+18 0 Inf
## data1$tabagisme1 2.603e+15 3.841e-16 0 Inf
## data1$Histology1 1.651e-15 6.056e+14 0 Inf
## data1$Histology2 3.350e-49 2.985e+48 0 Inf
## data1$ASIA_POD 3.030e-55 3.301e+54 0 Inf
## data1$ASIA_POE 4.454e-80 2.245e+79 0 Inf
## data1$ASIA_preopC 1.373e-08 7.284e+07 0 Inf
## data1$ASIA_preopD 1.926e-13 5.192e+12 0 Inf
## data1$ASIA_preopE 1.601e-06 6.246e+05 0 Inf
## data1$Tokuhashi_cat2 9.158e+23 1.092e-24 0 Inf
## data1$ambulation_preop1 6.349e+09 1.575e-10 0 Inf
## data1$ambulation_preop2 2.173e+20 4.601e-21 0 Inf
## data1$ambulationPO1 6.666e-02 1.500e+01 0 Inf
## data1$ambulationPO2 1.073e+01 9.324e-02 0 Inf
## data1$chimiotxPO1 4.541e+00 2.202e-01 0 Inf
##
## Concordance= 1 (se = 0 )
## Likelihood ratio test= 24.14 on 17 df, p=0.1
## Wald test = 0.01 on 17 df, p=1
## Score (logrank) test = 15.5 on 17 df, p=0.6