For this project, standard Cox lasso with cross-validation has an
unacceptable computational burden, while fitting the full unpenalized
Cox model with coxph() is feasible. These are established
project assumptions. The challenge is the repeated optimization needed
to examine different penalty strengths across the cross-validation
training samples.
The scientific questions remain the familiar ones addressed by Cox lasso: which measured patient characteristics should be retained in a model of time to an event, and which amount of shrinkage gives the best predictions for held-out patients? In claims data, the candidates might include diagnoses, medications, and healthcare utilization. Adaptive LSA addresses these questions by producing a selected set of predictors and their Cox coefficients. Selection here concerns associations and prediction, rather than establishing causal effects.
The computational shortcut starts with the full Cox fit. Its coefficients and non-robust covariance matrix provide a compact summary from which adaptive LSA constructs an approximation to the original fitting objective. Candidate models at different penalty strengths are then obtained from this summary. This avoids repeatedly optimizing the original Cox objective throughout the penalty path. Each cross-validation training sample needs one full Cox fit; the existing full-data fit is reused for the final model. Validation still evaluates predictions on the actual held-out patients.
The result serves the same variable-selection and prediction purposes as Cox lasso, using a more practical fitting workflow for this project. Because it uses an approximation and adaptive penalties, it need not produce exactly the same coefficients or selected variables as a direct Cox lasso fit.
Implement unstratified Cox regression + adaptive LSA + model-based covariance + cross-validation for an analysis-ready RWE claims cohort with more than one million patients and approximately 100 candidate coefficients.
A full unpenalized coxph() fit is feasible and is the
starting point. Use the existing full-data fit for the final LSA model.
Use separate training-fold fits during cross-validation.
| Step | R function or object |
|---|---|
| Fit the full candidate Cox model | survival::coxph() |
| Extract initial coefficients and non-robust covariance | coef(fit) and vcov(fit), with
robust=FALSE |
| Solve the adaptive LSA objective | glmnet::glmnet(), configured below |
| Evaluate held-out Cox predictions | glmnet::coxnet.deviance() |
| Select the tuning parameter | A custom five-fold loop around these functions |
Use the same tie setting as the accepted full-data Cox fit for all
training fits and validation evaluations. Install a glmnet
version whose coxnet.deviance() accepts the corresponding
cox.ties argument. Record the package versions.
For any training dataset \(A\), let \(m_A\) be its number of patients, \(b_A\) the full-model Cox coefficient estimate, and \(V_A\) its non-robust covariance matrix. Define
\[ H_A=\frac{V_A^{-1}}{m_A},\qquad w_{Aj}=\frac{1}{\max(|b_{Aj}|,\epsilon_A)},\qquad q_{Aj}=\frac{w_{Aj}}{p^{-1}\sum_{\ell=1}^{p}w_{A\ell}}. \]
Use the usual adaptive exponent \(\gamma=1\). Set \(\epsilon_A=1/m_A\) as a small numerical floor and record whether it is used. Keep coefficients, covariance, and predictions in the same covariate parameterization.
For each candidate \(\lambda\), compute
\[ \widehat\beta_A(\lambda)= \arg\min_\beta \left[ \frac12(\beta-b_A)^\top H_A(\beta-b_A) +\lambda\sum_{j=1}^{p}q_{Aj}|\beta_j| \right]. \tag{1} \]
This is the adaptive LSA objective associated with Wang and Leng (2007). The normalization of \(q_A\) changes the numerical scale of \(\lambda\), but not the available models along the complete path. All candidate coefficients are penalized.
Interpretation: the quadratic term keeps the estimate near the initial Cox fit, while the adaptive penalty shrinks coefficients and sets some to zero. Coefficients that are small in the initial fit receive stronger penalties.
Use the patient-level survival outcome, candidate covariates, and patient identifier from the agreed claims-data cohort. Preserve the candidate-column names and order.
Create five folds, balancing event counts and using a fixed random seed. Each patient belongs to one fold. Estimate any data-dependent preprocessing within the training fold and apply it to that fold’s validation patients.
Use these defaults:
| Setting | Value |
|---|---|
| Number of folds | 5 |
| Adaptive exponent | 1 |
| Candidate tuning fractions | 100 logarithmically spaced values from \(1\) to \(10^{-6}\) |
| Additional candidate | The unpenalized fit, represented by fraction \(0\) |
| Selection rule | Lowest mean held-out validation loss |
For fold \(k\), let \(A_k\) contain the other four folds and \(F_k\) the held-out patients.
Call survival::coxph() on \(A_k\) using the complete candidate formula,
robust=FALSE, singular.ok=FALSE, and the
agreed tie setting. Extract \(b_{A_k}\)
and \(V_{A_k}\) and construct Equation
(1).
Require convergence, finite coefficients, and a positive-definite covariance matrix. Retain the complete covariance matrix.
Recompute the initial coefficients, covariance, and adaptive weights in every training fold. Do not derive these quantities from the full-data fit for use inside cross-validation.
With an existing full-data fit, this requires five additional Cox fits. Cache each training-fold summary so that changing the tuning grid does not require repeating those fits.
Use a common fraction \(r\) to identify corresponding candidates:
\[ \lambda_{\max,A} = \max_j\frac{|(H_A b_A)_j|}{q_{Aj}}, \qquad \lambda_A(r)=r\lambda_{\max,A}. \tag{2} \]
Thus, a fraction identifies the same relative penalty level in every fold, while each fold has its own actual \(\lambda\) values. Use the same fraction grid for final full-data fitting.
At \(r=0\), use \(b_A\) directly. If \(\lambda_{\max,A}=0\), the path is the all-zero coefficient vector.
Prepare the inputs for glmnet():
R = chol(H_A).Xstar = sqrt(p) * R.ystar = drop(Xstar %*% b_A).Use the following arguments:
| Argument | Value |
|---|---|
x, y |
Xstar, ystar |
family |
"gaussian" |
alpha |
1 |
intercept |
FALSE |
standardize |
FALSE |
penalty.factor |
\(q_A\) |
lambda |
The explicit sequence \(\lambda_A(r)\) |
These settings solve Equation (1). Keep the exact matrix orientation and scaling above.
Retain the returned coefficient matrix, preserving coefficient names. Ensure that the solver returns every requested positive tuning value; adjust its documented path-stopping controls if necessary. Append \(b_A\) for \(r=0\).
For each candidate fraction, obtain the validation linear predictors
\[ \eta_k(r)=X_{F_k}\widehat\beta_{A_k}\{\lambda_{A_k}(r)\}. \]
Here \(X_{F_k}\) uses the preprocessing learned from \(A_k\).
Call glmnet::coxnet.deviance() with the validation
predictions, the validation Surv response, the agreed
cox.ties setting, and std.weights=FALSE.
Evaluate the null prediction vector once on the same validation
data.
Let \(D_k(\eta)\) denote the returned deviance and \(d_k\) the validation event count. Calculate
\[ L_k(r)=\frac{D_k\{\eta_k(r)\}-D_k(0)}{d_k}. \tag{3} \]
Lower values are better. The null subtraction does not change which candidate minimizes the validation loss.
This uses held-out-fold partial likelihood: pass the entire validation fold to the scoring function. Every fold must contain enough events for meaningful validation. If predictions are evaluated in batches, batch candidate models while retaining all validation patients in each scoring call.
Use the existing scoring function; no custom risk-set or tie-handling implementation is required.
Calculate
\[ \overline L(r)=\frac15\sum_{k=1}^{5}L_k(r), \qquad r_*=\arg\min_r\overline L(r). \tag{4} \]
For numerical ties, choose the larger fraction. If the minimum is at the smallest positive fraction, refine the grid near that endpoint using the cached training summaries.
The result may be \(r_*=0\). In that case, cross-validation has selected the unpenalized model.
Cross-validation must split patients and repeat Steps 2–5. Do not
call cv.glmnet() on Xstar and
ystar: those constructed rows are not patients.
From the accepted full-data coxph() fit, extract the
full-data coefficients and non-robust covariance. Recalculate \(H\), \(q\), and \(\lambda_{\max}\) using all analysis
patients.
Set \(\lambda_* = r_*\lambda_{\max}\) and solve Equation (1) using the same solver configuration. If \(r_*=0\), return the initial full-model coefficients.
The final model is this penalized adaptive-LSA estimate. Its selected variables are the coefficients estimated as nonzero. Return coefficients in the original covariate units, together with the preprocessing needed for prediction.
Return:
Confirm that all pilot fits used non-robust covariance, every validation fold was excluded from its training fit and adaptive weights, and the final tuning value was selected by Equation (4).
Wang, H., and Leng, C. (2007). Unified LASSO Estimation by Least Squares Approximation. Journal of the American Statistical Association, 102(479), 1039–1048. DOI.
Used for the adaptive LSA objective in Equation (1). Under suitable model assumptions and penalty choices, adaptive LSA can identify the truly relevant variables and estimate their effects as accurately, in large samples, as if the relevant variables had been known beforehand. Cross-validation selects the penalty for prediction; it does not automatically guarantee this theoretical variable-selection result.
Dai, B., and Breheny, P. (2024). Cross-Validation Approaches for Penalized Cox Regression. Statistical Methods in Medical Research, 33(4). DOI.
Used for held-out-fold partial-likelihood validation. Equations (3)–(4) specify this project’s event normalization and fold averaging; Equation (2) specifies how its adaptive-LSA paths are aligned.
Function documentation: