I am a student conducting a survival analysis in R. I am using 6,000 genes from 249 patients each, and am testing each gene separately by putting them in an individual Cox regression model. Here is the code that I have written in order to be able to test them all at the same time instead of creating 6,000 different models (for reference, EFST is the event free survival time and Status is the vital state after the first event):
library(survival) dataOne <- read.delim("/Users/menon/OneDrive/Desktop/csrsef files/FirstSplitDataFrame.txt") myFinalData <- dataOne[, -1] EFST <- myFinalData[ , c("EFST")] Status <- myFinalData[, c("Status")] formulas <- sapply(names(myFinalData)[3:ncol(myFinalData)], function(x) as.formula(paste ('Surv(EFST, Status) ~ ',x))) models <- lapply(formulas, function(x) coxph(x, data=myFinalData)) res <- lapply(models, function(x) return(cbind(HR=exp(coef(x)), exp(confint(x))))) res
This code works, however this only gives me the hazard ratios and the confidence intervals, whereas I would like the p-values as well. I have tried multiple different things to include the p-value, but to no avail.
Additionally; if I wanted to add other factors (for ex., factor(age)), how would I modify this code in order to do that? I am not able to make it a multivariate regression model whenever I try.
Any help would be extremely appreciated!