Lili Wang 2019-08-06
This R package is to analyze bivariate survival outcomes. Here specifically, we implement this method to analyze alternating recurrent events using a bivariate correlated frailty model. Both the regression parameters and the variance-covariance matrix will be estimated. This R package allows data to have many clusters. The estimation procedure is similar to the traditional PPL method as established in coxme, while we do not require knowing the correlation direction between the correlated two events.
This R package will be improved and upgraded in the near future.
To install the R package from Github, you will need to install another R package "devtools". Please uncomment the codes to install them.
# install.packages("devtools") # library(devtools) # install_github("lilywang1988/BivPPL") library(BivPPL) #> Loading required package: survival #> Loading required package: Matrix #> Loading required package: mvtnormGenerate alternating recurrent event data using gen.data, beta1 and beta2 are regression parameters for the two events, theta is defining the 3 different entries of the variance-covariance matrix of the two correlated frailties. N is the number of clusters or subjects in the framework of alternating recurrent events.
set.seed(100) N<-100 # number of clusters beta1 <- c(0.5,-0.3,0.5) # regression parameters for event type 1 beta2 <- c(0.8,-0.2,0.3) # regression parameters for event type 2 beta <- c(beta1,beta2) theta <- c(0.25,0.25,-0.125) # variance-covariance matrix for the bivariate frailty (denoted as D), it is a vector (D[1,1],D[2,2], D[1,2]) lambda01 <- 1 lambda02 <- 1 cen <- 10 # maximum censoring time centype <- TRUE # fixed censoring time at cen data <- gen.data(N,beta1,beta2,theta,lambda01,lambda02,c=cen,Ctype=centype) ptm<-proc.time() res <- BivPPL(data) proc.time() - ptm #> user system elapsed #> 20.534 0.917 21.703 res$beta_hat #> [1] 0.52126635 -0.26535133 0.45676406 0.74664981 -0.06705498 0.24962866 res$beta_ASE #> [1] 0.05652019 0.05266173 0.05729070 0.06239355 0.05456852 0.05988814 res$D_hat #> r1_hat r2_hat #> r1_hat 0.3791955 -0.2053474 #> r2_hat -0.2053474 0.2104612 # fit the model assuming the independence between the two frailties ptm<-proc.time() res2 <- BivPPL(data,independence=T) proc.time() - ptm #> user system elapsed #> 9.938 0.525 10.718 res2$beta_hat #> [1] 0.51974445 -0.26578661 0.44356750 0.73519612 -0.06245907 0.23199444 res2$beta_ASE #> [1] 0.05700419 0.05297851 0.05777039 0.06270781 0.05549611 0.06043905 res2$D_hat #> [,1] [,2] #> [1,] 0.3742242 0.0000000 #> [2,] 0.0000000 0.1866389 # A likelihood ratio test LRT<-2*(res$LogMargProb-res2$LogMargProb) print(round(pchisq(abs(LRT),df=1,lower.tail = F),3)) #> [,1] #> [1,] 0.002 # sparsen the Hessian matrix ptm<-proc.time() res3 <- BivPPL(data,huge=TRUE) proc.time() - ptm #> user system elapsed #> 13.449 0.110 13.773 res3$beta_hat #> z11 z12 z13 z21 z22 z23 #> 0.52042943 -0.26522250 0.45584999 0.74594886 -0.06716247 0.24938898 res3$beta_ASE #> [1] 0.05649124 0.05267874 0.05726610 0.06227382 0.05454575 0.05985640 res3$D_hat #> r1_hat r2_hat #> r1_hat 0.3728382 -0.2018795 #> r2_hat -0.2018795 0.2062161 # Compare with coxme # install.packages("coxme") library(coxme) #> Loading required package: bdsmatrix #> #> Attaching package: 'bdsmatrix' #> The following object is masked from 'package:base': #> #> backsolve data_coxme <- tocoxme(data) # assume that we do not know the two events are negatively associated and we transform the data to be positively associated ptm<-proc.time() res_coxme <- coxme(Surv(data_coxme$time,data_coxme$delta)~data_coxme$Z1+data_coxme$Z2+(1|data_coxme$b0)+(1|data_coxme$b1)+(1|data_coxme$b2)+strata(data_coxme$joint)) # disable the Hessian matrix sparsening and likelihood refining proc.time() - ptm #> user system elapsed #> 48.536 0.224 48.857 res_coxme$coefficients #> data_coxme$Z11 data_coxme$Z12 data_coxme$Z13 data_coxme$Z21 data_coxme$Z22 #> 0.52065085 -0.26593336 0.44450391 0.73457314 -0.06255747 #> data_coxme$Z23 #> 0.23176550 sqrt(diag(vcov(res_coxme))) #> [1] 0.05707456 0.05301526 0.05782401 0.06267015 0.05545393 0.06039687 assemble(as.vector(unlist(res_coxme$vcoef))) #> [,1] [,2] #> [1,] 0.3809510825 0.0003215251 #> [2,] 0.0003215251 0.1838460080 data_coxme_n <- tocoxme_n(data) # assume that we know the two events are negatively associated ptm<-proc.time() res_coxme_n<- coxme(Surv(data_coxme_n$time,data_coxme_n$delta)~data_coxme_n$Z1+data_coxme_n$Z2+(data_coxme_n$Z0|data_coxme_n$b0)+(1|data_coxme_n$b1)+(1|data_coxme_n$b2)+strata(data_coxme_n$joint)) proc.time() - ptm #> user system elapsed #> 57.075 0.256 57.395 res_coxme_n$coefficients #> data_coxme_n$Z11 data_coxme_n$Z12 data_coxme_n$Z13 data_coxme_n$Z21 #> 0.52257408 -0.26556525 0.45768187 0.74671591 #> data_coxme_n$Z22 data_coxme_n$Z23 #> -0.06678617 0.24903043 sqrt(diag(vcov(res_coxme_n))) #> [1] 0.05664816 0.05273497 0.05739250 0.06244051 0.05465062 0.05995059 assemble_n(as.vector(unlist(res_coxme_n$vcoef))) #> [,1] [,2] #> [1,] 0.3861649 -0.2004567 #> [2,] -0.2004567 0.2103093