pcoxtime: Fit penalized Cox model

Description

Fits a Cox model with either lasso, ridge or elasticnet penalty for both time-independent and time-dependent (varying) covariates survival data.

Usage

pcoxtime(
  formula,
  data,
  alpha = 1,
  lambda = 1,
  maxiter = 1e+05,
  tol = 1e-08,
  quietly = FALSE,
  lambmax = FALSE,
  origin_scale = TRUE,
  contrasts.arg = NULL,
  xlevs = NULL,
  na.action = na.omit,
  ...
)

Arguments

formula

object of class formula describing the model. The response is specified similar to Surv function from package survival. The terms (predictors) are specified on the right of "~" in the formula.

data

optional data frame containing variables specified in the formula.

alpha

elasticnet mixing parameter, with $0 \le\alpha\le 1$. See details

lambda

tuning parameter for the lasso penalization, with $\lambda \ge 0$. $\lambda = 0$ fits unpenalized Cox model. See details

maxiter

maximum number of iterations to convergence. Default is $1e4$. Consider increasing it if the model does not converge.

tol

convergence threshold for proximal gradient gradient descent. Each proximal update continues until the relative change in all the coefficients (i.e. $\sqrt{\sum(\beta_{k+1} - \beta_k)^2}$/stepsize) is less than tol. The default value is $1e-8$.

quietly

logical. If TRUE, iteration progress printed.

lambmax

logical. Sufficiently large, $\lambda_{\max}$, that sets $\beta = 0$ for regularization path. If TRUE, $\lambda_{\max}$ is returned.

origin_scale

logical. If TRUE (default), the estimated coefficients are returned on the original covariate scale. Otherwise, FALSE, coefficients are standardized.

contrasts.arg

an optional list. See the contrasts.arg of model.matrix.default.

xlevs

a named list of character vectors giving the full set of levels to be assumed for each factor. See model.frame.

na.action

a function which indicates what should happen when the data contain NAs. See model.frame.

...

additional arguments not implemented.

Value

An S3 object of class pcoxtime:

coef

a named vector of coefficients. If any of the coefficients violates KKT conditions, the model will print a warning but still return coefficient estimates.

min.nloglik

estimated log-likelihood at convergence.

min.dev

the deviation satisfying the tol stopping creteria.

iter.dev

deviations between previous and current coefficient estimate at each iteration.

convergence

convergence message containing the number of iterations

the number of observations used in the fit.

n.risk

the number of individuals at risk at time t.

n.event

the number of events that occur at time t.

n.censor

the number of subjects who exit the risk set, without an event, at time t.

time

time points defined by the risk set.

Surv object defining the event times and event status.

data

data frame used.

timevarlabel, eventvarlabel

time and event variables, respectively.

predictors

a vector of predictors/covariates in the model.

lambda, alpha

lambda and alpha used, respectively.

formula

model formula used in the fit.

means

vector of column means of the X matrix. Subsequent survival curves are adjusted to this value.

assign, xlevels, terms

See model.frame for assign, xlevels, contrasts and terms.

Details

The algorithm estimates the coefficients based on observed survival data, with either time-independent or time-dependent covariates, through penalized partial log-likelihood

$$\mathrm{pen} ~ \mathit{\ell(\beta)_{\alpha, \lambda}} = -\mathit{\ell(\beta)} + P_{\alpha, \lambda}(\beta)$$

using elasticnet (which combines both lasso and ridge) penalty

$$\lambda\left(\alpha\sum_{i=1}^p|\beta_i| + 0.5(1 - \alpha)\sum_{i=1}^p\beta_i^2 \right)$$.

alpha = 1 ($\alpha$) is the lasso penalty, and alpha = 0 is the ridge penalty. lambda = 0 fits the standard Cox proportional hazard model.

User can provide a particular lambda. Typical usage is to use the pcoxtimecv to select the optimal lambda first.

The routine to handle time-dependent covariates is similar to that implemented in coxph: if there are tied event times, Breslow approximation is used.

Examples

Run this code

# NOT RUN {
# Time-independent covariates
if (packageVersion("survival")>="3.2.9") {
	data(cancer, package="survival")
} else {
	data(veteran, package="survival")
}
## Fit unpenalized Cox using pcoxtime
lam <- 0 # Should fit unpenalized Cox model
pfit1 <- pcoxtime(Surv(time, status) ~ factor(trt) + karno + diagtime + age + prior
	, data = veteran
	, lambda = lam 
	, alpha = 1
)
print(pfit1)

## fit survival::coxph
cfit1 <- coxph(Surv(time, status) ~ factor(trt) + karno + diagtime + age + prior
	, data = veteran
	, method = 'breslow'
	, ties = "breslow"
)
print(cfit1)

## Penalized Cox model (pcoxtime) 
lam <- 0.1
alp <- 0.5
pfit2 <- pcoxtime(Surv(time, status) ~ factor(trt) + karno + diagtime + age + prior
	, data = veteran
	, lambda = lam 
	, alpha = alp
)
print(pfit2)

# Time-varying covariates
data(heart, package="survival")
lam <- 0.1
alp <- 0.8
pfit2 <- pcoxtime(Surv(start, stop, event) ~ age + year + surgery + transplant
	, data = heart
	, lambda = lam
	, alpha = alp
)
print(pfit2)

# }