R/surv_data_dc.R

Defines functions surv_data_dc

Documented in surv_data_dc

#' @title Generate a sample of time to event dataset with dependent right censoring under an Archimedean copula
#'
#' @description Generate a sample of time to event dataset with, dependent right censoring based on one of the Archimedean copulas the given Kendall's tau, sample size n and covariates matrix Z.
#'
#' @usage surv_data_dc(n, a, Z, lambda, betas, phis, cons7, cons9, tau, copula, distr.ev, distr.ce)
#'
#' @param n the sample size, or the number of the subjects in a sample.
#' @param a the shape parameter of baseline hazard for the event time \eqn{T}.
#' @param Z the covariate matrix with dimension of \eqn{n} by \eqn{p}, where \eqn{p} is the number of covariates.
#' @param lambda the scale parameter of baseline hazard for event time \eqn{T}.
#' @param betas the regression coefficient vector of proportional hazard model for the event time \eqn{T} with dimenion of \eqn{p} by \eqn{1}.
#' @param phis the regression coefficient vector of proportional hazard model for dependent censoring time \eqn{C} with dimenion of \eqn{p} by \eqn{1}.
#' @param cons7 the parameter of baseline hazard for the dependent censoring time \eqn{C} if assuming an exponential distribution.
#' @param cons9 the upper limit parameter of uniform distribution for the independent censoring time \eqn{A},
#' i.e. \eqn{A}~U(0, cons9).
#' @param tau the Kendall's correlation coefficient between \eqn{T} and \eqn{C}.
#' @param copula the Archemedean copula that captures the dependence between \eqn{T} and \eqn{C}, a characteristc value,
#' i.e. 'independent', 'clayton', 'gumbel' or 'frank'.
#' @param distr.ev the distribution of the event time, a characteristc value, i.e. 'weibull' or 'log logit'.
#' @param distr.ce the distribution of the dependent censoring time, a characteristc value, i.e. 'exponential' or 'weibull'.
#'
#' @details surv_data_dc allows to generate a survival dataset under dependent right censoring, at sample size \code{n}, based on one of the Archimedean \code{copula},
#' Kendall's \code{tau}, and covariates matrix \code{Z} with dimension of \eqn{n} by \eqn{p}. For example, at \code{p=2}, we have \code{Z=cbind(Z1, Z2)},
#' where \code{Z1} is treatment generated by distribution of bernoulli(0.5), i.e. 1 represents treatment group and 0 represents control group; \code{Z2} is the age
#' generated by distribution of U(-10, 10).
#'
#' The generated dataset includes three varaibles, which are \eqn{X_i}, \eqn{\delta_i} and \eqn{\eta_i}, i.e.
#' \eqn{X_i=min(T_i, C_i, A_i)}, \eqn{\delta_i=I(X_i=T_i)} and \eqn{\eta_i=I(X_i=C_i)}, for \eqn{i=1,\ldots,n}.
#' 'T' represents the event time, whose hazard function is \deqn{h_T(x)=h_{0T}(x)exp(Z^{\top}\beta)},
#' where the baseline hazard can take weibull form, i.e. \eqn{h_{0T}(x) = ax^{a-1} / \lambda^a}, or log logistic form, i.e.
#' \deqn{ h_{0T}(x) = \frac{ \frac{ 1 }{ a exp( \lambda ) } ( \frac{ x }{ exp( \lambda ) } )^{1/a -1 } }{ 1 + ( \frac{ x }{ exp( \lambda ) } )^{1/a}  } }.
#' 'C' represents the dependent censoring time, whose hazard function is \eqn{ h_{C}(x) = h_{0C}(x)exp( Z^{\top}\phi) }, where the baseline hazard can take exponential form, i.e.
#'  \eqn{h_{0C}(x)=cons7}, or weibull form, i.e. \eqn{h_{0C}(x) = ax^{a-1} / \lambda^a}.'A' represents the administrative or independent censoring time, where A~U(0, cons9).
#'
#' @return A sample of time to event dataset under dependent right censoring, which includes observed time \eqn{X}, event indicator \eqn{\delta} and dependent censoring indicator \eqn{\eta}.
#'
#' @author Jing Xu, Jun Ma, Thomas Fung
#'
#' @references Xu J, Ma J, Connors MH, Brodaty H. (2018). \emph{"Proportional hazard model estimation under dependent censoring using copulas and penalized likelihood"}.
#' Statistics in Medicine 37, 2238–2251.
#'
#' @seealso \code{\link{coxph_mpl_dc}}
#'
#' @importFrom copula indepCopula claytonCopula gumbelCopula frankCopula iTau rCopula
#' @importFrom stats pnorm qexp qweibull runif
#'
#' @examples
#'  ##-- Copula types
#'  copula3 <- 'frank'
#'
#'  ##-- Marginal distribution for T, C, and A
#'  a <- 2
#'  lambda <- 2
#'  cons7 <- 0.2
#'  cons9 <- 10
#'  tau <- 0.8
#'  betas <- c(-0.5, 0.1)
#'  phis <- c(0.3, 0.2)
#'  distr.ev <- 'weibull'
#'  distr.ce <- 'exponential'
#'
#'  ##-- Sample size
#'  n <- 200
#'
#'  ##-- One sample Monte Carlo dataset
#'  cova <- cbind(rbinom(n, 1, 0.5), runif(n, min=-10, max=10))
#'  surv <- surv_data_dc(n, a, cova, lambda, betas, phis, cons7, cons9,
#'                      tau, copula3, distr.ev, distr.ce)
#'  n <- nrow(cova)
#'  p <- ncol(cova)
#'  ##-- event and dependent censoring proportions
#'  colSums(surv)[c(2,3)]/n
#'  X <- surv[,1] # Observed time
#'  del<-surv[,2] # failure status
#'  eta<-surv[,3] # dependent censoring status
#'
#'
#' @export

surv_data_dc<-function(n, a, Z, lambda, betas, phis, cons7, cons9, tau, copula, distr.ev, distr.ce)
{
  {
    if(copula=='independent')
    {
      cop<-indepCopula(2)
    }
    else if(copula=='clayton')
    {
      alpha<-iTau(claytonCopula(100), tau)
      cop <- claytonCopula(alpha, dim = 2)
    }
    else if(copula=='gumbel')
    {
      alpha<-iTau(gumbelCopula(100), tau)
      cop <- gumbelCopula(alpha, dim = 2)
    }
    else if(copula=='frank')
    {
      alpha<-iTau(frankCopula(100), tau)
      cop <- frankCopula(alpha, dim = 2)
    }
    else (stop("Error! Define a proper copula function"))
  }
  x<-rep(0, n)
  del<-rep(0, n)
  eta<-rep(0, n)
  expbetT <- exp(Z%*%betas)
  lamT<-lambda*(exp(-Z%*%betas))^(1/a)
  lamCe<-cons7*exp(Z%*%phis)
  lamCw<-lambda*(exp(-Z%*%phis))^(1/a)
  for (i in (1:n))
  {
    cuv<-rCopula(1, cop)      #generate the values of marginal suvival functions
    if( distr.ev == 'weibull' )
    {
      ti<-qweibull((1-cuv[1,1]), a, lamT[i])   #weibull distribution
    }
    else if ( distr.ev == 'log logit' )
    {
      ti <- exp( lambda )*( cuv[1,1]^( -1/expbetT[i] ) - 1 )^( a )     #log logit distribution
    }
    else ( stop("Error! Define a proper distribution type for the event time") )
    if( distr.ce=='exponential' )
    {
      ci<-qexp((1-cuv[1,2]), lamCe[i]) #exponential distribution
    }
    else if ( distr.ce=='weibull' ) {
      ci<-qweibull((1-cuv[1,2]), a, lamCw[i])   #weibull distribution
    }
    else ( stop("Error! Define a proper distribution type for the dependent censoring time") )
    ai<-runif(1, min=0, max=cons9)
    if(ti<ci & ti<ai){
      del[i]<-1
      x[i]<-ti
    }
    else if (ci<ti & ci <ai)
    {
      eta[i]<-1
      x[i]<-ci
    }
    else{
      del[i]<-0
      eta[i]<-0
      x[i]<-ai
    }
  }
  dat<-cbind(x, del, eta)
  return(dat)
}

Try the survivalMPLdc package in your browser

Any scripts or data that you put into this service are public.

survivalMPLdc documentation built on Oct. 24, 2020, 5:06 p.m.