cmprsk-note
Li-Hsin Chien
2023-02-03
cmprsk
References: Fine JP and Gray RJ (1999) A proportional hazards model for the subdistribution of a competing risk. JASA 94:496-509.
Example: from “help(crr)”
#install.packages("cmprsk")
library(cmprsk)## Warning: 套件 'cmprsk' 是用 R 版本 4.2.2 來建造的## 載入需要的套件：survivalsimulated data to test
sample size = 200, 3 covariates \(Z=(x_1,x_2,x_3)\);
status = 0 (censored), 1(target event, ex.cancer death), 2(competing event, ex.death other than cancer)
set.seed(10)
ftime <- rexp(200) # observed survival time
fstatus <- sample(0:2,200,replace=TRUE) #0:censored, 1:target event
cov <- matrix(runif(600),nrow=200) # p=3, covariates
dimnames(cov)[[2]] <- c('x1','x2','x3')
head(cov)## x1 x2 x3
## [1,] 0.9060413 0.49313019 0.56228712
## [2,] 0.1501874 0.11733963 0.16163796
## [3,] 0.9916424 0.32289117 0.80716830
## [4,] 0.4817376 0.04328842 0.08930905
## [5,] 0.4729846 0.92961596 0.48674353
## [6,] 0.1389506 0.02522151 0.86876462\(\lambda_1(t;Z)=\lambda_{01}(t)exp(x_1\beta_1+x_2\beta_2+x_3\beta_3)\)
print(z <- crr(ftime,fstatus,cov)) ## convergence: TRUE
## coefficients:
## x1 x2 x3
## 0.26680 -0.05568 0.28050
## standard errors:
## [1] 0.4211 0.3812 0.3810
## two-sided p-values:
## x1 x2 x3
## 0.53 0.88 0.46summary(z) ## Competing Risks Regression
##
## Call:
## crr(ftime = ftime, fstatus = fstatus, cov1 = cov)
##
## coef exp(coef) se(coef) z p-value
## x1 0.2668 1.306 0.421 0.633 0.53
## x2 -0.0557 0.946 0.381 -0.146 0.88
## x3 0.2805 1.324 0.381 0.736 0.46
##
## exp(coef) exp(-coef) 2.5% 97.5%
## x1 1.306 0.766 0.572 2.98
## x2 0.946 1.057 0.448 2.00
## x3 1.324 0.755 0.627 2.79
##
## Num. cases = 200
## Pseudo Log-likelihood = -320
## Pseudo likelihood ratio test = 1.02 on 3 df,\((\hat{\beta}_1,\hat{\beta}_2,\hat{\beta}_3)=(0.2668,-0.0557,0.2805)\)
Suppose 2 individual with covariates \(Z_1=(x_{11},x_{12},x_{13})=(.1,.5,.8),Z_2=(.1,.5,.2)\). The cumulative incidence function of the target event (cancer death): for \(j=1,2\) \[\hat{F}_1(t;Z_j)=1-exp(-\hat{\Lambda}_1(t;Z_j)), \text{where } \hat{\Lambda}_1(t;Z_j)=\hat{\Lambda}_{10}(t)exp(x_{j1}\hat{\beta}_1+x_{j2}\hat{\beta}_2+x_{j3}\hat{\beta}_3)\]
z.p <- predict(z,rbind(c(.1,.5,.8),c(.1,.5,.2)))
plot(z.p,lty=1,color=2:3) 
head(z.p)## [,1] [,2] [,3]
## [1,] 0.01495641 0.004854511 0.004104118
## [2,] 0.03280299 0.009736033 0.008234200
## [3,] 0.09663135 0.014721413 0.012455409
## [4,] 0.13085151 0.019734443 0.016703367
## [5,] 0.14671290 0.024772259 0.020975720
## [6,] 0.15447993 0.029812790 0.025253797First column: vecotr of unique failure times, \(t\)
Second column: the cumulative incidence rate for the 1st individual,\(\hat{F}_1(t;Z_1)\)
third column: the cumulative incidence rate for the 2nd individual,\(\hat{F}_1(t;Z_2)\)
Z1<-c(.1,.5,.8)
H1.base<-cumsum(z$bfitj)
H1<-H1.base*exp(sum(Z1*z$coef))
cumu.inci.1<-1-exp(-H1)H1.base = \(\hat{\Lambda}_{10}(t)\)
H1 = \(\hat{\Lambda}_1(t;Z_1)=\hat{\Lambda}_{10}(t)exp(x_{11}\hat{\beta}_1+x_{12}\hat{\beta}_2+x_{13}\hat{\beta}_3)\)
cumu.inci.1 = \(\hat{F}_1(t;Z_1)=1-exp(-\hat{\Lambda}_1(t;Z_1))\)
head(cumu.inci.1)## [1] 0.004854511 0.009736033 0.014721413 0.019734443 0.024772259 0.029812790head(z.p)## [,1] [,2] [,3]
## [1,] 0.01495641 0.004854511 0.004104118
## [2,] 0.03280299 0.009736033 0.008234200
## [3,] 0.09663135 0.014721413 0.012455409
## [4,] 0.13085151 0.019734443 0.016703367
## [5,] 0.14671290 0.024772259 0.020975720
## [6,] 0.15447993 0.029812790 0.025253797The cumulative incidence rate when \(t=1\), \(\hat{F}_1(t=1;Z_j)\)
tmp<-sum(z.p[,1]<=1)
z.p[(tmp-1):(tmp+1),]## [,1] [,2] [,3]
## [1,] 0.9278297 0.2046290 0.1759169
## [2,] 0.9536605 0.2110094 0.1815072
## [3,] 1.0146017 0.2173996 0.1871130plot(z.p,lty=1,color=2:3)
abline(v=1,col="gray")
abline(h=0.2110094,col=2,lty=2)
abline(h=0.1815072,col=3,lty=2)
\(\hat{F}_1(t=1;Z_1)=0.2110094\)
\(\hat{F}_1(t=1;Z_2)=0.1815072\)
Confidence interval
\(\hat{\beta}\) ~ \(N(\beta, \hat{\Sigma}_{\hat{\beta}})\)
z$coef #beta_hat## x1 x2 x3
## 0.26675631 -0.05568438 0.28049017z$var #sigma_hat## [,1] [,2] [,3]
## [1,] 0.17736467 0.01485832 -0.04109189
## [2,] 0.01485832 0.14534173 -0.01394836
## [3,] -0.04109189 -0.01394836 0.14514930sqrt(diag(z$var)) #standard error of beta_hat## [1] 0.4211469 0.3812371 0.3809846summary(z)## Competing Risks Regression
##
## Call:
## crr(ftime = ftime, fstatus = fstatus, cov1 = cov)
##
## coef exp(coef) se(coef) z p-value
## x1 0.2668 1.306 0.421 0.633 0.53
## x2 -0.0557 0.946 0.381 -0.146 0.88
## x3 0.2805 1.324 0.381 0.736 0.46
##
## exp(coef) exp(-coef) 2.5% 97.5%
## x1 1.306 0.766 0.572 2.98
## x2 0.946 1.057 0.448 2.00
## x3 1.324 0.755 0.627 2.79
##
## Num. cases = 200
## Pseudo Log-likelihood = -320
## Pseudo likelihood ratio test = 1.02 on 3 df,Confidence interval for \(\hat{s}_1=x_{11}\hat{\beta}_1+x_{12}\hat{\beta}_2+x_{13}\hat{\beta}_3\) \[\hat{s}_1=Z_1^T\hat{\beta}\] \[sd(\hat{s}_1)=\sqrt{Var({\hat{s}_1})}=\sqrt{Z_1^T\hat{\Sigma}_{\hat{\beta}}Z_1)}\]
s_hat<-sum(Z1*z$coef)
s_hat_sd<-sqrt(t(as.matrix(Z1)) %*% z$var %*%(as.matrix(Z1)))
s.ci.u<-s_hat+qnorm(0.975)*s_hat_sd
s.ci.d<-s_hat-qnorm(0.975)*s_hat_sd
c(s.ci.d,s.ci.u) # 95% CI for s_hat, suppose fixed baseline hazard function## [1] -0.4407282 0.8871793Confidence interval for \(\hat{F}_1(t;Z_1)\) (given the baseline hazard function, \(\hat{\Lambda}_{10}(t)\), fixed)
Recall that \[\hat{F}_1(t;Z_1)=1-exp(-\hat{\Lambda}_1(t;Z_1)), \text{where } \hat{\Lambda}_1(t;Z_1)=\hat{\Lambda}_{10}(t)exp(x_{11}\hat{\beta}_1+x_{12}\hat{\beta}_2+x_{13}\hat{\beta}_3)=\hat{\Lambda}_{10}(t)exp(\hat{s}_1)\]
H1_u<-H1.base*exp(c(s.ci.u))
H1_d<-H1.base*exp(c(s.ci.d))
cumu.inci.1.u<-1-exp(-H1_u)
cumu.inci.1.d<-1-exp(-H1_d)
head(cbind(z.p[,1:2],cumu.inci.1.d,cumu.inci.1.u))## cumu.inci.1.d cumu.inci.1.u
## [1,] 0.01495641 0.004854511 0.002502110 0.009408107
## [2,] 0.03280299 0.009736033 0.005024121 0.018825045
## [3,] 0.09663135 0.014721413 0.007606023 0.028397284
## [4,] 0.13085151 0.019734443 0.010208645 0.037976697
## [5,] 0.14671290 0.024772259 0.012830649 0.047557071
## [6,] 0.15447993 0.029812790 0.015460650 0.057096040sfun<-stepfun(z.p[,1],c(0,z.p[,2]),right=F)
sfun.u<-stepfun(z.p[,1],c(0,cumu.inci.1.u),right=F)
sfun.d<-stepfun(z.p[,1],c(0,cumu.inci.1.d),right=F)
plot(sfun,do.point=F)
plot(sfun.u,do.point=F,add=T,col="gray")
plot(sfun.d,do.point=F,add=T,col="gray")