In lilywang1988/BivPPL: Bivariate frailty model estimation

knitr::opts_chunk$set(collapse = TRUE,
                     comment = "#>",
                     fig.width=12, fig.height=8,
                     fig.path = "man/figures/README-")

Purpose

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.

Install the package

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)

Vignettes

Generate 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
res$beta_hat
res$beta_ASE
res$D_hat


# fit the model assuming the independence between the two frailties
ptm<-proc.time()
res2 <- BivPPL(data,independence=T)
proc.time() - ptm
res2$beta_hat
res2$beta_ASE
res2$D_hat

# A likelihood ratio test
LRT<-2*(res$LogMargProb-res2$LogMargProb)
print(round(pchisq(abs(LRT),df=1,lower.tail = F),3))


# sparsen the Hessian matrix
ptm<-proc.time()
res3 <- BivPPL(data,huge=TRUE)
proc.time() - ptm
res3$beta_hat
res3$beta_ASE
res3$D_hat

# Compare with coxme
# install.packages("coxme")
library(coxme)
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
res_coxme$coefficients 
sqrt(diag(vcov(res_coxme)))
assemble(as.vector(unlist(res_coxme$vcoef)))

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
res_coxme_n$coefficients 
sqrt(diag(vcov(res_coxme_n)))
assemble_n(as.vector(unlist(res_coxme_n$vcoef)))