Type:	Package
Title:	Statistical Methods for Survival Data with Dependent Censoring
Version:	0.1.8
Maintainer:	Negera Wakgari Deresa <negera.deresa@gmail.com>
Description:	Several statistical methods for analyzing survival data under various forms of dependent censoring are implemented in the package. In addition to accounting for dependent censoring, it offers tools to adjust for unmeasured confounding factors. The implemented approaches allow users to estimate the dependency between survival time and dependent censoring time, based solely on observed survival data. For more details on the methods, refer to Deresa and Van Keilegom (2021) <doi:10.1093/biomet/asaa095>, Czado and Van Keilegom (2023) <doi:10.1093/biomet/asac067>, Crommen et al. (2024) <doi:10.1007/s11749-023-00903-9>, Deresa and Van Keilegom (2024) <doi:10.1080/01621459.2022.2161387>, Willems et al. (2025) <doi:10.48550/arXiv.2403.11860>, Ding and Van Keilegom (2025) and D'Haen et al. (2025) <doi:10.1007/s10985-025-09647-0>.
Imports:	survival, foreach, parallel, doParallel, pbivnorm, stats, MASS, nleqslv, methods, Matrix, EnvStats, mvtnorm, rafalib, rvinecopulib, matrixcalc, nloptr, numDeriv, copula, R6, lubridate, splines2
License:	GPL-3
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.3.2
Suggests:	testthat (≥ 3.0.0), rkriging, orthopolynom, OpenMx, stringr
Config/testthat/edition:	3
NeedsCompilation:	no
Packaged:	2025-11-06 15:46:53 UTC; Negera Deresa
Author:	Ilias Willems [aut], Gilles Crommen [aut], Negera Wakgari Deresa [aut, cre], Jie Ding [aut], Claudia Czado [aut], Myrthe D'Haen [aut], Ingrid Van Keilegom [aut]
Depends:	R (≥ 3.5.0)
Repository:	CRAN
Date/Publication:	2025-11-06 17:20:08 UTC

Compute bivariate survival probability

Description

This function calculates a bivariate survival probability based on multivariate normal distribution.

Usage

Bvprob(lx, ly, rho)

Arguments

Chronometer object

Description

R6 object that mimics a chronometer. It can be started, paused, record legs and stopped.

Methods

Public methods

• Chronometer$show()

• Chronometer$reset()

• Chronometer$start()

• Chronometer$stop()

• Chronometer$record.leg()

• Chronometer$get.chronometer.data()

• Chronometer$get.total.time()

• Chronometer$accumulate.legs()

• Chronometer$clone()

Method show()

Display the values stored in this chronometer object.

Usage

Chronometer$show()

Method reset()

Reset the chronometer.

Usage

Chronometer$reset()

Method start()

Start the chronometer

Usage

Chronometer$start()

Method stop()

Stop the chronometer

Usage

Chronometer$stop(leg.name = NULL)

Arguments

Method record.leg()

Record a leg time. The chronometer will continue running.

Usage

Chronometer$record.leg(leg.name = NULL)

Arguments

Method get.chronometer.data()

Like show method, but more rudimentary.

Usage

Chronometer$get.chronometer.data()

Method get.total.time()

Return the total time span between start and stop.

Usage

Chronometer$get.total.time(force = FALSE)

Arguments

Method accumulate.legs()

Return total time spent per leg category (using leg names).

Usage

Chronometer$accumulate.legs(force = FALSE)

Arguments

Method clone()

The objects of this class are cloneable with this method.

Usage

Chronometer$clone(deep = FALSE)

Arguments

Distance between vectors

Description

This function computes distance between two vectors based on L2-norm

Usage

Distance(b, a)

Arguments

Calculate the kernel function

Description

Calculate the kernel function

Usage

Kernel(u, name = "Gaussian")

Arguments

Calculate the likelihood function for the fully parametric joint distribution

Description

Calculate the likelihood function for the fully parametric joint distribution

Usage

Likelihood.Parametric(param, yobs, delta, copfam, margins, cure = FALSE)

Arguments

Calculate the profiled likelihood function with kernel smoothing

Description

Calculate the profiled likelihood function with kernel smoothing

Usage

Likelihood.Profile.Kernel(param, yobs, delta, copfam, margins, cure = FALSE)

Arguments

Solve the profiled likelihood function

Description

Solve the profiled likelihood function

Usage

Likelihood.Profile.Solve(
  yobs,
  delta,
  copfam,
  margins,
  ktau.init,
  parapar.init,
  cure,
  curerate.init,
  constraints,
  maxit,
  eps
)

Arguments

Calculate the semiparametric version of profiled likelihood function

Description

Calculate the semiparametric version of profiled likelihood function

Usage

Likelihood.Semiparametric(
  param,
  Syobs,
  yobs,
  delta,
  copfam,
  margins,
  cure = FALSE
)

Arguments

Change H to long format

Description

Change a nonparametric transformation function to long format

Usage

LongNPT(Z, T1, H)

Arguments

Fit a semiparametric transformation model for dependent censoring

Description

This function allows to estimate the dependency parameter along all other model parameters. First, estimates a non-parametric transformation function, and then at the second stage it estimates other model parameters assuming that the non-parametric function is known. The details for implementing the dependent censoring methodology can be found in Deresa and Van Keilegom (2021).

Usage

NonParTrans(
  resData,
  X,
  W,
  start = NULL,
  n.iter = 15,
  bootstrap = FALSE,
  n.boot = 50,
  eps = 0.001
)

Arguments

Value

This function returns a fit of a semiparametric transformation model; parameter estimates, estimate of the non-parametric transformation function, bootstrap standard errors for finite-dimensional parameters, the nonparametric cumulative hazard function, etc.

References

Deresa, N. and Van Keilegom, I. (2021). On semiparametric modelling, estimation and inference for survival data subject to dependent censoring, Biometrika, 108, 965–979.

Examples

# Toy data example to illustrate implementation
n = 300
beta = c(0.5, 1); eta = c(1,1.5); rho = 0.70
sigma = matrix(c(1,rho,rho,1),ncol=2)
err = MASS::mvrnorm(n, mu = c(0,0) , Sigma=sigma)
err1 = err[,1]; err2 = err[,2]
x1 = rbinom(n,1,0.5); x2 = runif(n,-1,1)
X = matrix(c(x1,x2),ncol=2,nrow=n); W = X   # data matrix
T1 = X%*%beta+err1
C =  W%*%eta+err2
T1 = exp(T1); C = exp(C)
A = runif(n,0,8); Y = pmin(T1,C,A)
d1 = as.numeric(Y==T1)
d2 = as.numeric(Y==C)
resData = data.frame("Z" = Y,"d1" = d1, "d2" = d2)   # should be data frame
colnames(X) = c("X1", "X2")
colnames(W) = c("W1","W2")

#  Bootstrap is false by default
output = NonParTrans(resData = resData, X = X, W = W, n.iter = 2)
output$parameterEstimates

Estimation of a parametric dependent censoring model without covariates.

Description

Note that it is not assumed that the association parameter of the copula function is known, unlike most other papers in the literature. The details for implementing the methodology can be found in Czado and Van Keilegom (2023).

Usage

ParamCop(Y, Delta, Copula, Dist.T, Dist.C, start = c(1, 1, 1, 1))

Arguments

Value

A table containing the minimized negative log-likelihood using the independence copula model, the estimated parameter values for the model with the independence copula, the minimized negative log-likelihood using the specified copula model and the estimated parameter values for the model with the specified copula.

References

Czado and Van Keilegom (2023). Dependent censoring based on parametric copulas. Biometrika, 110(3), 721-738.

Examples

tau = 0.75
Copula = "frank"
Dist.T = "weibull"
Dist.C = "lnorm"
par.T = c(2,1)
par.C = c(1,2)
n=1000

simdata<-TCsim(tau,Copula,Dist.T,Dist.C,par.T,par.C,n)

Y = simdata[[1]]
Delta = simdata[[2]]

ParamCop(Y,Delta,Copula,Dist.T,Dist.C)

Generate constraints of parameters

Description

Generate constraints of parameters

Usage

Parameters.Constraints(copfam, margins, cure)

Arguments

Estimate the model of D'Haen et al. (2025).

Description

This function estimates the parameters in the model of D'Haen et al. (2025).

Usage

QRdepCens(data, hp, var.estimate = FALSE, verbose = TRUE)

Arguments

Note

The variance estimation procedure could easily be paralelized. However, this is currently not implemented.

References

D'Haen, M., Van Keilegom, I. and Verhasselt, A. (2025). Quantile regression under dependent censoring with unknown association. Lifetime Data Analysis 31(2):253-299.

Examples

 # Load the data
 data(liver)

 # Give standard column names (required!)
 colnames(liver) <- c("patient", "Y", "Delta", "X1", "X2", "X3", "X4")
 liver <- liver[, c(-1, -6, -7)]

 # Run the model
 hp <- list(
   homoscedastic = FALSE,
   test_cop_name = "frank"
 )
 QRdepCens(liver, hp, var.estimate = FALSE)
 # Takes a while if var.estimate = TRUE...

Score equations of finite parameters

Description

This function computes the score vectors and the Jacobean matrix for finite model parameters.

Usage

ScoreEqn(theta, resData, X, W, H)

Arguments

Search function

Description

Function to indicate position of t in observed survival time

Usage

SearchIndicate(t, T1)

Arguments

Estimate a nonparametric transformation function

Description

This function estimates the nonparametric transformation function H when the survival time and censoring time are dependent given covariates. The estimating equation of H was derived based on the martingale ideas. More details about the derivation of a nonparmaetric estimator of H and its estimation algorithm can be found in Deresa and Van Keilegom (2021).

Usage

SolveH(theta, resData, X, W)

Arguments

Value

Returns the estimated transformation function H for a fixed value of parameters theta.

References

Deresa, N. and Van Keilegom, I. (2021). On semiparametric modelling, estimation and inference for survival data subject to dependent censoring, Biometrika, 108, 965–979.

Estimating equation for Ht1

Description

This function obtains an estimating equation of H at the first observed survival time t1.

Usage

SolveHt1(Ht1, Z, nu, t, X, W, theta)

Arguments

Cumulative hazard function of survival time under dependent censoring

Description

This function estimates the cumulative hazard function of survival time (T) under dependent censoring (C). The estimation makes use of the estimating equations derived based on martingale ideas.

Usage

SolveL(
  theta,
  resData,
  X,
  W,
  cop = c("Frank", "Gumbel", "Normal"),
  dist = c("Weibull", "lognormal")
)

Arguments

Value

This function returns an estimated hazard function, cumulative hazard function and distinct observed survival times;

Examples

n = 200
beta = c(0.5)
lambd = 0.35
eta = c(0.9,0.4)
X = cbind(rbinom(n,1,0.5))
W = cbind(rep(1,n),rbinom(n,1,0.5))
frank.cop <- copula::frankCopula(param = 5,dim = 2)
U = copula::rCopula(n,frank.cop)
T1 = (-log(1-U[,1]))/(lambd*exp(X*beta))         # Survival time'
T2 = (-log(1-U[,2]))^(1.1)*exp(W%*%eta)          # Censoring time
A = runif(n,0,15)                                # administrative censoring time
Z = pmin(T1,T2,A)
d1 = as.numeric(Z==T1)
d2 = as.numeric(Z==T2)
resData = data.frame("Z" = Z,"d1" = d1, "d2" = d2)
theta = c(0.3,1,0.3,1,2)

# Estimate cumulative hazard function
cumFit <- SolveL(theta, resData,X,W)
cumhaz = cumFit$cumhaz
time = cumFit$times

# plot hazard vs time

plot(time, cumhaz, type = "l",xlab = "Time",
ylab = "Estimated cumulative hazard function")

Cumulative hazard function of survival time under independent censoring

Description

This function estimates the cumulative hazard function of survival time (T) under the assumption of independent censoring. The estimating equation is derived based on martingale ideas.

Usage

SolveLI(theta, resData, X)

Arguments

Value

This function returns an estimated hazard function, cumulative hazard function and distinct observed survival times;

Examples

n = 200
beta = c(0.5)
lambd = 0.35
eta = c(0.9,0.4)
X = cbind(rbinom(n,1,0.5))
W = cbind(rep(1,n),rbinom(n,1,0.5))
frank.cop <- copula::frankCopula(param = 5,dim = 2)
U = copula::rCopula(n,frank.cop)
T1 = (-log(1-U[,1]))/(lambd*exp(X*beta))           # Survival time'
T2 = (-log(1-U[,2]))^(1.1)*exp(W%*%eta)            # Censoring time
A = runif(n,0,15)                                  # administrative censoring time
Z = pmin(T1,T2,A)
d1 = as.numeric(Z==T1)
d2 = as.numeric(Z==T2)
resData = data.frame("Z" = Z,"d1" = d1, "d2" = d2)
theta = c(0.3,1,0.3,1)

# Estimate cumulative hazard function

cumFit_ind <- SolveLI(theta, resData,X)

cumhaz = cumFit_ind$cumhaz
time = cumFit_ind$times

# plot hazard vs time

plot(time, cumhaz, type = "l",xlab = "Time",
ylab = "Estimated cumulative hazard function")

Estimate finite parameters based on score equations

Description

This function estimates the model parameters

Usage

SolveScore(theta, resData, X, W, H, eps = 0.001)

Arguments

Semiparametric Estimation of the Survival Function under Dependent Censoring

Description

Provide semiparametric approaches that can be used to model right-censored survival data under dependent censoring (without covariates). The copula-based approach is adopted and there is no need to explicitly specify the association parameter. One of the margins can be modeled nonparametrically. As a byproduct, both marginal distributions of survival and censoring times can be considered as fully parametric. The existence of a cured fraction concerning survival time can also be taken into consideration.

@references Czado and Van Keilegom (2023). Dependent censoring based on parametric copulas. Biometrika, 110(3), 721-738. @references Delhelle and Van Keilegom (2024). Copula based dependent censoring in cure models. TEST (to appear). @references Ding and Van Keilegom (2024). Semiparametric estimation of the survival function under dependent censoring (in preparation).

Usage

SurvDC(
  yobs,
  delta,
  tm = NULL,
  copfam = "frank",
  margins = list(survfam = NULL, censfam = "lnorm"),
  cure = FALSE,
  Var = list(do = TRUE, nboot = 200, level = 0.05),
  control = list(maxit = 300, eps = 1e-06, trace = TRUE, ktau.inits = NULL)
)

Arguments

Details

This unified function provide approaches that can be used to model right-censored survival data under dependent censoring (without covariates). Various specifications of marginal distributions can be considered by choosing different combinations of the provided arguments. Generally speaking, the following two scenarios are what we mainly focused on:

nonparametric survival margin and parametric censoring margin (without cure)

survfam = NULL, censfam is not NULL and cure = FALSE.

nonparametric survival margin and parametric censoring margin (with cure)

survfam = NULL, censfam is not NULL and cure = TRUE.

As byproducts, several other scenarios (the distribution of the underlying survival time is not nonparametric but fully parametric) can also be considered by this R function:

parametric survival and censoring margins (without cure)

both survfam and censfam are not NULL and cure = FALSE.

parametric survival and censoring margins (with cure)

both survfam and censfam are not NULL and cure = TRUE.

parametric survival margin and nonparametric censoring margin (without cure)

survfam is not NULL, censfam = NULL and cure = FALSE.

Furthermore, one might expect that a scenario with "parametric survival margin and nonparametric censoring margin (with cure)" can also be included. Indeed, it can be done based on: survfam is not NULL, censfam = NULL and cure = TRUE. However, from a theoretical perspective of view, whether this type of modeling is reasonable or not still needs further investigations.

We emphasize that the first scenario (in byproducts) has also be considered in another function of this package. Specifically, the scenario of "parametric survival margin and nonparametric censoring margin (without cure)" can be fitted based on ParamCop(). However, the default joint modeling of survival and censoring times are based on their joint survival function in line with the semiparametric case (instead of modeling joint distribution function directly as in Czado and Van Keilegom (2023) <doi:10.1093/biomet/asac067>), but the idea of estimation methodology are exactly the same.

Value

A list of fitted results is returned. Within this outputted list, the following elements can be found:

probs

survival probabilities of the survial margin at tm.

ktau

Kendall's tau.

parapar

estimation of all parameters (except Kendall's tau) contained in the parametric part.

GoF

goodness-of-test results.

curerate

cure rate. If cure = FALSE, it is NULL.

Examples

#----------------------------------------------------------#
# Basic preparations before running subsequent examples ####
#----------------------------------------------------------#

# library necessary packages

#------------------------------------------------------------------------#
# simulated data from Frank copula log-Normal margins (without cure)
#------------------------------------------------------------------------#

# generate the simulated data

# - the sample size of the generated data
n <- 1000

# information on the used copula
copfam.true <- "frank"
ktau.true <- 0.5
coppar.true <- 5.74

# parameters of the underlying log-normal marginal distributions
survpar.true <- c(2.20,1.00)
censpar.true <- c(2.20,0.25)

# - true underlying survival and censoring times
set.seed(1)
u.TC <- copula::rCopula(
  n        = n,
  copula   = copula::archmCopula(
    family = copfam.true,
    param  = coppar.true,
    dim    = 2
  )
)
yobs.T <- qlnorm(1-u.TC[,1],survpar.true[1],survpar.true[2])
yobs.C <- qlnorm(1-u.TC[,2],censpar.true[1],censpar.true[2])

# observations
yobs  <- pmin(yobs.T,yobs.C)
delta <- as.numeric(yobs.T<=yobs.C)
cat("censoring rate is", mean(1-delta))

# model the data under different scenarios

# scenario 1: nonparametric survival margin and parametric censoring margin
set.seed(1)
sol.scenario1 <- SurvDC(
  yobs    = yobs,
  delta   = delta,
  tm      = quantile(yobs, c(0.25,0.50,0.75)),
  copfam  = copfam.true,
  margins = list(survfam = NULL, censfam = "lnorm"),
  Var     = list(do = FALSE, nboot = 50)
)
sol.scenario1$probs
sol.scenario1$ktau
sol.scenario1$parapar

# scenario 2: parametric survival and censoring margins
set.seed(1)
sol.scenario2 <- SurvDC(
  yobs    = yobs,
  delta   = delta,
  tm      = quantile(yobs, c(0.25,0.50,0.75)),
  copfam  = copfam.true,
  margins = list(survfam = "lnorm", censfam = "lnorm"),
  Var     = list(do = FALSE, nboot = 50)
)
sol.scenario2$probs
sol.scenario2$ktau
sol.scenario2$parapar

# scenario 3: parametric survival margin and nonparametric censoring margin
set.seed(1)
sol.scenario3 <- SurvDC(
  yobs    = yobs,
  delta   = delta,
  tm      = quantile(yobs, c(0.25,0.50,0.75)),
  copfam  = copfam.true,
  margins = list(survfam = "lnorm", censfam = NULL),
  Var     = list(do = FALSE, nboot = 50)
)
sol.scenario3$probs
sol.scenario3$ktau
sol.scenario3$parapar

#------------------------------------------------------------------------
# simulated data from Frank copula log-Normal margins (with cure)
#------------------------------------------------------------------------

# generate the simulated data

#  true underlying cure rate
curerate.true <- 0.2

# true underlying survival and censoring times
set.seed(1)
u.TC <- copula::rCopula(
  n        = n,
  copula   = copula::archmCopula(
    family = copfam.true,
    param  = coppar.true,
    dim    = 2
  )
)
yobs.T <- sapply(u.TC[,1],function(uT){
  if(uT<=curerate.true){ val <- Inf }else{
    val <- EnvStats::qlnormTrunc((1-uT)/(1-curerate.true),survpar.true[1],survpar.true[2],0,15)
  }
  return(val)
})
yobs.C <- qlnorm(1-u.TC[,2],censpar.true[1],censpar.true[2])
cat("cure rate is",mean(yobs.T==Inf))

#  observations
yobs  <- pmin(yobs.T,yobs.C)
delta <- as.numeric(yobs.T<=yobs.C)
cat("censoring rate is",mean(1-delta))

# model the data under different scenarios (with cure)

# scenario 4: parametric survival and censoring margins
set.seed(1)
sol.scenario4 <- SurvDC(
  yobs    = yobs,
  delta   = delta,
  tm      = quantile(yobs, c(0.25,0.50,0.75)),
  copfam  = copfam.true,
  margins = list(survfam = "lnorm", censfam = "lnorm"),
  Var     = list(do = FALSE, nboot = 50),
  cure    = TRUE
)
sol.scenario4$probs
sol.scenario4$ktau
sol.scenario4$parapar
sol.scenario4$curerate

# scenario 5: nonparametric survival margin and parametric censoring margin
set.seed(1)
sol.scenario5 <- SurvDC(
  yobs    = yobs,
  delta   = delta,
  tm      = quantile(yobs, c(0.25,0.50,0.75)),
  copfam  = copfam.true,
  margins = list(survfam = NULL, censfam = "lnorm"),
  Var     = list(do = FALSE, nboot = 50),
  cure    = TRUE
)
sol.scenario5$probs
sol.scenario5$ktau
sol.scenario5$parapar
sol.scenario5$curerate

Calculate the goodness-of-fit test statistic

Description

Calculate the goodness-of-fit test statistic

Usage

SurvDC.GoF(
  yobs,
  delta,
  copfam,
  margins,
  ktau,
  parapar,
  cure = FALSE,
  curerate = NULL
)

Arguments

Estimated survival function based on copula-graphic estimator (Archimedean copula only)

Description

Estimated survival function based on copula-graphic estimator (Archimedean copula only)

Usage

SurvFunc.CG(tm = NULL, yobs, delta, copfam, ktau, coppar = NULL)

Arguments

Estimated survival function based on Kaplan-Meier estimator

Description

Estimated survival function based on Kaplan-Meier estimator

Usage

SurvFunc.KM(tm = NULL, yobs, delta, type = "right")

Arguments

Maximum likelihood estimator for a given parametric distribution

Description

Maximum likelihood estimator for a given parametric distribution

Usage

SurvMLE(
  yobs,
  delta,
  distribution,
  truncation = NULL,
  cure = FALSE,
  maxit = 300
)

Arguments

Likelihood for a given parametric distribution

Description

Likelihood for a given parametric distribution

Usage

SurvMLE.Likelihood(
  param,
  yobs,
  delta,
  distribution,
  truncation = NULL,
  cure = FALSE
)

Arguments

Function to simulate (Y,Delta) from the copula based model for (T,C).

Description

Generates the follow-up time and censoring indicator according to the specified model.

Usage

TCsim(
  tau = 0,
  Copula = "frank",
  Dist.T = "lnorm",
  Dist.C = "lnorm",
  par.T = c(0, 1),
  par.C = c(0, 1),
  n = 10000
)

Arguments

Value

A list containing the generated follow-up times and censoring indicators.

Examples

tau = 0.5
Copula = "gaussian"
Dist.T = "lnorm"
Dist.C = "lnorm"
par.T = c(1,1)
par.C = c(2,2)
n=1000

simdata <- TCsim(tau,Copula,Dist.T,Dist.C,par.T,par.C,n)
Y = simdata[[1]]
Delta = simdata[[2]]
hist(Y)
mean(Delta)

Nonparametric bootstrap approach for a Semiparametric transformation model under dependent censpring

Description

This function estimates the bootstrap standard errors for the finite-dimensional model parameters and for the non-parametric transformation function. Parallel computing using foreach has been used to speed up the estimation of standard errors.

Usage

boot.nonparTrans(init, resData, X, W, n.boot, n.iter, eps)

Arguments

Value

Bootstrap standard errors for parameter estimates and for estimated cumulative hazard function.

Prepare initial values within the control arguments

Description

Prepare initial values within the control arguments

Usage

control.arguments(maxit = 300, eps = 1e-06, trace = TRUE, ktau.inits = NULL)

Arguments

The distribution function of the Archimedean copula

Description

The distribution function of the Archimedean copula

Usage

copdist.Archimedean(x, copfam, ktau, coppar = NULL)

Arguments

The h-function of the copula

Description

The h-function of the copula

Usage

cophfunc(x, coppar, copfam, condvar = 1)

Arguments

Convert the copula parameter the Kendall's tau

Description

Convert the copula parameter the Kendall's tau

Usage

coppar.to.ktau(coppar, copfam)

Arguments

Estimate the competing risks model of Willems et al. (2025).

Description

This function estimates the parameters in the competing risks model described in Willems et al. (2025). Note that this model extends the model of Crommen, Beyhum and Van Keilegom (2024) and as such, this function also implements their methodology.

Usage

estimate.cmprsk(
  data,
  admin,
  conf,
  eoi.indicator.names = NULL,
  Zbin = NULL,
  inst = "cf",
  realV = NULL,
  compute.var = TRUE,
  eps = 0.001
)

Arguments

Value

A list of parameter estimates in the second stage of the estimation algorithm (hence omitting the estimates for the control function), as well as an estimate of their variance and confidence intervals.

References

Willems et al. (2025). Flexible control function approach under competing risks (submitted).

Crommen, G., Beyhum, J., and Van Keilegom, I. (2024). An instrumental variable approach under dependent censoring. Test, 33(2), 473-495.

Examples

n <- 200

# Set parameters
gamma <- c(1, 2, 1.5, -1)
theta <- c(0.5, 1.5)
eta1 <- c(1, -1, 2, -1.5, 0.5)
eta2 <- c(0.5, 1, 1, 3, 0)

# Generate exogenous covariates
x0 <- rep(1, n)
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)

# Generate confounder and instrument
w <- rnorm(n)
V <- rnorm(n, 0, 2)
z <- cbind(x0, x1, x2, w) %*% gamma + V
realV <- z - (cbind(x0, x1, x2, w) %*% gamma)

# Generate event times
err <- MASS::mvrnorm(n, mu = c(0, 0), Sigma =
matrix(c(3, 1, 1, 2), nrow = 2, byrow = TRUE))
bn <- cbind(x0, x1, x2, z, realV) %*% cbind(eta1, eta2) + err
Lambda_T1 <- bn[,1]; Lambda_T2 <- bn[,2]
x.ind = (Lambda_T1>0)
y.ind <- (Lambda_T2>0)
T1 <- rep(0,length(Lambda_T1))
T2 <- rep(0,length(Lambda_T2))
T1[x.ind] = ((theta[1]*Lambda_T1[x.ind]+1)^(1/theta[1])-1)
T1[!x.ind] = 1-(1-(2-theta[1])*Lambda_T1[!x.ind])^(1/(2-theta[1]))
T2[y.ind] = ((theta[2]*Lambda_T2[y.ind]+1)^(1/theta[2])-1)
T2[!y.ind] = 1-(1-(2-theta[2])*Lambda_T2[!y.ind])^(1/(2-theta[2]))
# Generate adminstrative censoring time
C <- runif(n, 0, 40)

# Create observed data set
y <- pmin(T1, T2, C)
delta1 <- as.numeric(T1 == y)
delta2 <- as.numeric(T2 == y)
da <- as.numeric(C == y)
data <- data.frame(cbind(y, delta1, delta2, da, x0, x1, x2, z, w))
colnames(data) <- c("y", "delta1", "delta2", "da", "x0", "x1", "x2", "z", "w")

# Estimate the model
admin <- TRUE                # There is administrative censoring in the data.
conf <- TRUE                 # There is confounding in the data (z)
eoi.indicator.names <- NULL  # We will not impose that T1 and T2 are independent
Zbin <- FALSE                # The confounding variable z is not binary
inst <- "cf"                 # Use the control function approach
compute.var <- TRUE          # Variance of estimates should be computed.
# Since we don't use the oracle estimator, this argument is ignored anyway
realV <- NULL
estimate.cmprsk(data, admin, conf, eoi.indicator.names, Zbin, inst, realV,
                compute.var)

Fit Dependent Censoring Models

Description

This function allows to estimate the dependency parameter along all other model parameters. First, estimates the cumulative hazard function, and then at the second stage it estimates other model parameters assuming that the cumulative hazard function is known. The details for implementing the dependent censoring methodology can be found in Deresa and Van Keilegom (2024).

Usage

fitDepCens(
  resData,
  X,
  W,
  cop = c("Frank", "Gumbel", "Normal"),
  dist = c("Weibull", "lognormal"),
  start = NULL,
  n.iter = 50,
  bootstrap = TRUE,
  n.boot = 150,
  ncore = 7,
  eps = 1e-04
)

Arguments

Value

This function returns a fit of dependent censoring model; parameter estimates, estimate of the cumulative hazard function, bootstrap standard errors for finite-dimensional parameters, the nonparametric cumulative hazard function, etc.

References

Deresa and Van Keilegom (2024). Copula based Cox proportional hazards models for dependent censoring, Journal of the American Statistical Association, 119:1044-1054.

Examples

# Toy data example to illustrate implementation
n = 300
beta = c(0.5)
lambd = 0.35
eta = c(0.9,0.4)
X = cbind(rbinom(n,1,0.5))
W = cbind(rep(1,n),rbinom(n,1,0.5))
# generate dependency structure from Frank
frank.cop <- copula::frankCopula(param = 5,dim = 2)
U = copula::rCopula(n,frank.cop)
T1 = (-log(1-U[,1]))/(lambd*exp(X*beta))           # Survival time
T2 = (-log(1-U[,2]))^(1.1)*exp(W%*%eta)            # Censoring time
A = runif(n,0,15)                                  # administrative censoring time
Z = pmin(T1,T2,A)
d1 = as.numeric(Z==T1)
d2 = as.numeric(Z==T2)
resData = data.frame("Z" = Z,"d1" = d1, "d2" = d2)   # should be data frame
colnames(W) <- c("ones","cov1")
colnames(X) <- "cov.surv"

# Fit dependent censoring model

fit <- fitDepCens(resData = resData, X = X, W = W, bootstrap = FALSE)

# parameter estimates

fit$parameterEstimates

# summary fit results
summary(fit)

# plot cumulative hazard vs time

plot(fit$observedTime, fit$cumhazardFunction, type = "l",xlab = "Time",
ylab = "Estimated cumulative hazard function")

Fit Independent Censoring Models

Description

This function allows to estimate all model parameters under the assumption of independent censoring. First, estimates the cumulative hazard function, and then at the second stage it estimates other model parameters assuming that the cumulative hazard is known.

Usage

fitIndepCens(
  resData,
  X,
  W,
  dist = c("Weibull", "lognormal"),
  start = NULL,
  n.iter = 50,
  bootstrap = TRUE,
  n.boot = 150,
  ncore = 7,
  eps = 1e-04
)

Arguments

Value

This function returns a fit of independent censoring model; parameter estimates, estimate of the cumulative hazard function, bootstrap standard errors for finite-dimensional parameters, the nonparametric cumulative hazard function, etc.

Examples

# Toy data example to illustrate implementation
n = 300
beta = c(0.5)
lambd = 0.35
eta = c(0.9,0.4)
X = cbind(rbinom(n,1,0.5))
W = cbind(rep(1,n),rbinom(n,1,0.5))
# generate dependency structure from Frank
frank.cop <- copula::frankCopula(param = 5,dim = 2)
U = copula::rCopula(n,frank.cop)
T1 = (-log(1-U[,1]))/(lambd*exp(X*beta))                  # Survival time'
T2 = (-log(1-U[,2]))^(1.1)*exp(W%*%eta)                   # Censoring time
A = runif(n,0,15)                                         # administrative censoring time
Z = pmin(T1,T2,A)
d1 = as.numeric(Z==T1)
d2 = as.numeric(Z==T2)
resData = data.frame("Z" = Z,"d1" = d1, "d2" = d2)     # should be data frame

colnames(W) <- c("ones","cov1")
colnames(X) <- "cov.surv"

# Fit independent censoring model

fitI <- fitIndepCens(resData = resData, X = X, W = W, bootstrap = FALSE)

# parameter estimates

fitI$parameterEstimates

# summary fit results
summary(fitI)

# plot cumulative hazard vs time

 plot(fitI$observedTime, fitI$cumhazardFunction, type = "l",xlab = "Time",
 ylab = "Estimated cumulative hazard function")

The generator function of the Archimedean copula

Description

The generator function of the Archimedean copula

Usage

generator.Archimedean(x, coppar, copfam, inverse = FALSE)

Arguments

Convert the Kendall's tau into the copula parameter

Description

Convert the Kendall's tau into the copula parameter

Usage

ktau.to.coppar(ktau, copfam)

Arguments

Liver cirrhosis data set.

Description

End stage liver disease data set, as for example analyzed by D'Haen et al. (2025).

Usage

liver

Format

A data frame with 281 rows and 7 variables:

patient

patient ID.

time

Survival time.

status

Censoring indicator.

age

Age of patient.

gender

Gender of patient (1 - male, 0 - female).

bmi

Body mass index.

ukeld

UK End-stage Liver Disease score. Higher scores indicate more severe disease state.

Log-likelihood function for the Clayton copula.

Description

This likelihood function is maximized to estimate the model parameters under the Clayton copula.

Usage

loglike.clayton.unconstrained(para, Y, Delta, Dist.T, Dist.C)

Arguments

Value

Maximized log-likelihood value.

Log-likelihood function for the Frank copula.

Description

This likelihood function is maximized to estimate the model parameters under the Frank copula.

Usage

loglike.frank.unconstrained(para, Y, Delta, Dist.T, Dist.C)

Arguments

Value

Maximized log-likelihood value.

Log-likelihood function for the Gaussian copula.

Description

This likelihood function is maximized to estimate the model parameters under the Gaussian copula.

Usage

loglike.gaussian.unconstrained(para, Y, Delta, Dist.T, Dist.C)

Arguments

Value

Maximized log-likelihood value.

Log-likelihood function for the Gumbel copula.

Description

This likelihood function is maximized to estimate the model parameters under the Gumbel copula.

Usage

loglike.gumbel.unconstrained(para, Y, Delta, Dist.T, Dist.C)

Arguments

Value

Maximized log-likelihood value.

Log-likelihood function for the independence copula.

Description

This likelihood function is maximized to estimate the model parameters under the independence copula.

Usage

loglike.indep.unconstrained(para, Y, Delta, Dist.T, Dist.C)

Arguments

Value

Maximized log-likelihood value.

Fit the dependent censoring models.

Description

Estimates the model parameters by maximizing the log-likelihood.

Usage

optimlikelihood(Y, Delta, Copula, Dist.T, Dist.C, start)

Arguments

Value

A list containing the minimized negative log-likelihood using the independence copula model, the estimated parameter values for the model with the independence copula, the minimized negative log-likelihood using the specified copula model and the estimated parameter values for the model with the specified copula.

Obtain the value of the density function

Description

Obtain the value of the density function

Usage

parafam.d(x, parameter, distribution, truncation = NULL)

Arguments

Obtain the value of the distribution function

Description

Obtain the value of the distribution function

Usage

parafam.p(x, parameter, distribution, truncation = NULL)

Arguments

Obtain the adjustment value of truncation

Description

Obtain the adjustment value of truncation

Usage

parafam.trunc(truncation, parameter, distribution)

Arguments

Estimate the model of Willems et al. (2025).

Description

This function estimates bounds on the coefficients the single- index model \Lambda(x^\top \beta(t)) for the conditional cumulative distribution function of the event time.

Usage

pi.surv(
  data,
  idx.param.of.interest,
  idxs.c,
  t,
  par.space,
  search.method = "GS",
  add.options = list(),
  verbose = 0,
  picturose = FALSE,
  parallel = FALSE
)

Arguments

Value

Matrix containing the identified intervals of the specified coefficients, as well as corresponding convergence information of the estimation algorithm.

References

Willems, I., Beyhum, J. and Van Keilegom, I. (2025). Partial identification for a class of survival models under dependent censoring. (Submitted).

Examples

  # Clear workspace
  rm(list = ls())

  # Load the survival package
  library(survival)

  # Set random seed
  set.seed(123)

  # Load and preprocess data
  data <- survival::lung
  data[, "intercept"] <- rep(1, nrow(data))
  data[, "status"] <- data[, "status"] - 1
  data <- data[, c("time", "status", "intercept", "age", "sex")]
  colnames(data) <- c("Y", "Delta", "X0", "X1", "X2")

  # Standardize age variable
  data[, "X1"] <- scale(data[, "X1"])

  ## Example:
  ## - Link function: AFT link function (default setting)
  ## - Number of IF: 5 IF per continuous covariate (default setting)
  ## - Search method: Binary search
  ## - Type of IF: Cubic spline functions for continuous covariate, indicator
  ##   function for discrete covariate (default setting).

  # Settings for main estimation function
  idx.param.of.interest <- 2 # Interest in effect of age
  idxs.c <- 1                # X1 (age) is continuous
  t <- 200                   # Model imposed at t = 200
  search.method <- "GS"      # Use binary search
  par.space <- matrix(rep(c(-10, 10), 3), nrow = 3, byrow = TRUE)
  add.options <- list()
  picturose <- TRUE
  parallel <- FALSE

  # Estimate the identified intervals
  pi.surv(data, idx.param.of.interest, idxs.c, t, par.space,
          search.method = search.method, add.options = add.options,
          picturose = picturose, parallel = parallel)

Summary of depCensoringFit object

Description

Summary of depCensoringFit object

Usage

## S3 method for class 'depFit'
summary(object, ...)

Arguments

Value

Summary of dependent censoring model fit in the form of table

Summary of indepCensoringFit object

Description

Summary of indepCensoringFit object

Usage

## S3 method for class 'indepFit'
summary(object, ...)

Arguments

Value

Summary of independent censoring model fit in the form of table