How to conduct stratified Cox regression with Lasso regularization and k-fold cross validation in glmnet in R?

I have been conducting Cox PH survival analysis and I am new to this. I am following the steps below:

Create surv object and predictor model matrix for variables I am interested in. Since my data is small (~700 obs, 200 explanatory variables, some are highly correlated >0.80), I am running 5-fold CV with L1 regularization and simulate over a sequence of lambda penalization term. Pick lambda minimizing cross validated mean error, and then pick the associated non-zero coefficient variables as meaningful/important variables. Run bi-directional step-wise model to pick a final model with statistical significance. Then I check if the model violates PH assumptions or not with cox.zph(). While doing these, I found out that one of my binary variables' chisq from cox.zph test comes out very small (X15 below).To my understanding, it means it is violating the assumption of proportional hazard due to some sort of interaction or time variant.

model1 <- coxph(Surv(time, status)~ X1+X2+X15, dat=df)

# Result of cox.zph for model 1
print(cox.zph(model1))
chisq   df    p
X1       2.74001  1 0.0979
X2       0.00073  1 0.9785
X15      6.70160  1 0.0096
GLOBAL   6.94504  3 0.0737

# model1 is a coxph object
AIC(model1)
3803.522

# concordance
model1$concordance[[6]] 0.579523  One of the solutions suggested is to use this particular variable as a strata. When I use and re-estimate model, my AIC decreases significantly (by about 485), and anova gives a very small p-value. This makes me believe that my model with strata is a better predictor. However, my concordance is decreased by a lot (0.58 to 0.55) when I use strata. What I mean is similar to the following: model2 <- coxph(Surv(time, status)~ X1+X2+strata(X15), dat=df) # Result of cox.zph for model 1 print(cox.zph(model2)) chisq df p X1 0.1926 1 0.66 X2 0.0957 1 0.76 GLOBAL 0.2603 2 0.88 # model2 is a coxph object AIC(model2) 3318.471 # concordance model2$concordance[[6]]
0.5507666

# anova
anova(model1, model2)

Analysis of Deviance Table
Cox model: response is  Surv(time, status)
Model 1: ~ X1 + X2 + X15
Model 2: ~ X1 + X2 + strata(X15)
loglik  Chisq Df P(>|Chi|)
1 -1898.8
2 -1657.2 483.05  1 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


I am not sure what exactly I should do as I discovered this after cross validation and lasso step or how I should move forward with this since AIC decreases and concordance as well decreasing but anova says model2 is better. I was wondering if the reason is because I do not have any strata in my cross-validation step. If so, I was wondering how I can add my binary variable as a strata to the k-fold cv with lasso. Also, is adding a strata variabla the same thing as stratified k-fold cross validation? I am a little confused with the term strata here and stratified cross validation. I have found out balancedFolds from c060 package in R but I am not sure if this is what I need.

I tried to generate a synthetic data similar to the data for a reproducible example for community as I cannot share the original due to confidentiality. However, any synthetic data I generated did not exhibit the exact similar situation I tried to display above, neither the data I tried to generate with runif/rnorm. Regardless, I am sharing a piece of my code regarding what I was trying to do.

# necessary libraries
library(survival)
library(glmnet)
# set seed for reproducible results
set.seed(12345)
# Step 1) Create Surv Object and predictor matrix
pred_vars <- as.matrix(dat[, 3:206])
surv_obj <- Surv(dat$$time, dat$$status)
# Step 2) Create lambda sequence to try for Lasso
lambdas_to_try <- 10^seq(-5,5,length.out = 500)
# Step 3) Set k-fold CV with 5 folds, alpha=1 for Lasso
cv <- cv.glmnet(pred_vars, surv_obj, alpha=1, family="cox", lambda=lambdas_to_try, standardize=T, nfolds=5)
plot(cv)
dim(cv$$glmnet.fit$$beta)
# Step 4) Pick best nonzero coefficients from min lambda
min_lambda <- cv$$lambda.min best_coefs <- cv$$glmnet.fit$$beta[1:dim(cv$$glmnet.fit$$beta)[[1]],which(cv$$lambda==cv$lambda.min)] nonzero_best_coefs <- abs(best_coefs)>0 best_coefs <- best_coefs[nonzero_best_coefs] # Step 5) Run stepwise from meaningful variables based on AIC to pick best model selected_vars <- pred_vars[,nonzero_best_coefs] new_dat <-as.data.frame(cbind(surv_obj, selected_vars)) reduced_mod <- step(coxph(Surv(time,status) ~ .,data=new_dat),direction="both") summary(reduced_mod) # Models m1 <- coxph(Surv(time, status) ~ X1 + X4, data=new_dat) AIC(m1) m1$concordance
print(cox.zph(m1))
plot(cox.zph(m1))

m2 <- coxph(Surv(time, status) ~ X1 + X4 + X15, data=new_dat)
AIC(m2)
m2$concordance print(cox.zph(m2)) plot(cox.zph(m2)) m3 <- coxph(Surv(time, status) ~ X1 + X4 + strata(X15), new_dat) AIC(m3) m3$concordance
print(cox.zph(m3))
plot(cox.zph(m3))