I am trying to find a variable signature associated with a characteristic. Particularly I am looking to get a prognostic model from multi-variable data for gene expression. I have the "Time (survival years)" and "Status (dead/alive)" data for individuals and I also have expression data for those individuals.
Based on my reading about statistical models and posts 1 , 2 and 3, I decided to use multivariate Cox propotional-hazard model. I was reading how to select few best markers, and I found "backward" selection and "LASSO" method to achieve that. Further reading suggested "LASSO" could be a good choice and can be implemented using glmnet package in R.
I divided my data into 60-40 (TRAIN-TEST) proportions, and using the 60% to find those variables and then use them on 40% to find there success as prognostic markers.
The training data has 457 subjects (patients) and 612 covariates (continuous log transformed values). The survivability data is right censored with ~15% Dead cases (Status=1). I don't have any idea about what would be a good number of predictors.
For the TRAIN data, I used the code below.
md #expression data mdsurv #surviablity data( years) and death (0 or 1) #creating table eventT<-mdsurv$OS obsLen<-5 # right censored 5 years censT<-rep(obsLen,365) obsT<-pmin(eventT,censT) status<-eventT<=censT mdsurv$Status<-as.numeric(status) keep<-mdsurv$Dead==0 & mdsurv$Status==1 mdsurv$Status[keep]<-0 cvfit<-cv.glmnet(mdt,mdsurv, family = "cox",alpha=1) plot(cvfit) c<-coef(cvfit ,s='lambda.min')
I am getting 16 genes(covariates) with lambda.min cutoff for coef() and zero with lambda.1se. When I do univariate cox analysis. HR value for some of those genes is less than one. Shall I still use those genes for my test subject or just make a model using the one with HR>=1.
Thanks. P.S. Edited with new information.More edits