,
Gilles Crommen [aut] ,
Negera Wakgari Deresa [aut, cre]
,
Jie Ding [aut] ,
Claudia Czado [aut] ,
Ingrid Van Keilegom [aut] | Title: | Statistical Methods for Survival Data with Dependent Censoring |
|---|---|
| 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>, Rutten et al. (2024+) <doi:10.48550/arXiv.2403.11860> and Ding and Van Keilegom (2024). |
| Authors: | Ilias Willems [aut] ,
Gilles Crommen [aut] ,
Negera Wakgari Deresa [aut, cre]
,
Jie Ding [aut] ,
Claudia Czado [aut] ,
Ingrid Van Keilegom [aut] |
| Maintainer: | Negera Wakgari Deresa <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1.7 |
| Built: | 2025-03-06 09:34:09 UTC |
| Source: | https://github.com/nago2020/depcensoring |
This function performs the approximation step in the EAM
algorithm. More specifically, it fits a Gaussian-process regression model
(Kriging) to the evaluated data points .
A_step(evaluations, verbose = 0)A_step(evaluations, verbose = 0)
evaluations |
Matrix containing each point that was already evaluated, alongside the corresponding test statistic and critical value, as its rows. |
verbose |
Verbosity parameter. |
Results of the A-step.
Package rkriging.
This function estimates the bootstrap standard errors for the finite-dimensional model parameters and for the non-parametric cumulative hazard function. Parallel computing using foreach has been used to speed up the estimation of standard errors.
boot.fun( init, resData, X, W, lhat, cumL, dist, k, lb, ub, Obs.time, cop, n.boot, n.iter, ncore, eps )boot.fun( init, resData, X, W, lhat, cumL, dist, k, lb, ub, Obs.time, cop, n.boot, n.iter, ncore, eps )
init |
Initial values for the finite dimensional parameters obtained from the fit of |
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T |
W |
Data matrix with covariates related to C. First column of W should be a vector of ones |
lhat |
Initial values for the hazard function obtained from the fit of |
cumL |
Initial values for the cumulative hazard function obtained from the fit of |
dist |
The distribution to be used for the dependent censoring time C. Only two distributions are allowed, i.e, Weibull
and lognormal distributions. With the value |
k |
Dimension of X |
lb |
lower boundary for finite dimensional parameters |
ub |
Upper boundary for finite dimensional parameters |
Obs.time |
Observed survival time, which is the minimum of T, C and A, where A is the administrative censoring time. |
cop |
Which copula should be computed to account for dependency between T and C. This argument can take
one of the values from |
n.boot |
Number of bootstraps to use in the estimation of bootstrap standard errors. |
n.iter |
Number of iterations; the default is |
ncore |
The number of cores to use for parallel computation is configurable, with the default |
eps |
Convergence error. This is set by the user in such away that the desired convergence is met; the default is |
Bootstrap standard errors for parameter estimates and for estimated cumulative hazard function.
This function estimates the bootstrap standard errors for the finite-dimensional model parameters and for the non-parametric cumulative hazard function under the assumption of independent censoring. Parallel computing using foreach has been used to speed up the computation.
boot.funI( init, resData, X, W, lhat, cumL, dist, k, lb, ub, Obs.time, n.boot, n.iter, ncore, eps )boot.funI( init, resData, X, W, lhat, cumL, dist, k, lb, ub, Obs.time, n.boot, n.iter, ncore, eps )
init |
Initial values for the finite dimensional parameters obtained from the fit of |
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T |
W |
Data matrix with covariates related to C. First column of W should be a vector of ones |
lhat |
Initial values for the hazard function obtained from the fit of |
cumL |
Initial values for the cumulative hazard function obtained from the fit of |
dist |
The distribution to be used for the dependent censoring time C. Only two distributions are allowed, i.e, Weibull
and lognormal distributions. With the value |
k |
Dimension of X |
lb |
lower boundary for finite dimensional parameters |
ub |
Upper boundary for finite dimensional parameters |
Obs.time |
Observed survival time, which is the minimum of T, C and A, where A is the administrative censoring time. |
n.boot |
Number of bootstraps to use in the estimation of bootstrap standard errors. |
n.iter |
Number of iterations; the default is |
ncore |
The number of cores to use for parallel computation is configurable, with the default |
eps |
Convergence error. This is set by the user in such away that the desired convergence is met; the default is |
Bootstrap standard errors for parameter estimates and for estimated cumulative hazard function.
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.
boot.nonparTrans(init, resData, X, W, n.boot, n.iter, eps)boot.nonparTrans(init, resData, X, W, n.boot, n.iter, eps)
init |
Initial values for the finite dimensional parameters obtained from the fit of |
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T |
W |
Data matrix with covariates related to C. |
n.boot |
Number of bootstraps to use in the estimation of bootstrap standard errors. |
n.iter |
Number of iterations; the default is |
eps |
Convergence error. This is set by the user in such away that the desired convergence is met; the default is |
Bootstrap standard errors for parameter estimates and for estimated cumulative hazard function.
This function evaluates the specified B-spline defined on the
unit interval, when considering n.if.per.cov B-splines. Currently, the
implementation is based on the one in Andrews, Shi 2013 (supplementary
materials).
Bspline.unit.interval(x, spline.index, n.if.per.cov, degree = 3)Bspline.unit.interval(x, spline.index, n.if.per.cov, degree = 3)
x |
value inside the unit interval at which to evaluate the spline. |
spline.index |
Index of the spline to evaluate. |
n.if.per.cov |
Number of B-splines to consider over the unit interval. |
degree |
Degree of the B-splines. Default is |
Andrews, D.W.K. and Shi, X. (2013). Inference based on confitional moment inequalities. Econometrica. 81(2):609-666.
This function calculates a bivariate survival probability based on multivariate normal distribution.
Bvprob(lx, ly, rho)Bvprob(lx, ly, rho)
lx |
The first lower bound of integration |
ly |
The second lower bound |
rho |
Association parameter |
This function combines a list of individual identified intervals to a single one based on majority vote. Note that the intersection of all intervals can be viewed as a majority vote as well, so that it is included as a special case.
cbMV(results.list, threshold)cbMV(results.list, threshold)
results.list |
List object containing the individual identified intervals. |
threshold |
Threshold proportion of identified intervals a given value
should be contained in in order for it to be included in the combined
identified interval. For intersection bounds, set this value to |
The combined identified interval.
This function checks whether the arguments supplied to the
main estimation function pi.surv are valid. When arguments are
invalid, the an exception is thrown.
check.args.pisurv( data, idx.param.of.interest, idxs.c, t, par.space, search.method, add.options )check.args.pisurv( data, idx.param.of.interest, idxs.c, t, par.space, search.method, add.options )
data |
Data frame containing the data on which to fit the model. The columns should be named as follows: 'Y' = observed timed, 'Delta' = censoring indicators, 'X0' = intercept column, 'X1' - 'Xp' = covariates. |
idx.param.of.interest |
Index of element in the covariate vector for
which the identified interval should be estimated. It can also be specified
as |
idxs.c |
Vector of indices of the continuous covariates. Suppose the
given data contains 5 covariates, of which 'X2' and 'X5' are continuous, this
argument should be specified as |
t |
Time point for which to estimate the identified set of
|
par.space |
Matrix containing bounds on the space of the parameters. The first column corresponds to lower bounds, the second to upper bounds. The i'th row corresponds to the bounds on the i'th element in the parameter vector. |
search.method |
The search method to be used to find the identified
interval. Default is |
add.options |
List of additional options to be specified to the method.
Notably, it can be used to select the link function |
This function transforms the parameters of the Cholesky de- composition to the covariance matrix, represented as a the row-wise con- catenation of the upper-triangular elements.
chol2par(par.chol1)chol2par(par.chol1)
par.chol1 |
The vector of Cholesky parameters. |
Covariance matrix corresponding to the provided Cholesky decomposition.
This function transforms the parameters of the Cholesky de- composition to a covariance matrix element. This function is used in chol2par.R.
chol2par.elem(a, b, par.chol1)chol2par.elem(a, b, par.chol1)
a |
The row index of the covariance matrix element to be computed. |
b |
The column index of the covariance matrix element to be computed. |
par.chol1 |
The vector of Cholesky parameters. |
Specified element of the covariance matrix resulting from the provided Cholesky decomposition.
R6 object that mimics a chronometer. It can be started, paused, record legs and stopped.
show()
Display the values stored in this chronometer object.
Chronometer$show()
reset()
Reset the chronometer.
Chronometer$reset()
start()
Start the chronometer
Chronometer$start()
stop()
Stop the chronometer
Chronometer$stop(leg.name = NULL)
leg.name(optional) Name for the stopped leg.
record.leg()
Record a leg time. The chronometer will continue running.
Chronometer$record.leg(leg.name = NULL)
leg.nameName for the recorded leg.
get.chronometer.data()
Like show method, but more rudimentary.
Chronometer$get.chronometer.data()
get.total.time()
Return the total time span between start and stop.
Chronometer$get.total.time(force = FALSE)
forceBoolean variable. If TRUE, avoids error when calling
this function while chronometer has not been stopped yet.
accumulate.legs()
Return total time spent per leg category (using leg names).
Chronometer$accumulate.legs(force = FALSE)
forceforce Boolean variable. If TRUE, avoids error when calling
this function while chronometer has not been stopped yet.
clone()
The objects of this class are cloneable with this method.
Chronometer$clone(deep = FALSE)
deepWhether to make a deep clone.
This function clears the plotting window
clear.plt.wdw()clear.plt.wdw()
This function estimates phi function at fixed time point t
CompC(theta, t, X, W, ld, cop, dist)CompC(theta, t, X, W, ld, cop, dist)
theta |
Estimated parameter values/initial values for finite dimensional parameters |
t |
A fixed time point |
X |
Data matrix with covariates related to T |
W |
Data matrix with covariates related to C. First column of W should be ones |
ld |
Output of |
cop |
Which copula should be computed to account for dependency between T and C. This argument can take
one of the values from |
dist |
The distribution to be used for the dependent censoring C. Only two distributions are allowed, i.e, Weibull
and lognormal distributions. With the value |
Prepare initial values within the control arguments
control.arguments(maxit = 300, eps = 1e-06, trace = TRUE, ktau.inits = NULL)control.arguments(maxit = 300, eps = 1e-06, trace = TRUE, ktau.inits = NULL)
maxit |
a positive integer that denotes the maximum iteration number in optimization. |
eps |
a positive small numeric value that denotes the tolerance for convergence. |
trace |
a logical value that judges whereh the tracing information on the progress of the model fitting should be produced. The default value if |
ktau.inits |
a numeric vector that contains initial values of the Kendall's tau. |
The distribution function of the Archimedean copula
copdist.Archimedean(x, copfam, ktau, coppar = NULL)copdist.Archimedean(x, copfam, ktau, coppar = NULL)
x |
the value at which the distribution function will be calculated at. |
copfam |
a character string that denotes the copula family. |
ktau |
a numeric value that denotes the Kendall's tau. |
coppar |
a numeric value that denotes the copula parameter. |
The h-function of the copula
cophfunc(x, coppar, copfam, condvar = 1)cophfunc(x, coppar, copfam, condvar = 1)
x |
the value at which the h-function will be calculated at. |
coppar |
a numeric value that denotes the copula parameter. |
copfam |
a character string that denotes the copula family. |
condvar |
a numeric value that specifies the type of the h-function. |
Convert the copula parameter the Kendall's tau
coppar.to.ktau(coppar, copfam)coppar.to.ktau(coppar, copfam)
coppar |
a numeric value that denotes the copula parameter. |
copfam |
a character string that denotes the copula family. |
This function implements the second step likelihood function of the competing risk model defined in Willems et al. (2024+).
cr.lik( n, s, Y, admin, cens.inds, M, Sigma, beta.mat, sigma.vct, rho.mat, theta.vct )cr.lik( n, s, Y, admin, cens.inds, M, Sigma, beta.mat, sigma.vct, rho.mat, theta.vct )
n |
The sample size. |
s |
The number of competing risks. |
Y |
The observed times. |
admin |
Boolean value indicating whether or not administrative censoring should be taken into account. |
cens.inds |
matrix of censoring indicators (each row corresponds to a single observation). |
M |
Design matrix, consisting of [intercept, exo.cov, Z, cf]. Note that
|
Sigma |
The covariance matrix. |
beta.mat |
Matrix containing all of the covariate effects. |
sigma.vct |
Vector of standard deviations. Should be equal to
|
rho.mat |
The correlation matrix. |
theta.vct |
Vector containing the parameters of the Yeo-Johnsontrans- formations. |
Evaluation of the log-likelihood function
Willems et al. (2024+). Flexible control function approach under competing risks (in preparation).
This function computes the diagonal matrix of the sample variance-covariance matrix.
D.hat(input, beta = NULL, t = NULL, hp = NULL, m.avg = NULL, mi.mat = NULL)D.hat(input, beta = NULL, t = NULL, hp = NULL, m.avg = NULL, mi.mat = NULL)
input |
Can either be the variance-covariance matrix obtained from the function Sigma.hat, or the data frame. |
beta |
The coefficient vector. Only needs to be supplied when the
argument for |
t |
The time point of interest. Only needs to be supplied when the
argument for |
hp |
List of hyperparameters. Only needs to be supplied when the
argument for |
m.avg |
See documentation of |
mi.mat |
See documentation of |
This function generates competing risk data that can be used in simulation studies.
dat.sim.reg.comp.risks( n, par, iseed, s, conf, Zbin, Wbin, type.cov, A.upper = 15 )dat.sim.reg.comp.risks( n, par, iseed, s, conf, Zbin, Wbin, type.cov, A.upper = 15 )
n |
The sample size of the generated data set. |
par |
List of parameter vectors for each component of the transformation model. |
iseed |
Random seed. |
s |
The number of competing risks. Note that the given parameter vector could specify the parameters for the model with more than s competing risks, but in this case only the first s sets of parameters will be considered. |
conf |
Boolean value indicating whether the data set should contain confounding. |
Zbin |
Indicator whether the confounded variable is binary
|
Wbin |
Indicator whether the instrument is binary ( |
type.cov |
Vector of characters "c" and "b", indicating which exogenous
covariates should be continuous |
A.upper |
The upper bound on the support of the administrative censoring
distribution. This can be used to control for the amount of administrative
censoring in the data. Default is |
A generated data set
This function defines the derivative of the transformation function that maps Cholesky parameters to the full covariance matrix.
dchol2par(par.chol1)dchol2par(par.chol1)
par.chol1 |
The vector of Cholesky parameters. |
Derivative of the function that transforms the cholesky parameters to the full covariance matrix.
This function defines the derivative of the transformation function that maps Cholesky parameters to a covariance matrix element. This function is used in dchol2par.R.
dchol2par.elem(k, q, a, b, par.chol1)dchol2par.elem(k, q, a, b, par.chol1)
k |
The row index of the parameter with respect to which to take the derivative. |
q |
the column index of the parameter with respect to which to take the derivative. |
a |
The row index of the covariance matrix element to be computed. |
b |
The column index of the covariance matrix element to be computed. |
par.chol1 |
The vector of Cholesky parameters. |
Derivative of the function that transforms the cholesky parameters to the specified element of the covariance matrix, evaluated at the specified arguments.
This function computes the matrix of sample derivatives of the sample variances.
dD.hat( data, beta, t, hp, mi.mat = NULL, m.avg = NULL, dm.avg = NULL, dmi.tens = NULL )dD.hat( data, beta, t, hp, mi.mat = NULL, m.avg = NULL, dm.avg = NULL, dmi.tens = NULL )
data |
Data frame. |
beta |
Vector of coefficients. |
t |
Time point at which to evaluate the (derivatives of) the moment functions. Also allowed to be a vector of time points (used in estimating the model under assumed time- independent coefficients). |
hp |
List of hyperparamerers. |
mi.mat |
A precomputed matrix of moment function evaluations at each
observation. If supplied, some computations can be skipped. Default is
|
m.avg |
A precomputed vector of the sample average of the moment
functions. If not supplied, this vector is computed. Default is
|
dm.avg |
Matrix of precomputed sample averages of the derivatives of the
moment functions. Default is |
dmi.tens |
3D tensor of precomputed evaluations of the derivatives of
the moment functions. Default is |
A matrix containing the partial derivatives of the variances of the moment functions. Each row corresponds to a moment function, each column corresponds to a covariate.
This function computes distance between two vectors based on L2-norm
This function computes distance between two vectors based on L2-norm
Distance(b, a) Distance(b, a)Distance(b, a) Distance(b, a)
b |
Second vector |
a |
First vector |
This function defines the derivative of the AFT link function.
dLambda_AFT_ll(t)dLambda_AFT_ll(t)
t |
time parameter. |
This function defines the derivative of the Cox PH link function.
dLambda_Cox_wb(t)dLambda_Cox_wb(t)
t |
time parameter. |
This function computes the matrix containing the sample average of the partial derivatives of the moment functions.
dm.bar(data, beta, t, hp, dmi.tens = NULL)dm.bar(data, beta, t, hp, dmi.tens = NULL)
data |
Data frame. |
beta |
Vector of coefficients. |
t |
Time point at which to compute the derivative of the moment functions. Also allowed to be a vector of time points (used in estimating the model under assumed time- independent coefficients). |
hp |
List of hyperparameters. |
dmi.tens |
Tensor of derivative moment function evaluations. Can be used
to avoid some computation. Default is |
A matrix containing the sample average of the partial derivatives of the moment functions. Each row corresponds to a moment function, each column corresponds to a coefficient.
This function finds the point for which the expected improvement is optimal, based on a given set of starting values. (M_step.R)
do.optimization.Mstep(start.vals, EI.Mstep, hyperparams)do.optimization.Mstep(start.vals, EI.Mstep, hyperparams)
start.vals |
Starting values for optimization. |
EI.Mstep |
Function to compute expected improvements. |
hyperparams |
List of hyperparameters. |
Maximum of the expected imrpovement function.
Used in the M-step (get.starting.values.R). ToDo: Adapt this code so as to also perform sample space contractions as in the MatLab implementation of Bei (2024).
draw.sv.init(theta.hash, dir, hyperparams, EI.fun)draw.sv.init(theta.hash, dir, hyperparams, EI.fun)
theta.hash |
Tentative optimal value for theta, i.e., the largest or smallest feasible value for theta (if dir = 1 or dir = -1, respectively). A 'feasible value' is one that satisfies all moment restrictions. |
dir |
Search direction. |
hyperparams |
List of hyperparameters. |
EI.fun |
Function used to compute the expected improvement. See also
|
Initial set of starting values.
Bei, X. (2024). Local linearieation based subvector inference in moment inequality models. Journal of Econometrics. 238:105549-
Evaluates the derivative of the Yeo-Johnson transformation at the provided argument.
DYJtrans(y, theta)DYJtrans(y, theta)
y |
The argument to be supplied to the derivative of the Yeo-Johnson transformation. |
theta |
The parameter of the Yeo-Johnson transformation. This should be a number in the range [0,2]. |
The transformed value of y.
This function performs the estimation step in the EAM algorithm.
E_step(thetas, test.fun, dir, evaluations, verbose)E_step(thetas, test.fun, dir, evaluations, verbose)
thetas |
Points at which to perform the E-step. Usually the result of the M-step. |
test.fun |
Function returning the test statistic, as well as the critical value. |
dir |
Direction in which to optimize. For finding upper bounds, set
|
evaluations |
Matrix containing each point that was already evaluated, alongside the corresponding test statistic and critical value, as its rows. |
verbose |
Verbosity parameter. |
Results of the E-step.
This function implements the EAM search strategy as described in
Kaido, Molinari and Stoye (2019). This is a generic function, requiring the
specification of a test function (test.fun), as well as the
specification of the parameter space (hyperparams).
EAM( dir, test.fun, hyperparams, evaluations = NULL, time.run.duration = FALSE, verbose = 0, picturose = FALSE )EAM( dir, test.fun, hyperparams, evaluations = NULL, time.run.duration = FALSE, verbose = 0, picturose = FALSE )
dir |
The direction of the test. |
test.fun |
The test function to be inverted in order to obtain the
identified set. It should take a scalar parameter as argument (i.e. the
specified value of a component of the full parameter vector) and return a
list with named elements |
hyperparams |
A list of hyperparameters needed in order for this method
to run (see |
evaluations |
Matrix of already evaluated points, of which at least one
is feasible. When |
time.run.duration |
Boolean value indicating whether to time each step
in the EAM algorithm. Requires |
verbose |
Boolean value indicating whether or not to print run time
updates to the console. Default is |
picturose |
Boolean value indicating whether or not to visualize the
identified set search. Default is |
List containing the evaluations of the test statistic and critical values, convergence information, and run times.
Kaido et al. (2019). Confidence intervals for projections of partially identified parameters. Econometrica. 87(4):1397-1432.
This function checks the convergence of the EAM algorithm. ToDo: Get rid of hard coding stop of algorithm when no more improvement of theta values (maybe related to parameter space contraction, since the problem is that the given points to check in the E-step of the following iteration can always be the same and always be rejected (except of course for the randomly chosen one), while the most promising theta value continues to be the same, infeasible value. In this way, it is possible that theta.hash - mp.theta.next at some point will never decrease).
EAM.converged( opt.val.prev, evaluations, mp.theta.next, iter.nbr, dir, hyperparams, verbose )EAM.converged( opt.val.prev, evaluations, mp.theta.next, iter.nbr, dir, hyperparams, verbose )
opt.val.prev |
Previous optimal theta value. |
evaluations |
Matrix of violation curve evaluations. |
mp.theta.next |
Most promising value of theta for which to run the E-step in the following iteration |
iter.nbr |
Number of iterations of the EAM algorithm run so far. |
dir |
Search direction. |
hyperparams |
List of hyperparameters used in the EAM algorithm. |
verbose |
Verbosity parameter. |
Boolean value whether or not algorithm has converged.
Used in the M-step (M_step.R). Note: predict(fit.krige, ...) has weird beheviour when making predictions for a single value in terms of standard error. We work around this issue in this implementation.
EI(theta, test.fun, fit.krige, theta.hash, dir)EI(theta, test.fun, fit.krige, theta.hash, dir)
theta |
Vector of coefficients. |
test.fun |
Test function (cf. |
fit.krige |
Fitted Kriging model. |
theta.hash |
Tentative optimal value for theta, i.e., the largest or smallest feasible value for theta (if dir = 1 or dir = -1, respectively). A 'feasible value' is one that satisfies all moment restrictions. |
dir |
Search direction. |
The expected improvement.
This function estimates the control function for the endogenous variable based on the provided covariates. This function is called inside estimate.cmprsk.R.
estimate.cf(XandW, Z, Zbin, gammaest = NULL)estimate.cf(XandW, Z, Zbin, gammaest = NULL)
XandW |
Design matrix of exogenous covariates. |
Z |
Endogenous covariate. |
Zbin |
Boolean value indicating whether endogenous covariate is binary. |
gammaest |
Vector of pre-estimated parameter vector. If |
List containing the vector of values for the control function and the regression parameters of the first step.
This function estimates the parameters in the competing risks model described in Willems et al. (2024+). Note that this model extends the model of Crommen, Beyhum and Van Keilegom (2024) and as such, this function also implements their methodology.
estimate.cmprsk( data, admin, conf, eoi.indicator.names = NULL, Zbin = NULL, inst = "cf", realV = NULL, compute.var = TRUE, eps = 0.001 )estimate.cmprsk( data, admin, conf, eoi.indicator.names = NULL, Zbin = NULL, inst = "cf", realV = NULL, compute.var = TRUE, eps = 0.001 )
data |
A data frame, adhering to the following formatting rules:
|
admin |
Boolean value indicating whether the data contains administrative censoring. |
conf |
Boolean value indicating whether the data contains confounding
and hence indicating the presence of |
eoi.indicator.names |
Vector of names of the censoring indicator columns
pertaining to events of interest. Events of interest will be modeled allowing
dependence between them, whereas all censoring events (corresponding to
indicator columns not listed in |
Zbin |
Indicator value indicating whether ( |
inst |
Variable encoding which approach should be used for dealing with
the confounding. |
realV |
Vector of numerics with length equal to the number of rows in
|
compute.var |
Boolean value indicating whether the variance of the
parameter estimates should be computed as well (this can be very
computationally intensive, so may want to be disabled). Default is
|
eps |
Value that will be added to the diagonal of the covariance matrix during estimation in order to ensure strictly positive variances. |
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.
Willems et al. (2024+). Flexible control function approach under competing risks (in preparation).
Crommen, G., Beyhum, J., and Van Keilegom, I. (2024). An instrumental variable approach under dependent censoring. Test, 33(2), 473-495.
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)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)
Also called the 'initialization' step in KMS19, this method tries to find at least one initial feasible point, which is required to run the EAM algorithm. ToDo: Investigate whether the feasible point search of Bei (2024) is better. If so, implement it.
feasible_point_search( test.fun, hyperparams, verbose, picturose = FALSE, parallel = FALSE )feasible_point_search( test.fun, hyperparams, verbose, picturose = FALSE, parallel = FALSE )
test.fun |
Function that takes a parameter vector as a first argument and returns the test statistic, as well as the critical value. |
hyperparams |
List of hyperparameters. |
verbose |
Verbosity parameter. |
picturose |
Picturosity flag. If |
parallel |
Flag for whether or not parallel computing should be used.
Default is |
Results of the initial feasible point search.
Kaido et al. (2019). Confidence intervals for projections of partially identified parameters. Econometrica. 87(4):1397-1432.
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).
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 )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 )
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T. |
W |
Data matrix with covariates related to C. First column of W should be a vector of ones. |
cop |
Which copula should be computed to account for dependency between T and C. This argument can take
one of the values from |
dist |
The distribution to be used for the censoring time C. Only two distributions are allowed, i.e, Weibull
and lognormal distributions. With the value |
start |
Initial values for the finite dimensional parameters. If |
n.iter |
Number of iterations; the default is |
bootstrap |
A boolean indicating whether to compute bootstrap standard errors for making inferences. |
n.boot |
Number of bootstrap samples to use in the estimation of bootstrap standard errors if |
ncore |
The number of cores to use for parallel computation in bootstrapping, with the default |
eps |
Convergence error. This is set by the user in such away that the desired convergence is met; the default is |
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.
Deresa and Van Keilegom (2024). Copula based Cox proportional hazards models for dependent censoring, Journal of the American Statistical Association, 119:1044-1054.
# 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")# 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")
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.
fitIndepCens( resData, X, W, dist = c("Weibull", "lognormal"), start = NULL, n.iter = 50, bootstrap = TRUE, n.boot = 150, ncore = 7, eps = 1e-04 )fitIndepCens( resData, X, W, dist = c("Weibull", "lognormal"), start = NULL, n.iter = 50, bootstrap = TRUE, n.boot = 150, ncore = 7, eps = 1e-04 )
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T. |
W |
Data matrix with covariates related to C. First column of W should be ones. |
dist |
The distribution to be used for the censoring time C. Only two distributions are allowed, i.e, Weibull
and lognormal distributions. With the value |
start |
Initial values for the finite dimensional parameters. If |
n.iter |
Number of iterations; the default is |
bootstrap |
A boolean indicating whether to compute bootstrap standard errors for making inferences. |
n.boot |
Number of bootstrap samples to use in the estimation of bootstrap standard errors if |
ncore |
The number of cores to use for parallel computation is configurable, with the default |
eps |
Convergence error. This is set by the user in such away that the desired convergence is met; the default is |
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.
# 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")# 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")
This function defined the class of box functions as defined in Willems et al. (2024+).
G.box( x, g.idx, data, n.box.per.cov, norm.func, cov.ranges = NULL, norm.cov.out = NULL, ... )G.box( x, g.idx, data, n.box.per.cov, norm.func, cov.ranges = NULL, norm.cov.out = NULL, ... )
x |
Vector of covariates to be normalized alongside the data. Default is
|
g.idx |
Index of the instrumental function, in {1, ..., n.inst.func}. |
data |
Data frame. |
n.box.per.cov |
Number of box functions to consider per continuous covariate. |
norm.func |
Function to be used to normalize the covariates. |
cov.ranges |
Matrix of ranges of the covariates. Used for normalizing
the data to the unit interval before applying the instrumental functions.
Default is |
norm.cov.out |
Output of a preliminary call to the supplied covariate normalization function. |
... |
Additional arguments will be ignored. Useful for allowing
compatibility with the implementations of other instrument function families.
Specifically, it allows to ignore the |
The function normalizes the continuous covariates to lie in the unit interval and then evaluates the subvector of continuous covariates on the specified family of instrumental function. For the discrete elements, indicator functions are used for each level.
G.cd( x, g.idx, data, n.if.per.cov, idxs.c, G.c, norm.func, discrete.covariate.levels = NULL, cov.ranges = NULL, norm.cov.out = NULL, degree = 3 )G.cd( x, g.idx, data, n.if.per.cov, idxs.c, G.c, norm.func, discrete.covariate.levels = NULL, cov.ranges = NULL, norm.cov.out = NULL, degree = 3 )
x |
The vector of covariates at which to evaluate the B-splines |
g.idx |
The index of the instrumental function. |
data |
Data frame containing the data. |
n.if.per.cov |
Number of instrumental functions per continuous covariate. |
idxs.c |
Vector of indices of the continuous elements in the vector of covariates. |
G.c |
Family of instrumental functions to use for the subvector of continuous covariates. |
norm.func |
Function to be used to normalize the covariates. |
discrete.covariate.levels |
Matrix containing as rows all possible
'combined' levels of the discrete covariates. Default is
|
cov.ranges |
Matrix containing as its rows the lower and upper bounds
for each continuous covariate. Default is |
norm.cov.out |
Output of a preliminary call to a covariate normalization
function (defined above). This is used to speed up computations. Note that
this argument should only apply to continuous covariates!! Default is
|
degree |
Degree of the spline functions to be used as instrumental
functions for the continuous covariates (if applicable). Default is
|
This function defines the family of discrete/continuous instrumental functions in the case of many covariates. It does so by considering a instrumental functions for each pair of entries in the given covariate vector.
G.cd.mc( x, g.idx, data, n.if.per.cov, idxs.c, G.c, norm.func, info.manycov = NULL, cov.ranges = NULL, norm.cov.out = NULL, degree = 3, ... )G.cd.mc( x, g.idx, data, n.if.per.cov, idxs.c, G.c, norm.func, info.manycov = NULL, cov.ranges = NULL, norm.cov.out = NULL, degree = 3, ... )
x |
The vector of covariates at which to evaluate the B-splines |
g.idx |
The index of the instrumental function. |
data |
Data frame containing the data. |
n.if.per.cov |
Number of instrumental functions per continuous covariate. |
idxs.c |
Vector of indices of the continuous elements in the vector of covariates. |
G.c |
Family of instrumental functions to use for the subvector of continuous covariates. |
norm.func |
Function to be used to normalize the covariates. |
info.manycov |
Data frame containing some information about the global
structure of the instrumental functions of this class. If
|
cov.ranges |
Matrix containing as its rows the lower and upper bounds
for each continuous covariate. Default is |
norm.cov.out |
Output of function that normalizes the covariates. |
degree |
Degree of the spline functions to be used as instrumental
functions for the continuous covariates (if applicable). Default is
|
... |
Arguments specified here will be ignored. Used for compatibility with other instrumental function classes. |
Compute the Gn matrix in step 3b of Bei (2024).
G.hat( data, beta, t, hp, mi.mat = NULL, m.avg = NULL, dm.avg = NULL, dmi.tens = NULL, D = NULL )G.hat( data, beta, t, hp, mi.mat = NULL, m.avg = NULL, dm.avg = NULL, dmi.tens = NULL, D = NULL )
data |
Data frame. |
beta |
Vector of coefficients. |
t |
Time point at which to evaluate the (derivatives of) the moment functions. |
hp |
List of hyperparamerers. |
mi.mat |
A precomputed matrix of moment function evaluations at each
observation. If supplied, some computations can be skipped. Default is
|
m.avg |
A precomputed vector of the sample average of the moment
functions. If not supplied, this vector is computed. Default is
|
dm.avg |
Matrix of precomputed sample averages of the derivatives of the
moment functions. Default is |
dmi.tens |
3D tensor of precomputed evaluations of the derivatives of
the moment functions. Default is |
D |
Diagonal of D-matrix. |
A matrix containing the partial derivatives of the variances of the moment functions. Each row corresponds to a moment function, each column corresponds to a covariate.
Bei, X. (2024). Local linearieation based subvector inference in moment inequality models. Journal of Econometrics. 238:105549-
This function normalizes the covariates to lie in the unit interval and then evaluates each B-spline at each observation, multiplying together the results per observation.
G.spline( x, g.idx, data, n.if.per.cov, norm.func, cov.ranges = NULL, norm.cov.out = NULL, degree = 3 )G.spline( x, g.idx, data, n.if.per.cov, norm.func, cov.ranges = NULL, norm.cov.out = NULL, degree = 3 )
x |
The vector of covariates at which to evaluate the B-splines |
g.idx |
The index of the instrumental function. Note that g.idx ranges between 1 and n.if.per.cov^n.cov, as an instrumental function is the product of the appropriate B-spline evaluation for each element in the covariate vector. |
data |
Data frame containing the data. |
n.if.per.cov |
Number of instrumental variables to be used per covariate. |
norm.func |
Function to be used to normalize the covariates. |
cov.ranges |
Matrix of ranges of the covariates. Used for normalizing
the covariates. If |
norm.cov.out |
Output of a preliminary call to the given covariate
normalization function. Default is |
degree |
Degree of B-splines to use. Default value is |
The generator function of the Archimedean copula
generator.Archimedean(x, coppar, copfam, inverse = FALSE)generator.Archimedean(x, coppar, copfam, inverse = FALSE)
x |
the value at which the generator function will be calculated at. |
coppar |
a numeric value that denotes the copula parameter. |
copfam |
a character string that denotes the copula family. |
inverse |
a logical value that specifies whether the inverse function will be used. |
The points returned by this function can be used as corner points in the family of box functions, or as knots in the family of B-spline functions.
get.anchor.points(data, n.if.per.cov, normalized = FALSE)get.anchor.points(data, n.if.per.cov, normalized = FALSE)
data |
Data set. |
n.if.per.cov |
Number of instrumental functions to use per continuous covariate. |
normalized |
Boolean value indicating whether the covariates in the
given data frame have been normalized. Default is |
This function computes the 1(Y <= t) - Lambda(X^T beta(t)) and Lambda(X^T beta(t)) - 1(Y <= t, Delta = 1) parts of the moment functions. (Used in function get.mi.mat.R)
get.cond.moment.evals(data, beta, t, hp)get.cond.moment.evals(data, beta, t, hp)
data |
Data frame. |
beta |
Vector of coefficients. |
t |
Time point of interest. |
hp |
List of hyperparameters. |
A vector of 2n elements containing in the first n positions the evaluations of 1(Y <= t) - Lambda(X^T beta(t)) and in the last n positions the evaluations of Lambda(X^T beta(t)) - 1(Y <= t, Delta = 1).
This function computes the critical value following the algorithm of Section 4.3 in Bei (2024).
get.cvLLn( BetaI.r, data, t, hp, c, r, par.space, inst.func.evals = NULL, alpha = 0.95 )get.cvLLn( BetaI.r, data, t, hp, c, r, par.space, inst.func.evals = NULL, alpha = 0.95 )
BetaI.r |
Matrix containing in its columns the minimizers of the S-function leading to the test statistic. |
data |
Data frame. |
t |
Time point of interest. Also allowed to be a vector of time points (used in estimating the model under assumed time- independent coefficients). |
hp |
List of hyperparameters. |
c |
Projection vector. |
r |
Result of projection of parameter vector onto |
par.space |
Bounds on the parameter space. |
inst.func.evals |
Matrix of precomputed instrumental function
evaluations for each observation in the data set. If |
alpha |
Confidence level. |
The critical value for the test statistic.
Bei, X. (2024). Local linearieation based subvector inference in moment inequality models. Journal of Econometrics. 238:105549-
This function evaluates the derivatives of the conditional moment function at each observation. Used in get.dmi.tens.R
get.deriv.mom.func(data, beta, t, hp)get.deriv.mom.func(data, beta, t, hp)
data |
Data frame. |
beta |
Parameter vector. |
t |
Time point of interest. |
hp |
List of hyperparameters. |
This function provides a faster implementation of obtaining the evaluations of the derivative of the moment functions at each observation (wrt the previous implementation using 'dm.comp' and 'dm.R'). Used in the function G.hat.R
get.dmi.tens(data, beta, t, hp, inst.func.evals = NULL)get.dmi.tens(data, beta, t, hp, inst.func.evals = NULL)
data |
Data frame. |
beta |
Vector of coefficients. |
t |
Time point of interest. Also allowed to be a vector of time points (used in estimating the model under assumed time- independent coefficients). |
hp |
List of hyperparameters. |
inst.func.evals |
Precomputed matrix of instrumental function
evaluations. Defaults is |
Function used to obtain extra theta values to be supplied to the
E-step in the next iteration (M_step.R). Note: this function should be
changed when implementing the sample space contractions (see comment made in
documentation of M_step).
get.extra.Estep.points(dir, theta.hash, maxviol.hash, hyperparams)get.extra.Estep.points(dir, theta.hash, maxviol.hash, hyperparams)
dir |
Search direction. |
theta.hash |
Tentative optimal value for theta, i.e., the largest or smallest feasible value for theta (if dir = 1 or dir = -1, respectively). A 'feasible value' is one that satisfies all moment restrictions. |
maxviol.hash |
Violation curve evaluated at |
hyperparams |
List of hyperparameters. |
Points to evaluate in E-step.
Obtain the evaluations of each observation on each of the instrumental functions. (Used in function get.mi.mat.R)
get.instrumental.function.evals(data, hp)get.instrumental.function.evals(data, hp)
data |
Data frame. |
hp |
List of hyperparameters. Notably, it contains the instrumental
function to be used in an element named |
This function obtains the moment function evaluations.
get.mi.mat(data, beta, t, hp, inst.func.evals = NULL)get.mi.mat(data, beta, t, hp, inst.func.evals = NULL)
data |
Data frame. |
beta |
Vector of coefficients. |
t |
Time point at which to compute the moment function. Also allowed to be a vector of time points (used in estimating the model under assumed time- independent coefficients). |
hp |
List of hyperparameters |
inst.func.evals |
Matrix of instrumental function evaluations. If
|
Function to obtain the next point to evaluate in the search for a feasible point. This function evaluates the entire parameter space of the component of theta as evenly as possible. Used in the initialization step (feasible_point_search.R)
get.next.point(evaluations, lb.theta, ub.theta)get.next.point(evaluations, lb.theta, ub.theta)
evaluations |
Matrix of evaluations performed so far. |
lb.theta |
Lower bound on the parameter of interest. |
ub.theta |
Upper bound on the parameter of interest. |
Next point in the feasible point search.
Obtain starting values used in the M-step (M_step.R).
get.starting.values(theta.hash, dir, EI.Mstep, hyperparams)get.starting.values(theta.hash, dir, EI.Mstep, hyperparams)
theta.hash |
Tentative optimal value for theta, i.e., the largest or smallest feasible value for theta (if dir = 1 or dir = -1, respectively). A 'feasible value' is one that satisfies all moment restrictions. |
dir |
Search direction. |
EI.Mstep |
Function to compute expected improvements. |
hyperparams |
List of hyperparameters. |
Vector of starting values
Obtain the test statistic by minimizing the S-function over the
feasible region .
get.test.statistic( beta.init, data, par.space, t, hp, c, r, inst.func.evals = NULL )get.test.statistic( beta.init, data, par.space, t, hp, c, r, inst.func.evals = NULL )
beta.init |
Starting value of minimization algorithm. |
data |
Data frame. |
par.space |
Matrix containing the bounds on the parameter space. |
t |
Time point at which to evaluate beta. Also allowed to be a vector of time points (used in estimating the model under assumed time- independent coefficients). |
hp |
List of hyperparameters. |
c |
Projection vector |
r |
hypothesised value of the projection. |
inst.func.evals |
Matrix of precomputed instrumental function
evaluations for each observation in the data set. If |
A list containing the value of the test statistic and the parameter at which this value was attained.
This function implements the gridsearch and binary search algorithms used to compute the roots of the violation curve and hence in estimating the identified intervals.
gridSearch( dir, test.fun, hyperparams, evaluations = NULL, time.run.duration = FALSE, verbose = 0, picturose = FALSE )gridSearch( dir, test.fun, hyperparams, evaluations = NULL, time.run.duration = FALSE, verbose = 0, picturose = FALSE )
dir |
Search direction. |
test.fun |
The test function to be inverted in order to obtain the
identified set. It should take a scalar parameter as argument (i.e. the
specified value of a component of the full parameter vector) and return a
list with named elements |
hyperparams |
List of hyperparameters. |
evaluations |
Matrix of already evaluated points, of which at least one
is feasible. When |
time.run.duration |
Boolean value indicating whether to time each step
in the EAM algorithm. Requires |
verbose |
Boolean value indicating whether or not to print run time
updates to the console. Default is |
picturose |
Boolean value indicating whether or not to visualize the
identified set search. Default is |
List containing the evaluations of the test statistic and critical values, convergence information, and run times.
This function implements a rudimentary, bidirectional search algorithm. It works by expanding a grid with given step.size in both directions, starting from an initial feasible point.
gs.algo.bidir(test.results, max.iter, step.size)gs.algo.bidir(test.results, max.iter, step.size)
test.results |
Matrix containing the evaluations of the test statistic and critical value. |
max.iter |
Maximum number of iterations. |
step.size |
Step size based on which the grid is constructed. |
The next point to evaluate in the grid search.
This function implements the binary search algorithm, that starts from a given feasible point and looks in the given direction for the root of the violation curve.
gs.binary(evaluations, dir, iter.nbr, hp)gs.binary(evaluations, dir, iter.nbr, hp)
evaluations |
Matrix of evaluated test statistics and critical values. |
dir |
Search direction. |
iter.nbr |
Iteration number. |
hp |
List of hyperparameters. |
The next point to evaluate.
This function implements the interpolation search algorithm, that starts from a given feasible point and looks in the given direction for the root of the violation curve.
gs.interpolation(evaluations, dir, iter.nbr, hp)gs.interpolation(evaluations, dir, iter.nbr, hp)
evaluations |
Matrix of evaluated test statistics and critical values. |
dir |
Search direction. |
iter.nbr |
Iteration number. |
hp |
List of hyperparameters. |
The next point to evaluate.
This function implements a unidirectional grid search, that works by expanding a grid starting from a given feasible point in the given direction.
gs.regular(evaluations, dir, iter.nbr, hp)gs.regular(evaluations, dir, iter.nbr, hp)
evaluations |
Matrix of evaluated test statistics and critical values. |
dir |
Search direction. |
iter.nbr |
Iteration number. |
hp |
List of hyperparameters. |
Next point to evaluate in the search algorithm.
Used in initalization step (feasible_point_search.R).
insert.row(evaluations, row, idx.after)insert.row(evaluations, row, idx.after)
evaluations |
Matrix of violation function evaluations. |
row |
Row (evaluations) to be added to the evaluation matrix. |
idx.after |
Index of the row of |
Evaluation matrix.
Computes the inverse Yeo-Johnson transformation of the provided argument.
IYJtrans(y, theta)IYJtrans(y, theta)
y |
The argument to be supplied to the inverse Yeo-Johnson transformation. |
theta |
The parameter of the inverted Yeo-Johnson transformation. This should be a number in the range [0,2]. |
The transformed value of y.
Calculate the kernel function
Kernel(u, name = "Gaussian")Kernel(u, name = "Gaussian")
u |
the value in which the kernel function will be calculated at. |
name |
a character used to specify the type of kernel function. |
Convert the Kendall's tau into the copula parameter
ktau.to.coppar(ktau, copfam)ktau.to.coppar(ktau, copfam)
ktau |
a numeric value that denotes the Kendall's tau. |
copfam |
a character string that denotes the copula family. |
This function defines the AFT link function.
Lambda_AFT_ll(t)Lambda_AFT_ll(t)
t |
time parameter. |
This function defines the Cox PH link function.
Lambda_Cox_wb(t)Lambda_Cox_wb(t)
t |
time parameter. |
This function defines the inverse of the AFT link function.
Lambda_inverse_AFT_ll(p)Lambda_inverse_AFT_ll(p)
p |
probability. |
This function defines the inverse of the Cox PH link function.
Lambda_inverse_Cox_wb(p)Lambda_inverse_Cox_wb(p)
p |
probability. |
This function defines the loss function used in computing the penalized local linear approximation of the test statistic in order to construct the bootstrap distribution of the test statistic.
lf.delta.beta1( Delta.sub, vnb, phi, Gn, Omegan, beta, c, r, data, par.space, epsilon.n, lambda.n )lf.delta.beta1( Delta.sub, vnb, phi, Gn, Omegan, beta, c, r, data, par.space, epsilon.n, lambda.n )
Delta.sub |
Subvector of Delta. |
vnb |
Bootstrapped stochastic process. |
phi |
Moment selection functions. |
Gn |
First-order approximation matrix. |
Omegan |
Correlation matrix of sample moment functions. |
beta |
Coefficient vector. |
c |
Projection vector. |
r |
Value of projected coefficient vector. |
data |
Data frame. |
par.space |
Matrix containing the bounds on the parameter space. |
epsilon.n |
Parameter used in constructing the feasible region as in Example 4.1 in Bei (2024). Not used in this function. |
lambda.n |
Weight of penalty term. |
Loss function evaluation evaluated at the given Delta.
Bei, X. (2024). Local linearieation based subvector inference in moment inequality models. Journal of Econometrics. 238:105549-
This function implements the loss function used in computing the test statistic.
lf.ts(beta.sub, data, t, hp, c, r, inst.func.evals = NULL)lf.ts(beta.sub, data, t, hp, c, r, inst.func.evals = NULL)
beta.sub |
Subvector of coefficient vector. |
data |
Data frame. |
t |
Time point of interest. Also allowed to be a vector of time points (used in estimating the model under assumed time- independent coefficients). |
hp |
List of hyperparameters. |
c |
Unit vector containing unity at the location of the parameter of interest. |
r |
Value of the parameter of interest that is tested. |
inst.func.evals |
Pre-computed matrix of insturmental function evaluations. If not supplied, it will be computed during execution of this function. |
S-functions evaluation for the specified parameter vector.
Loglikehood function under independent censoring
LikCopInd(theta, resData, X, W, lhat, cumL, dist)LikCopInd(theta, resData, X, W, lhat, cumL, dist)
theta |
Estimated parameter values/initial values for finite dimensional parameters |
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T |
W |
Data matrix with covariates related to C. First column of W should be ones |
lhat |
The estimated hazard function obtained from the output of |
cumL |
The estimated cumulative hazard function from the output of |
dist |
The distribution to be used for the dependent censoring C. Only two distributions are allowed, i.e, Weibull
and lognormal distributions. With the value |
Maximized log-likelihood value
Calculate the likelihood function for the fully parametric joint distribution
Likelihood.Parametric(param, yobs, delta, copfam, margins, cure = FALSE)Likelihood.Parametric(param, yobs, delta, copfam, margins, cure = FALSE)
param |
a vector contains all parametric parameters. |
yobs |
a numeric vector that indicated the observed survival times. |
delta |
a numeric vector that stores the right-censoring indicators. |
copfam |
a character string that specifies the copula family. |
margins |
a list used to define the distribution structures of both the survival and censoring margins. |
cure |
a logical value that indicates whether the existence of a cured fraction should be considered. |
Calculate the profiled likelihood function with kernel smoothing
Likelihood.Profile.Kernel(param, yobs, delta, copfam, margins, cure = FALSE)Likelihood.Profile.Kernel(param, yobs, delta, copfam, margins, cure = FALSE)
param |
a vector contains all parametric parameters. |
yobs |
a numeric vector that indicated the observed survival times. |
delta |
a numeric vector that stores the right-censoring indicators. |
copfam |
a character string that specifies the copula family. |
margins |
a list used to define the distribution structures of both the survival and censoring margins. |
cure |
a logical value that indicates whether the existence of a cured fraction should be considered. |
Solve the profiled likelihood function
Likelihood.Profile.Solve( yobs, delta, copfam, margins, ktau.init, parapar.init, cure, curerate.init, constraints, maxit, eps )Likelihood.Profile.Solve( yobs, delta, copfam, margins, ktau.init, parapar.init, cure, curerate.init, constraints, maxit, eps )
yobs |
a numeric vector that indicated the observed survival times. |
delta |
a numeric vector that stores the right-censoring indicators. |
copfam |
a character string that specifies the copula family. |
margins |
a list used to define the distribution structures of both the survival and censoring margins. |
ktau.init |
initial value of Kendall's tau. |
parapar.init |
initial value of parametric parameters. |
cure |
a logical value that indicates whether the existence of a cured fraction should be considered. |
curerate.init |
initial value of cure rate. |
constraints |
constraints of parameters. |
maxit |
a positive integer that denotes the maximum iteration number in optimization. |
eps |
a positive small numeric value that denotes the tolerance for convergence. |
Calculate the semiparametric version of profiled likelihood function
Likelihood.Semiparametric( param, Syobs, yobs, delta, copfam, margins, cure = FALSE )Likelihood.Semiparametric( param, Syobs, yobs, delta, copfam, margins, cure = FALSE )
param |
a vector contains all parametric parameters. |
Syobs |
values of survival function at observed time points. |
yobs |
a numeric vector that indicated the observed survival times. |
delta |
a numeric vector that stores the right-censoring indicators. |
copfam |
a character string that specifies the copula family. |
margins |
a list used to define the distribution structures of both the survival and censoring margins. |
cure |
a logical value that indicates whether the existence of a cured fraction should be considered. |
This function defines the log-likelihood used to estimate the second step in the competing risks extension of the model described in Willems et al. (2024+).
LikF.cmprsk(par, data, admin, conf, cf)LikF.cmprsk(par, data, admin, conf, cf)
par |
Vector of all second step model parameters, consisting of the regression parameters, variance-covariance matrix elements and transformation parameters. |
data |
Data frame resulting from the 'uniformize.data.R' function. |
admin |
Boolean value indicating whether the data contains administrative censoring. |
conf |
Boolean value indicating whether the data contains confounding and hence indicating the presence of Z and W. |
cf |
"Control function" to be used. This can either be the (i) estimated
control function, (ii) the true control function, (iii) the instrumental
variable, or (iv) nothing ( |
Log-likelihood evaluation of the second step.
Willems et al. (2024+). Flexible control function approach under competing risks (in preparation).
This function parametrizes the covariance matrix using its Cholesky decomposition, so that optimization of the likelihood can be done based on this parametrization, and positive-definiteness of the covariance matrix is guaranteed at each step of the optimization algorithm.
likF.cmprsk.Cholesky(par.chol, data, admin, conf, cf, eps = 0.001)likF.cmprsk.Cholesky(par.chol, data, admin, conf, cf, eps = 0.001)
par.chol |
Vector of all second step model parameters, consisting of the regression parameters, Cholesky decomposition of the variance-covariance matrix elements and transformation parameters. |
data |
Data frame resulting from the 'uniformize.data.R' function. |
admin |
Boolean value indicating whether the data contains administrative censoring. |
conf |
Boolean value indicating whether the data contains confounding and hence indicating the presence of Z and W. |
cf |
"Control function" to be used. This can either be the (i) estimated
control function, (ii) the true control function, (iii) the instrumental
variable, or (iv) nothing ( |
eps |
Minimum value for the diagonal elements in the covariance matrix.
Default is |
Log-likelihood evaluation of the second step.
This function defines the maximum likelihood used to estimate the control function in the case of a continuous endogenous variable.
LikGamma1(gamma, Z, M)LikGamma1(gamma, Z, M)
gamma |
Vector of coefficients in the linear model used to estimate the control function. |
Z |
Endogenous covariate. |
M |
Design matrix, containing an intercept, the exogenous covariates and the instrumental variable. |
Returns the log-likelihood function corresponding to the data,
evaluated at the point gamma.
This function defines the maximum likelihood used to estimate the control function in the case of a binary endogenous variable.
LikGamma2(gamma, Z, M)LikGamma2(gamma, Z, M)
gamma |
Vector of coefficients in the logistic model used to estimate the control function. |
Z |
Endogenous covariate. |
M |
Design matrix, containing an intercept, the exogenous covariates and the instrumental variable. |
Returns the log-likelihood function corresponding to the data,
evaluated at the point gamma.
This function defines the log-likelihood used in estimating the second step in the competing risks extension of the model described in Willems et al. (2024+). The results of this function will serve as starting values for subsequent optimizations (LikI.comprsk.R and LikF.cmprsk.R)
LikI.bis(par, data, admin, conf, cf)LikI.bis(par, data, admin, conf, cf)
par |
Vector of all second step model parameters, consisting of the regression parameters, variance-covariance matrix elements and transformation parameters. |
data |
Data frame resulting from the 'uniformize.data.R' function. |
admin |
Boolean value indicating whether the data contains administrative censoring. |
conf |
Boolean value indicating whether the data contains confounding and hence indicating the presence of Z and W |
cf |
"Control function" to be used. This can either be the (i) estimated
control function, (ii) the true control function, (iii) the instrumental
variable, or (iv) nothing ( |
Starting values for subsequent optimization function used in the second step of the estimation procedure.
Willems et al. (2024+). Flexible control function approach under competing risks (in preparation).
This function defines the log-likelihood used to estimate the second step in the competing risks extension assuming independence of some of the competing risks in the model described in Willems et al. (2024+).
LikI.cmprsk(par, data, eoi.indicator.names, admin, conf, cf)LikI.cmprsk(par, data, eoi.indicator.names, admin, conf, cf)
par |
Vector of all second step model parameters, consisting of the regression parameters, variance-covariance matrix elements and transformation parameters. |
data |
Data frame resulting from the 'uniformize.data.R' function. |
eoi.indicator.names |
Vector of names of the censoring indicator columns
pertaining to events of interest. Events of interest will be modeled allowing
dependence between them, whereas all censoring events (corresponding to
indicator columns not listed in |
admin |
Boolean value indicating whether the data contains administrative censoring. |
conf |
Boolean value indicating whether the data contains confounding and hence indicating the presence of Z and W |
cf |
"Control function" to be used. This can either be the (i) estimated
control function, (ii) the true control function, (iii) the instrumental
variable, or (iv) nothing ( |
Log-likelihood evaluation for the second step in the esimation procedure.
Willems et al. (2024+). Flexible control function approach under competing risks (in preparation).
This function does the same as LikI.cmprsk (in fact, it even calls said function), but it parametrizes the covariance matrix using its Cholesky decomposition in order to guarantee positive definiteness. This function is never used, might not work and could be deleted.
LikI.cmprsk.Cholesky( par.chol, data, eoi.indicator.names, admin, conf, cf, eps = 0.001 )LikI.cmprsk.Cholesky( par.chol, data, eoi.indicator.names, admin, conf, cf, eps = 0.001 )
par.chol |
Vector of all second step model parameters, consisting of the regression parameters, Cholesky decomposition of the variance-covariance matrix elements and transformation parameters. |
data |
Data frame resulting from the 'uniformize.data.R' function. |
eoi.indicator.names |
Vector of names of the censoring indicator columns
pertaining to events of interest. Events of interest will be modeled allowing
dependence between them, whereas all censoring events (corresponding to
indicator columns not listed in |
admin |
Boolean value indicating whether the data contains administrative censoring. |
conf |
Boolean value indicating whether the data contains confounding and hence indicating the presence of Z and W. |
cf |
"Control function" to be used. This can either be the (i) estimated
control function, (ii) the true control function, (iii) the instrumental
variable, or (iv) nothing ( |
eps |
Minimum value for the diagonal elements in the covariance matrix.
Default is |
Log-likelihood evaluation for the second step in the estimation procedure.
This function defines the 'full' likelihood of the model. Specifically, it includes the estimation of the control function in the computation of the likelihood. This function is used in the estimation of the variance of the estimates (variance.cmprsk.R).
likIFG.cmprsk.Cholesky( parhatG, data, eoi.indicator.names, admin, conf, Zbin, inst )likIFG.cmprsk.Cholesky( parhatG, data, eoi.indicator.names, admin, conf, Zbin, inst )
parhatG |
The full parameter vector. |
data |
Data frame. |
eoi.indicator.names |
Vector of names of the censoring indicator columns
pertaining to events of interest. Events of interest will be modeled allowing
dependence between them, whereas all censoring events (corresponding to
indicator columns not listed in |
admin |
Boolean value indicating whether the data contains administrative censoring. |
conf |
Boolean value indicating whether the data contains confounding and hence indicating the presence of Z and W. |
Zbin |
Boolean value indicating whether the confounding variable is binary. |
inst |
Type of instrumental function to be used. |
Full model log-likelihood evaluation.
Computes the logarithm of a number.
log_transform(y)log_transform(y)
y |
Numerical value of which the logarithm is computed. |
This function returns the logarithm of the provided argument y if it is greater than zero. If y is smaller than zero, it will return 0.
This likelihood function is maximized to estimate the model parameters under the Clayton copula.
loglike.clayton.unconstrained(para, Y, Delta, Dist.T, Dist.C)loglike.clayton.unconstrained(para, Y, Delta, Dist.T, Dist.C)
para |
Estimated parameter values/initial values. |
Y |
Follow-up time. |
Delta |
Censoring indicator. |
Dist.T |
The distribution to be used for the survival time T. This argument can take one of the values from |
Dist.C |
The distribution to be used for the censoring time C. This argument can take one of the values from |
Maximized log-likelihood value.
This likelihood function is maximized to estimate the model parameters under the Frank copula.
loglike.frank.unconstrained(para, Y, Delta, Dist.T, Dist.C)loglike.frank.unconstrained(para, Y, Delta, Dist.T, Dist.C)
para |
Estimated parameter values/initial values. |
Y |
Follow-up time. |
Delta |
Censoring indicator. |
Dist.T |
The distribution to be used for the survival time T. This argument can take one of the values from |
Dist.C |
The distribution to be used for the censoring time C. This argument can take one of the values from |
Maximized log-likelihood value.
This likelihood function is maximized to estimate the model parameters under the Gaussian copula.
loglike.gaussian.unconstrained(para, Y, Delta, Dist.T, Dist.C)loglike.gaussian.unconstrained(para, Y, Delta, Dist.T, Dist.C)
para |
Estimated parameter values/initial values. |
Y |
Follow-up time. |
Delta |
Censoring indicator. |
Dist.T |
The distribution to be used for the survival time T. This argument can only the value |
Dist.C |
The distribution to be used for the censoring time C. This argument can only the value |
Maximized log-likelihood value.
This likelihood function is maximized to estimate the model parameters under the Gumbel copula.
loglike.gumbel.unconstrained(para, Y, Delta, Dist.T, Dist.C)loglike.gumbel.unconstrained(para, Y, Delta, Dist.T, Dist.C)
para |
Estimated parameter values/initial values. |
Y |
Follow-up time. |
Delta |
Censoring indicator. |
Dist.T |
The distribution to be used for the survival time T. This argument can take one of the values from |
Dist.C |
The distribution to be used for the censoring time C. This argument can take one of the values from |
Maximized log-likelihood value.
This likelihood function is maximized to estimate the model parameters under the independence copula.
loglike.indep.unconstrained(para, Y, Delta, Dist.T, Dist.C)loglike.indep.unconstrained(para, Y, Delta, Dist.T, Dist.C)
para |
Estimated parameter values/initial values. |
Y |
Follow-up time. |
Delta |
Censoring indicator. |
Dist.T |
The distribution to be used for the survival time T. This argument can take one of the values from |
Dist.C |
The distribution to be used for the censoring time C. This argument can take one of the values from |
Maximized log-likelihood value.
Change hazard and cumulative hazard to long format
Longfun(Z, T1, lhat, Lhat)Longfun(Z, T1, lhat, Lhat)
Z |
Observed survival time, which is the minimum of T, C and A, where A is the administrative censoring time. |
T1 |
Distinct observed survival time |
lhat |
Hazard function estimate |
Lhat |
Cumulative hazard function estimate |
Change a nonparametric transformation function to long format
LongNPT(Z, T1, H)LongNPT(Z, T1, H)
Z |
Observed survival time, which is the minimum of T, C and A, where A is the administrative censoring time. |
T1 |
Distinct observed survival time |
H |
Nonparametric transformation function estimate |
This function performs the maximization step in the EAM
algorithm. More specifically, it maximizes the expected improvement.
ToDo: implement sample space contractions (see comment made in documentation
of draw.sv.init).
M_step( dir, evaluations, theta.hash, fit.krige, test.fun, c, par.space, hyperparams, verbose )M_step( dir, evaluations, theta.hash, fit.krige, test.fun, c, par.space, hyperparams, verbose )
dir |
Direction to search in. |
evaluations |
Matrix containing each point that was already evaluated, alongside the corresponding test statistic and critical value, as its rows. |
theta.hash |
Tentative best value of theta. Obtained from the E-step. |
fit.krige |
Kriging model obtained from the A-step. |
test.fun |
The test function to be inverted in order to obtain the identified set. |
c |
Projection vector. |
par.space |
Bounds of the parameter space. |
hyperparams |
Parameters used in obtaining initial values
for the maximization algorithm. If |
verbose |
Verbosity parameter. |
This function obtains the vector of sample averages of each moment function.
m.bar(data, beta, t, hp, mi.mat = NULL)m.bar(data, beta, t, hp, mi.mat = NULL)
data |
Data frame. |
beta |
Vector of coefficients. |
t |
Time point at which to compute the moment functions. Also allowed to be a vector of time points (used in estimating the model under assumed time- independent coefficients). |
hp |
List of hyperparameters. |
mi.mat |
Matrix of moment function evaluations. Can be used to avoid
some computation. Default is |
Create starting values for EI maximization. Used in the M-step (get.starting.values.R).
MSpoint(draws.init)MSpoint(draws.init)
draws.init |
Initial draws. |
Bei, X. (2024). Local linearieation based subvector inference in moment inequality models. Journal of Econometrics. 238:105549-
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).
NonParTrans( resData, X, W, start = NULL, n.iter = 15, bootstrap = FALSE, n.boot = 50, eps = 0.001 )NonParTrans( resData, X, W, start = NULL, n.iter = 15, bootstrap = FALSE, n.boot = 50, eps = 0.001 )
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T |
W |
Data matrix with covariates related to C |
start |
Initial values for the finite dimensional parameters. If |
n.iter |
Number of iterations; the default is |
bootstrap |
A boolean indicating whether to compute bootstrap standard errors for making inferences. |
n.boot |
Number of bootstrap samples to use in the estimation of bootstrap standard errors if |
eps |
Convergence error. This is set by the user in such away that the desired convergence is met; the default is |
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.
Deresa, N. and Van Keilegom, I. (2021). On semiparametric modelling, estimation and inference for survival data subject to dependent censoring, Biometrika, 108, 965–979.
# 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# 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
This function normalized the covariates in the data to lie in
the unit interval based on either the empirical or known ranges of the
covariates. It is useful to perform this step when defining the instrumental
functions later on. This function is used in G.box.R, G.spline.R
and by extension in G.cd.R.
normalize.covariates( data = NULL, x = NULL, cov.ranges = NULL, idxs.c = "all", norm.cov.out = NULL, ... )normalize.covariates( data = NULL, x = NULL, cov.ranges = NULL, idxs.c = "all", norm.cov.out = NULL, ... )
data |
(optional) Data set to be used to construct the normalizing
transformation. Default is |
x |
(optional) Vector of covariates to be normalized alongside the data.
Default is |
cov.ranges |
(optional) Matrix that specifies the range of each of the
covariates in the data set. Each column corresponds to a covariate. The first
row contains the lower bound, the second row contains the upper bound.
If not supplied, the data will be normalized based on the minimum and maximum
detected values. If supplied, the non data-dependent transformation function
listed in the appendix of Andrews, Shi 2013 will be used. Default is
|
idxs.c |
(optional) Vector of indices of covariates that are continuous.
Note that that indices are relative to the covariate vector, not the full
data set. Default value is |
norm.cov.out |
(optional) The output of a previous call to this function.
Can be used to speed up computation. If both |
... |
Allows easier interchangeability between covariate normalization
functions. All arguments specified under |
Andrews, D.W.K. and Shi, X. (2013). Inference based on confitional moment inequalities. Econometrica. 81(2):609-666.
This function normalized the covariates in the data to lie in
the unit interval based on a principal component analysis. It is useful to
perform this step when defining the instrumental functions later on. This
function is used in G.box, G.spline and by extension G.cd.
normalize.covariates2( data = NULL, x = NULL, idxs.c = "all", norm.cov.out = NULL, ... )normalize.covariates2( data = NULL, x = NULL, idxs.c = "all", norm.cov.out = NULL, ... )
data |
(optional) Data set to be used to construct the normalizing
transformation. Default is |
x |
(optional) Vector of covariates to be normalized alongside the data.
Default is |
idxs.c |
(optional) Vector of indices of covariates that are continuous.
Note that that indices are relative to the covariate vector, not the full
data set. Default value is |
norm.cov.out |
(optional) The output of a previous call to this function.
Can be used to speed up computation. If both |
... |
Allows easier interchangeability between covariate normalization
functions. All arguments specified under |
This function computes the correlation matrix corresponding to
the variance-covariance matrix as returned by Sigma.hat.R
Omega.hat(Sigma)Omega.hat(Sigma)
Sigma |
The output of the function Sigma.hat |
Estimates the model parameters by maximizing the log-likelihood.
optimlikelihood(Y, Delta, Copula, Dist.T, Dist.C, start)optimlikelihood(Y, Delta, Copula, Dist.T, Dist.C, start)
Y |
Follow-up time. |
Delta |
Censoring indicator. |
Copula |
The copula family. This argument can take values from |
Dist.T |
The distribution to be used for the survival time T. This argument can take one of the values from |
Dist.C |
The distribution to be used for the censoring time C. This argument can take one of the values from |
start |
Starting values |
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
parafam.d(x, parameter, distribution, truncation = NULL)parafam.d(x, parameter, distribution, truncation = NULL)
x |
the value in which the density function will be calculated at. |
parameter |
the parameter of the specified distribution |
distribution |
the specified distribution function. |
truncation |
a positive numeric value thats denotes the value of truncation for the assumed distribution. |
Obtain the value of the distribution function
parafam.p(x, parameter, distribution, truncation = NULL)parafam.p(x, parameter, distribution, truncation = NULL)
x |
the value in which the distribution function will be calculated at. |
parameter |
the parameter of the specified distribution |
distribution |
the specified distribution function. |
truncation |
a positive numeric value thats denotes the value of truncation for the assumed distribution. |
Obtain the adjustment value of truncation
parafam.trunc(truncation, parameter, distribution)parafam.trunc(truncation, parameter, distribution)
truncation |
a positive numeric value thats denotes the value of truncation for the assumed distribution. |
parameter |
the parameter of the specified distribution |
distribution |
the specified distribution function. |
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).
ParamCop(Y, Delta, Copula, Dist.T, Dist.C, start = c(1, 1, 1, 1))ParamCop(Y, Delta, Copula, Dist.T, Dist.C, start = c(1, 1, 1, 1))
Y |
Follow-up time. |
Delta |
Censoring indicator. |
Copula |
The copula family. This argument can take values from |
Dist.T |
The distribution to be used for the survival time T. This argument can take one of the values from |
Dist.C |
The distribution to be used for the censoring time C. This argument can take one of the values from |
start |
Starting values |
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.
Czado and Van Keilegom (2023). Dependent censoring based on parametric copulas. Biometrika, 110(3), 721-738.
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)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
Parameters.Constraints(copfam, margins, cure)Parameters.Constraints(copfam, margins, cure)
copfam |
a character string that specifies the copula family. |
margins |
a list used to define the distribution structures of both the survival and censoring margins. |
cure |
a logical value that indicates whether the existence of a cured fraction should be considered. |
This function estimates bounds on the coefficients the single-
index model for the conditional cumulative
distribution function of the event time.
pi.surv( data, idx.param.of.interest, idxs.c, t, par.space, search.method = "GS", add.options = list(), verbose = 0, picturose = FALSE, parallel = FALSE )pi.surv( data, idx.param.of.interest, idxs.c, t, par.space, search.method = "GS", add.options = list(), verbose = 0, picturose = FALSE, parallel = FALSE )
data |
Data frame containing the data on which to fit the model. The columns should be named as follows: 'Y' = observed timed, 'Delta' = censoring indicators, 'X0' = intercept column, 'X1' - 'Xp' = covariates. |
idx.param.of.interest |
Index of element in the covariate vector for
which the identified interval should be estimated. It can also be specified
as |
idxs.c |
Vector of indices of the continuous covariates. Suppose the
given data contains 5 covariates, of which 'X2' and 'X5' are continuous, this
argument should be specified as |
t |
Time point for which to estimate the identified set of
|
par.space |
Matrix containing bounds on the space of the parameters. The first column corresponds to lower bounds, the second to upper bounds. The i'th row corresponds to the bounds on the i'th element in the parameter vector. |
search.method |
The search method to be used to find the identified
interval. Default is |
add.options |
List of additional options to be specified to the method.
Notably, it can be used to select the link function |
verbose |
Verbosity level. The higher the value, the more verbose the
method will be. Default is |
picturose |
Picturosity flag. If |
parallel |
Flag for whether or not parallel computing should be used.
Default is |
Matrix containing the identified intervals of the specified coefficients, as well as corresponding convergence information of the estimation algorithm.
Willems, I., Beyhum, J. and Van Keilegom, I. (2024+). Partial identification for a class of survival models under dependent censoring. (In preparation).
# 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)# 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)
This function draws the points to be evaluated.
plot_addpte(pte, col = "orange")plot_addpte(pte, col = "orange")
pte |
Vector of points to be evaluated. |
col |
Color of the points. |
This function draws evaluated points. Feasible points are indicated in green, red points correspond to infeasible points.
plot_addpte.eval(evaluations)plot_addpte.eval(evaluations)
evaluations |
Matrix of evaluations to be drawn. |
This functon draws the base plot, used when
picturose = TRUE.
plot_base(c, hp)plot_base(c, hp)
c |
Projection vector |
hp |
List of hyperparameters |
Computes a given power of a number.
power_transform(y, pw)power_transform(y, pw)
y |
The number which one wants to raise to a certain power |
pw |
The power to which to raise |
This function returns the result of raising y to the power
pw when y > 0. Otherwise, it will return 1.
The PseudoL function is maximized in order to
estimate the finite dimensional model parameters, including the dependency parameter.
This function assumes that the cumulative hazard function is known.
PseudoL(theta, resData, X, W, lhat, cumL, cop, dist)PseudoL(theta, resData, X, W, lhat, cumL, cop, dist)
theta |
Estimated parameter values/initial values for finite dimensional parameters |
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T |
W |
Data matrix with covariates related to C. First column of W should be ones |
lhat |
The estimated hazard function obtained from the output of |
cumL |
The estimated cumulative hazard function from the output of |
cop |
Which copula should be computed to account for dependency between T and C. This argument can take
one of the values from |
dist |
The distribution to be used for the dependent censoring C. Only two distributions are allowed, i.e, Weibull
and lognormal distributions. With the value |
maximized log-likelihood value
This function computes the loss function at a given point.
S.func(m, Sigma)S.func(m, Sigma)
m |
Vector of averages of moment functions. |
Sigma |
Sample variance-covariance matrix of moment functions. |
S(m, Sigma).
This function computes the score vectors and the Jacobean matrix for finite model parameters.
ScoreEqn(theta, resData, X, W, H)ScoreEqn(theta, resData, X, W, H)
theta |
Vector of parameters in the semiparametric transformation model. |
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T. |
W |
Data matrix with covariates related to C. |
H |
The estimated non-parametric transformation function for a given value of theta |
Function to indicate position of t in observed survival time
SearchIndicate(t, T1)SearchIndicate(t, T1)
t |
fixed time t |
T1 |
distinct observed survival time |
This function returns a list with the (default) hyperparameters used in the EAM algorithm
set.EAM.hyperparameters(options)set.EAM.hyperparameters(options)
options |
A list of user-specified values for (some of) the hyperparameters. These hyperparameters can include:
|
List of hyperparameters for the EAM algotithm.
This function returns a list with the (default) hyperparameters used in the grid search algorithm
set.GS.hyperparameters(options)set.GS.hyperparameters(options)
options |
A list of user-specified values for (some of) the hyperparameters. These hyperparameters could include:
|
List of hyperparameters for the gridsearch and binary search algorithms.
This function defines all the necessary hyperparameters used to run the methodology.
set.hyperparameters(data, par.space, c, search.method, options)set.hyperparameters(data, par.space, c, search.method, options)
data |
Data frame. |
par.space |
Bounds on the parameter space. |
c |
Projection vector. |
search.method |
Search method to use ( |
options |
List of user specified hyperparameters that will substitute the corresponding default values. This list can contain the entries:
Other (hidden) options can also be overwritten, though we highly discourage
this. If necessary, you can consult the source code of this functions to
find the names of the desired parameters and add their name alongside their
desired value as an entry in |
The list of hyperparameters.
This function comptutes the empricical variance-covariance matrix of the moment functions.
Sigma.hat(data, beta, t, hp, m.avg = NULL, mi.mat = NULL)Sigma.hat(data, beta, t, hp, m.avg = NULL, mi.mat = NULL)
data |
Data frame. |
beta |
Coefficient vector. |
t |
Time point of interest. |
hp |
List of hyperparameters. |
m.avg |
A precomputed vector of the sample average of the moment
functions. If not supplied, this vector is computed. Default is
|
mi.mat |
A precomputed matrix of moment function evaluations at each
observation. If supplied, some computations can be skipped. Default is
|
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).
SolveH(theta, resData, X, W)SolveH(theta, resData, X, W)
theta |
Vector of parameters in the semiparametric transformation model. |
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T. |
W |
Data matrix with covariates related to C. |
Returns the estimated transformation function H for a fixed value of parameters theta.
Deresa, N. and Van Keilegom, I. (2021). On semiparametric modelling, estimation and inference for survival data subject to dependent censoring, Biometrika, 108, 965–979.
This function obtains an estimating equation of H at the first observed survival time t1.
SolveHt1(Ht1, Z, nu, t, X, W, theta)SolveHt1(Ht1, Z, nu, t, X, W, theta)
Ht1 |
The solver solves for an optimal value of Ht1 by equating the estimating equation to zero. |
Z |
The observed survival time, which is the minimum of T, C and A. |
nu |
The censoring indicator for T or C |
t |
A fixed time point |
X |
Data matrix with covariates related to T. |
W |
Data matrix with covariates related to C. |
theta |
Vector of parameters |
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.
SolveL( theta, resData, X, W, cop = c("Frank", "Gumbel", "Normal"), dist = c("Weibull", "lognormal") )SolveL( theta, resData, X, W, cop = c("Frank", "Gumbel", "Normal"), dist = c("Weibull", "lognormal") )
theta |
Estimated parameter values/initial values for finite dimensional parameters |
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T |
W |
Data matrix with covariates related to C. First column of W should be ones |
cop |
Which copula should be computed to account for dependency between T and C. This argument can take
one of the values from |
dist |
The distribution to be used for the dependent censoring C. Only two distributions are allowed, i.e, Weibull
and lognormal distributions. With the value |
This function returns an estimated hazard function, cumulative hazard function and distinct observed survival times;
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")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")
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.
SolveLI(theta, resData, X)SolveLI(theta, resData, X)
theta |
Estimated parameter values/initial values for finite dimensional parameters |
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T |
This function returns an estimated hazard function, cumulative hazard function and distinct observed survival times;
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")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")
This function estimates the model parameters
SolveScore(theta, resData, X, W, H, eps = 0.001)SolveScore(theta, resData, X, W, H, eps = 0.001)
theta |
Vector of parameters in the semiparametric transformation model. |
resData |
Data matrix with three columns; Z = the observed survival time, d1 = the censoring indicator of T and d2 = the censoring indicator of C. |
X |
Data matrix with covariates related to T. |
W |
Data matrix with covariates related to C. |
H |
The estimated non-parametric transformation function for a given value of theta. |
eps |
Convergence error. |
depCensoringFit objectSummary of depCensoringFit object
## S3 method for class 'depFit' summary(object, ...)## S3 method for class 'depFit' summary(object, ...)
object |
Output of |
... |
Further arguments |
Summary of dependent censoring model fit in the form of table
indepCensoringFit objectSummary of indepCensoringFit object
## S3 method for class 'indepFit' summary(object, ...)## S3 method for class 'indepFit' summary(object, ...)
object |
Output of |
... |
Further arguments |
Summary of independent censoring model fit in the form of table
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.
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) )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) )
yobs |
a numeric vector that indicated the observed survival times. |
delta |
a numeric vector that stores the right-censoring indicators. |
tm |
a numeric vector that contains interested non-negative time points at which the survival probabilities will be evluated.
Note that if we omit the definition of this argument (the default value becomes |
copfam |
a character string that specifies the copula family.
Currently, it supports Archimedean copula families, including |
margins |
a list used to define the distribution structures of both the survival and censoring margins. Specifically, it contains the following elements:
Note if one of the marginal distributions should be modeled nonparametrically, one can let the corresponding argument to be |
cure |
a logical value that indicates whether the existence of a cured fraction should be considered. |
Var |
a list that controls the execution of the bootstrap for variance estimation,
and it contains two elements:
|
control |
indicates more detailed control of the underlying model fitting procedures. It is a list of the following three arguments:
|
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.
@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).
A list of fitted results is returned. Within this outputted list, the following elements can be found:
probssurvival probabilities of the survial margin at tm.
ktauKendall's tau.
paraparestimation of all parameters (except Kendall's tau) contained in the parametric part.
GoFgoodness-of-test results.
cureratecure rate. If cure = FALSE, it is NULL.
#----------------------------------------------------------# # 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#----------------------------------------------------------# # 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
SurvDC.GoF( yobs, delta, copfam, margins, ktau, parapar, cure = FALSE, curerate = NULL )SurvDC.GoF( yobs, delta, copfam, margins, ktau, parapar, cure = FALSE, curerate = NULL )
yobs |
a numeric vector that indicated the observed survival times. |
delta |
a numeric vector that stores the right-censoring indicators. |
copfam |
a character string that specifies the copula family. |
margins |
a list used to define the distribution structures of both the survival and censoring margins. |
ktau |
Kendall's tau. |
parapar |
parametric parameters. |
cure |
a logical value that indicates whether the existence of a cured fraction should be considered. |
curerate |
value of cure rate. |
Estimated survival function based on copula-graphic estimator (Archimedean copula only)
SurvFunc.CG(tm = NULL, yobs, delta, copfam, ktau, coppar = NULL)SurvFunc.CG(tm = NULL, yobs, delta, copfam, ktau, coppar = NULL)
tm |
a vector contains all time points that the survival function will be calculated at. |
yobs |
a numeric vector that indicated the observed survival times. |
delta |
a numeric vector that stores the right-censoring indicators. |
copfam |
a character string that denotes the copula family. |
ktau |
a numeric value that denotes the Kendall's tau. |
coppar |
a numeric value that denotes the copula parameter. |
Estimated survival function based on Kaplan-Meier estimator
SurvFunc.KM(tm = NULL, yobs, delta, type = "right")SurvFunc.KM(tm = NULL, yobs, delta, type = "right")
tm |
a vector contains all time points that the survival function will be calculated at. |
yobs |
a numeric vector that indicated the observed survival times. |
delta |
a numeric vector that stores the right-censoring indicators. |
type |
a character string that specifies the type of the step function. If |
Maximum likelihood estimator for a given parametric distribution
SurvMLE( yobs, delta, distribution, truncation = NULL, cure = FALSE, maxit = 300 )SurvMLE( yobs, delta, distribution, truncation = NULL, cure = FALSE, maxit = 300 )
yobs |
a numeric vector that indicated the observed survival times. |
delta |
a numeric vector that stores the right-censoring indicators. |
distribution |
the specified distribution function. |
truncation |
a positive numeric value thats denotes the value of truncation for the assumed distribution. |
cure |
a logical value that indicates whether the existence of a cured fraction should be considered. |
maxit |
a positive integer that denotes the maximum iteration number in optimization. |
Likelihood for a given parametric distribution
SurvMLE.Likelihood( param, yobs, delta, distribution, truncation = NULL, cure = FALSE )SurvMLE.Likelihood( param, yobs, delta, distribution, truncation = NULL, cure = FALSE )
param |
a vector contains all parametric parameters. |
yobs |
a numeric vector that indicated the observed survival times. |
delta |
a numeric vector that stores the right-censoring indicators. |
distribution |
the specified distribution function. |
truncation |
a positive numeric value thats denotes the value of truncation for the assumed distribution. |
cure |
a logical value that indicates whether the existence of a cured fraction should be considered. |
Generates the follow-up time and censoring indicator according to the specified model.
TCsim( tau = 0, Copula = "frank", Dist.T = "lnorm", Dist.C = "lnorm", par.T = c(0, 1), par.C = c(0, 1), n = 10000 )TCsim( tau = 0, Copula = "frank", Dist.T = "lnorm", Dist.C = "lnorm", par.T = c(0, 1), par.C = c(0, 1), n = 10000 )
tau |
Value of Kendall's tau for (T,C). The default value is 0. |
Copula |
The copula family. This argument can take values from |
Dist.T |
Distribution of the survival time T. This argument can take one of the values from |
Dist.C |
Distribution of the censoring time C. This argument can take one of the values from |
par.T |
Parameter values for the distribution of T. |
par.C |
Parameter values for the distribution of C. |
n |
Sample size. |
A list containing the generated follow-up times and censoring indicators.
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)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)
This function performs the unconditional moment restriction test as described in Bei (2024).
test.point_Bei( r, c, t, par.space, data, hp, verbose = FALSE, inst.func.evals = NULL, alpha = 0.95, parallel = FALSE )test.point_Bei( r, c, t, par.space, data, hp, verbose = FALSE, inst.func.evals = NULL, alpha = 0.95, parallel = FALSE )
r |
Result of the projection for which the test should be carried out. |
c |
The projection matrix. For now, c is restricted to being an elementary vector, i.e. c = (0, ...,0, 1, 0, ..., 0). |
t |
The time point at which to evaluate theta. |
par.space |
Matrix containing 2 columns and |
data |
Data frame on which to base the test. |
hp |
List of hyperparameters needed. |
verbose |
Boolean variable indicating whether to print updates of the estimation process to the console. |
inst.func.evals |
Matrix of precomputed instrumental function
evaluations for each observation in the data set. Used to speed up the
simulations. If |
alpha |
The significance level at which to perform the test. Default is
|
parallel |
Flag for whether or not parallel computing should be used.
Default is |
Bei, X. (2024). Local linearization based subvector inference in moment inequality models. Journal of Econometrics, 238(1), 105549-. https://doi.org/10.1016/j.jeconom.2023.10554
This function performs the unconditional moment restriction test
as described in Bei (2024). This function directly extends
test.point_Bei by allowing for pairs of moment restrictions over a
grid of time points.
test.point_Bei_MT( r, c, t, par.space, data, hp, verbose = FALSE, inst.func.evals = NULL, alpha = 0.95, parallel = FALSE )test.point_Bei_MT( r, c, t, par.space, data, hp, verbose = FALSE, inst.func.evals = NULL, alpha = 0.95, parallel = FALSE )
r |
Result of the projection for which the test should be carried out. |
c |
The projection matrix. For now, c is restricted to being an elementary vector, i.e. c = (0, ...,0, 1, 0, ..., 0). |
t |
The time point at which to evaluate theta. Also allowed to be a vector of time points (used in estimating the model under assumed time- independent coefficients). |
par.space |
Matrix containing 2 columns and |
data |
Data frame on which to base the test. |
hp |
List of hyperparameters needed. |
verbose |
Boolean variable indicating whether to print updates of the estimation process to the console. |
inst.func.evals |
Matrix of precomputed instrumental function
evaluations for each observation in the data set. Used to speed up the
simulations. If |
alpha |
The significance level at which to perform the test. Default is
|
parallel |
Flag for whether or not parallel computing should be used.
Default is |
Bei, X. (2024). Local linearization based subvector inference in moment inequality models. Journal of Econometrics, 238(1), 105549-. https://doi.org/10.1016/j.jeconom.2023.10554
Checks the required preconditions of the data and possibly restructures the data.
uniformize.data( data, admin = FALSE, conf = FALSE, comp.risks = FALSE, Zbin = NULL, Wbin = NULL )uniformize.data( data, admin = FALSE, conf = FALSE, comp.risks = FALSE, Zbin = NULL, Wbin = NULL )
data |
A data frame that should contain columns named |
admin |
Boolean value indicating whether the provided data frame contains
administrative (i.e. independent) censoring on top of the dependent censoring
(in the column named |
conf |
Boolean value indicating whether the provided data frame contains
a confounded variable and a corresponding instrument. If |
comp.risks |
Boolean value indicating whether the provided data frame
contains competing risks. If |
Zbin |
Boolean or integer value (0, 1) indicating whether the confounded
variable is binary. |
Wbin |
Boolean or integer value (0, 1) indicating whether the instrument
is binary. |
Returns the uniformized data set.
This function computes the variance of the estimates computed by the 'estimate.cmprsk.R' function.
variance.cmprsk( parhatc, gammaest, data, admin, conf, inst, cf, eoi.indicator.names, Zbin, use.chol, n.trans, totparl )variance.cmprsk( parhatc, gammaest, data, admin, conf, inst, cf, eoi.indicator.names, Zbin, use.chol, n.trans, totparl )
parhatc |
Vector of estimated parameters, computed in the first part of
|
gammaest |
Vector of estimated parameters in the regression model for the control function. |
data |
A data frame. |
admin |
Boolean value indicating whether the data contains administrative censoring. |
conf |
Boolean value indicating whether the data contains confounding
and hence indicating the presence of |
inst |
Variable encoding which approach should be used for dealing with
the confounding. |
cf |
The control function used to estimate the second step. |
eoi.indicator.names |
Vector of names of the censoring indicator columns
pertaining to events of interest. Events of interest will be modeled allowing
dependence between them, whereas all censoring events (corresponding to
indicator columns not listed in |
Zbin |
Indicator value indicating whether ( |
use.chol |
Boolean value indicating whether the cholesky decomposition was used in estimating the covariance matrix. |
n.trans |
Number of competing risks in the model (and hence, number of transformation models). |
totparl |
Total number of covariate effects (including intercepts) in all of the transformation models combined. |
Variance estimates of the provided vector of estimated parameters.
Computes the Yeo-Johnson transformation of the provided argument.
YJtrans(y, theta)YJtrans(y, theta)
y |
The argument to be supplied to the Yeo-Johnson transformation. |
theta |
The parameter of the Yeo-Johnson transformation. This should be a number in the range [0,2]. |
The transformed value of y.