#' @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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.