WeightedSurv/hccweightedcode.R at main · Ghoshlab/WeightedSurv · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
# Weighted code


hcctcga <- read.csv2("livercancerTCGAclinical.csv",header=T,sep=",")
hccbclc <- read.csv("HCCBCLC.csv",header=T)
names(hccbclc)[8] <- "age"
names(hccbclc)[9] <- "gender"
hccbclc$gender <- ifelse(hccbclc$gender == 1,"FEMALE","MALE")
hcctcga$gender <- as.factor(hcctcga$gender)
hccbclc$gender <- as.factor(hccbclc$gender)
hccbclc$futime <- hccbclc$OS
hccbclc$fustat <- hccbclc$D5
library(survival)
library(randomForestSRC)
library(survminer)

fit <- survfit(Surv(futime,fustat)~gender,data=hcctcga)
ggsurvplot(fit,risk.table=T,tables.height=0.3)

fit <- survfit(Surv(futime,fustat)~gender,data=hccbclc)
ggsurvplot(fit,risk.table=T,tables.height=0.3)


cph1 <- coxph(Surv(futime,fustat)~age+gender,data=hcctcga)
wts <- predict(cph1,newdata=hccbclc,type="risk")

cph2 <- coxph(Surv(OS,D5)~BCLC+ER+LR,data=hccbclc,weights = wts)
# Now with no weights
cph3 <- coxph(Surv(OS,D5)~BCLC+ER+LR,data=hccbclc)
cph4 <- coxph(Surv(OS,D5)~BCLC+ER+LR+tmp,data=hccbclc)
n <- length(wts)
# Try new approach
B <- 1000
coefs <- matrix(0,nrow=B,ncol=3)
for (i in 1:B) {
  set.seed(520+i)
  simy <- rexp(n)
  pwts <- simy*wts
  tmpcph <-  coxph(Surv(OS,D5)~BCLC+ER+LR,data=hccbclc,weights = pwts)
  coefs[i,] <- coef(tmpcph)
  cat(i,"\n")
}


# Now try random forests

obj <- rfsrc(Surv(futime,fustat)~age+gender,data=hcctcga,ntree=1000)
wts <- predict(obj,hccbclc)$chf[,120]

cph2 <- coxph(Surv(OS,D5)~BCLC+ER+LR,data=hccbclc,weights = wts)
# Try new approach
B <- 1000
coefs <- matrix(0,nrow=B,ncol=3)
n <- length(wts)
for (i in 1:B) {
  set.seed(520+i)
  simy <- rexp(n)
  pwts <- simy*wts
  tmpcph <-  coxph(Surv(OS,D5)~BCLC+ER+LR,data=hccbclc,weights = pwts)
  coefs[i,] <- coef(tmpcph)
  cat(i,"\n")
}

# next, additive hazards with proportional hazards weights
ah1 <- aalen(Surv(OS,D5)~const(BCLC)+const(ER)+const(LR),data=hccbclc)

cph1 <- coxph(Surv(futime,fustat)~age+gender,data=hcctcga)
wts <- predict(cph1,newdata=hccbclc,type="risk")

ah2 <- aalen(Surv(OS,D5)~const(BCLC)+const(ER)+const(LR),data=hccbclc,weights=wts)
# Try new approach
B <- 1000
coefs <- matrix(0,nrow=B,ncol=3)
n <- length(wts)
for (i in 1:B) {
  set.seed(520+i)
  simy <- rexp(n)
  pwts <- simy*wts
  tmpah<-  aalen(Surv(OS,D5)~const(BCLC)+const(ER)+const(LR),
                 data=hccbclc,weights=pwts)
  coefs[i,] <- coef(tmpah)[,1]
  cat(i,"\n")
}

# Next, additive hazards with random forest weights
obj <- rfsrc(Surv(futime,fustat)~age+gender,data=hcctcga,ntree=1000)
wts <- predict(obj,hccbclc)$chf[,120]
ah3 <- aalen(Surv(OS,D5)~const(BCLC)+const(ER)+const(LR),data=hccbclc,weights=wts)
# Try new approach
B <- 1000
coefs <- matrix(0,nrow=B,ncol=3)
n <- length(wts)
for (i in 1:B) {
  set.seed(520+i)
  simy <- rexp(n)
  pwts <- simy*wts
  tmpah<-  aalen(Surv(OS,D5)~const(BCLC)+const(ER)+const(LR),
                 data=hccbclc,weights=pwts)
  coefs[i,] <- coef(tmpah)[,1]
  cat(i,"\n")
}