1
$\begingroup$

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)
# 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)[[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))
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.