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.

#An an example:

> 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
[1] 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[1] & shortlist$gene_8>a[2] & shortlist$gene_9>a[3] & shortlist$gene_12<a[4] & shortlist$gene_18<a[5] & shortlist$gene_26>a[6] & shortlist$gene_27>a[7], "high_risk", "low_risk")

However, when I do that, there is only one patient that fits in that criteria:

>sum(shortlist$gene_1>a[1]&shortlist$gene_8>a[2]&shortlist$gene_9>a[3]&shortlist$gene_12<a[4]&shortlist$gene_18<a[5]&shortlist$gene_26>a[6]&shortlist$gene_27>a[7])
    [1] 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?

For example:

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!