I am performing a multivariate Cox regression analysis, and would like to find what combination of those covariates best predict my outcome.
Say I have a list of candidate genes whose expressions showed (1) to be associated with overall survival (OS) (Cox regression), and (2) also associated among themselves (multivariate linear modeling). For example, high levels of gene_1 AND low levels of gene_2 AND low levels of gene_3 are associated with poor prognosis. So, if I am using only these 3 variables and a decent number of patients (e.g. n=200), it is not hard to find those patients with that combination of genes outcome.
However, if for instance I have a list comprising 7 of such candidates, the chances of finding a patient that fits this criterion (now for these 7 hits) is nearly none.
So my question is: is there a way to do some sort of permutation/combinatorial analysis coupled with Cox regression to find the combination of those 7 targets that best associates with OS? Considering that a satisfactory combination of factors is represented say in at least 25% of the patient population.
> shortlist[1:5,1:10] Age PFS PFS_codex OS OS_codex gene_1 gene_2 gene_3 gene_4 gene_5 Sample_1 67.9 117.23 0 115.69 0 9.451046 5.572303 7.260597 8.492154 4.010582 Sample_2 61 69.27 0 72.30 1 9.520935 9.956700 8.370941 6.854242 4.638455 Sample_3 69.1 1.23 1 1.08 1 9.691664 8.713712 8.840432 7.891189 3.707268 Sample_4 72.2 15.27 1 69.63 1 9.490668 9.015255 8.601908 9.584230 4.277126 Sample_5 40.7 61.43 1 78.41 1 9.439942 7.769697 7.337121 7.222432 4.843211
This shortlist is already with only those candidate that passed a Cox Regression univariate analysis. So now, I run Cox again, with the exception that this time all candidates #from the shortlist are put up together:
multicox <- coxph(Surv(OS, OS_codex) ~ gene_1 + gene_2 + gene_3 + ... + gene_27, data=shortlist)
The results are the following (because of space, only those with significant p-values are listed):
> summary(multicox) Call: coxph(formula = Surv(OS, OS_codex) ~ gene_1 + gene_2 + gene_3 + ... + gene_27, data = shortlist) n= 196, number of events= 133 coef exp(coef) se(coef) z Pr(>|z|) gene_1 0.80167 2.22926 0.22501 3.563 0.000367 *** gene_8 0.15332 1.16570 0.06417 2.389 0.016888 * gene_9 0.76781 2.15505 0.24603 3.121 0.001803 ** gene_12 -0.84846 0.42807 0.30868 -2.749 0.005984 ** gene_18 -0.60773 0.54459 0.13789 -4.407 1.05e-05 *** gene_26 0.35992 1.43321 0.14591 2.467 0.013634 * gene_27 0.41905 1.52052 0.14972 2.799 0.005127 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
There are 7 candidates showing association to the OS (considering p<0.05). Based on that, and on the HR (to know whether upregulation or downregulation is associated to OS) of each one of them, I would like to test the hypothesis that the combination of all these candidates are indeed good predictors for poor prognosis.
For that, I would like to artificially create 2 groups of patients: 'high risk' and 'low risk'. The 'high risk' group is the one which the up- or down-regulation for each candidate showed to be associated with OS, respectively. The 'low risk' is the rest. (OBS: here, I am considering up- or down-regulation all values that are above and below the median for each gene, respectively)
First, I store the median values for all the 7 genes:
> a <- data.frame(gene_1=median(shortlist$gene_1), gene_8=median(shortlist$gene_8), gene_9=median(shortlist$gene_9), gene_12=median(shortlist$gene_12), gene_18=median(shortlist$gene_18), gene_26=median(shortlist$gene_26), gene_27=median(shortlist$gene_27)) > a <- as.numeric(a) > a  8.999678 5.681134 5.907599 8.420542 6.158107 3.279144 7.020527
So the classification now comprises the following:
shortlist$risk_class <- ifelse(shortlist$gene_1>a & shortlist$gene_8>a & shortlist$gene_9>a & shortlist$gene_12<a & shortlist$gene_18<a & shortlist$gene_26>a & shortlist$gene_27>a, "high_risk", "low_risk")
However, when I do that, there is only one patient that fits in that criteria:
>sum(shortlist$gene_1>a&shortlist$gene_8>a&shortlist$gene_9>a&shortlist$gene_12<a&shortlist$gene_18<a&shortlist$gene_26>a&shortlist$gene_27>a)  1
So, my question is, would there be a way to test for the best combination of those 7 candidates, with a reasonable number of patients (i.e. >30)? I was thinking about some sort of permutation/combinatorial analysis coupled with Cox regression to find the combination of those targets that best associates with OS?
gene_1 + gene_9 + gene_18 + gene_26 OR
gene_1 + gene_18 + gene_27 OR
gene_9 + gene_18 + gene_26 + gene_27 OR
gene_18 + gene_27 ... and so forth.
And finally, would there be any package in R that can perform it? Thanks a lot!