Load the library
library(survC1)
read the sample data (PBC data in survival package)
D=CompCase(pbc[1:200,c(2:4,10:14)])
D[,2]=as.numeric(D[,2]==2)
Lisiting a part of the data (outcome and predictors)
head(D)
## time status trt edema bili chol albumin copper
## 1 400 1 1 1.0 14.5 261 2.60 156
## 2 4500 0 1 0.0 1.1 302 4.14 54
## 3 1012 1 1 0.5 1.4 176 3.48 210
## 4 1925 1 1 0.5 1.8 244 2.54 64
## 5 1504 0 2 0.0 3.4 279 3.53 143
## 6 2503 1 2 0.0 0.8 248 3.98 50
see KM plots of event time and censoring time
par(mfrow=c(1,2))
plot(survfit(Surv(time/365.25, status)~1, data=D), main="Survival time")
plot(survfit(Surv(time/365.25, status==0)~1, data=D), main="Censoring time")
Inference for C –> Inf.Cval() function
tau=365.25*8
C=Inf.Cval(D, tau, itr=200)
C$ft
## Call:
## coxph(formula = Surv(mydata[, 1], mydata[, 2]) ~ covs)
##
##
## coef exp(coef) se(coef) z p
## covstrt -0.014813 0.985 0.213216 -0.0695 9.4e-01
## covsedema 0.835274 2.305 0.400275 2.0867 3.7e-02
## covsbili 0.081894 1.085 0.021980 3.7258 1.9e-04
## covschol 0.000101 1.000 0.000449 0.2253 8.2e-01
## covsalbumin -1.031226 0.357 0.311795 -3.3074 9.4e-04
## covscopper 0.004857 1.005 0.001030 4.7170 2.4e-06
##
## Likelihood ratio test=99.7 on 6 df, p=0 n= 177, number of events= 94
round(c(C$Dhat, C$se, C$low95, C$upp95), digits=3)
## [1] 0.776 0.027 0.722 0.830
Inference of Delta C between 2 models –> Inf.Cval.Delta() function
model0<-D[,c(1:2,4:5)] ;
model1<-D
covs1<-as.matrix(model1[,c(-1,-2)])
covs0<-as.matrix(model0[,c(-1,-2)])
head(covs1)
## trt edema bili chol albumin copper
## 1 1 1.0 14.5 261 2.60 156
## 2 1 0.0 1.1 302 4.14 54
## 3 1 0.5 1.4 176 3.48 210
## 4 1 0.5 1.8 244 2.54 64
## 5 2 0.0 3.4 279 3.53 143
## 6 2 0.0 0.8 248 3.98 50
head(covs0)
## edema bili
## 1 1.0 14.5
## 2 0.0 1.1
## 3 0.5 1.4
## 4 0.5 1.8
## 5 0.0 3.4
## 6 0.0 0.8
Delta=Inf.Cval.Delta(model0[,1:2], covs0, covs1, tau, itr=200)
round(Delta, digits=3)
## Est SE Lower95 Upper95
## Model1 0.776 0.026 0.725 0.827
## Model0 0.775 0.027 0.723 0.827
## Delta 0.000 0.021 -0.040 0.041
Point estimation via cross-validation –> cvC() function
head(model0)
## time status edema bili
## 1 400 1 1.0 14.5
## 2 4500 0 0.0 1.1
## 3 1012 1 0.5 1.4
## 4 1925 1 0.5 1.8
## 5 1504 0 0.0 3.4
## 6 2503 1 0.0 0.8
cvC(model0,tau,cvK=2,Rep=100)
## [1] 0.7748626