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[] 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[] 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) # Read data dat <- read.table("./dat.txt") # 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)[],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))