There's no right answer to this, but it felt too long as a comment. Anyway, here are a few thoughts, hopefully supplemented by the answers of others.

Guides I have found online recommend either combining correlated features together, or dropping features entirely. Given that the VIF does not indicate with which other features a particular feature is correlated, it is difficult to automate 'merging' of correlated features across all 31,867 RNA splicing events. I have also read that dimensionality reduction techniques such as PCA (Principal Component Analysis) will produce uncorrelated output features that can be used for downstream analysis, but I am unsure of how to determine how many principal components to form, or which of the features making up each principal component are responsible for causing certain principal components to exhibit a significant relationship with survival in the Cox model.

Finally, I am not sure what the accepted 'gold standard' is when selecting 'top' features from a random forest classifier. I am currently taking the top 0.1% of features, but this was a largely arbitrary choice. Is there a particular method to determine this threshold, or a convincing argument for a commonly used threshold?

Agnostic, eigenvector-based decomposition algorithms You probably know a lot of this, but I am writing for a more general audience so please bear with me. Principal Component Analysis (PCA) is a form of eigenvector-based decomposition. It is a descriptive statistical technique; it does not ask questions about what variables might be collinear with what; rather, it simply takes the data, finds the longest axis of variation, then records scores for the degree to which each subject "belongs to" that principal component. Then, the algorithm successively finds the next longest axis of variation, with the stipulation that it must be {orthogonal to, independent of, not correlated to, not related to} the 1st one. Then, it goes on and on and on down the way successively producing the next largest axis of variation and the next.

There is really no number of PCs that is fundamentally meaningful a priori. For other eigenvector-based decomposition algorithms (there are many), you might have to specify a number of eigenvectors ahead of time, so the question is more meaningful in those cases. But for PCA, you can continue going till you've decomposed ALL the variation in the data into orthogonal axes. If that's what you want to do, that algo can absolutely do it for you.

But I don't recommend that, although ironically PCA could be used for this project in the following way:

Quality Control

- Run a PCA on your data.
- Correlate the top 10, or 20, or 50, 100, etc. (doesnt really matter) to all of your covariates (be sure to choose the right kind of correlation coefficient). The goal of this is to identify if any of the PCs correlate with any labels (sex, age, BATCH (i.e., can scan for batch effect), treatment status, etc. 
- The purpose of this is to identify if an **UNWANTED** variables creating collinearity in your data, that do not correspond to what you want to measure (i.e., over and above what should be present biologically). Consider reading Price *et al.* (2006) to get a sense for how this was done for GWAS studies.

Once you are satisfied about QC, batch effect, PCA, and your data, you can move on to biological analysis of a dataset that is, hopefully, more free of unwanted covariance structures.

Biological Investigation: The genome is an orchestra, not a guitar solo:

Given that the aim of this research project necessitates model interpretation in order to determine the RNA splicing events/isoforms that have a significant relationship with survival, my understanding is that the multicollinearity within my dataset must be removed.

In my opinion, this thought is incredibly damaging to an understanding of biology. Almost no single datum is important on its own. Take BRCA1 positivity. I have the most famous breast cancer gene, so I will get breast cancer. Right? Wrong. Penetrance varies from 8-80% depending on the rest of your genetic background - in other words depending upon covariance.

Rather, it is important that the above QC process is done well, because biological data themselves are and ought to be heavily collinear. Why? Because genes don't function on their own, they organize into higher and higher and higher order structures to accomplish much larger scale objectives. As a simple example, think of how non-sensical it would be for a bacteria to upregulate one gene that makes up a flagellar subunit, without any others. What would it do with just one piece? It couldn't make a flagellum, but the protein is only useful insofar as it is part of a flagellum...so would be pointless.

So, the trick here is to be able to relate many changes made in a coordinate fashion to some other known biological entity... For instance, let's say you have a bunch of splice variants, but you can map them all back to interferon gamma signaling. You might also find that TLR3 signaling and Stat1 signaling also relate; for that dataset all those splicing signatures would be collinear and that would not be in any way problematic, rather, it is evidence that the data are meaningful.

In other words, I am suggesting that you derive meaning from your data not by resolving it into uncorrelated axes, but rather by trying to understand what the relationships that drive large-scale differences in splicing are.

I don't know much about analyzing splicing data (never been asked to do it) but if it were me, I would probably get a list of all gene products that actually carry out splicing, then get lists of every gene those splicing proteins are thought to act on. I would then see what signatures match my data. If you find a strong match, this could also give you a starting point inasmuch as it would tell you what gene complexes are doing the splicing. That plus cell type plus literature review might be enough to begin to speculate on what is happening.

Spectral Decomp and a recent manuscript Consider reading this preprint. The authors use spectral decomposition on the proteosome of several thousand bacteria. They identify layers upon layers of nested, hierarchical structures of covariance by analyzing spectral components of different magnitudes.

Notice that the top few SCs recapitulate kingdom, phylum, class, order, family, genus, species (1-40), then it goes subcellular. By the time you get to SCs in the 1000s, you are looking at minute covariance structures, perhaps such as a few dozen genes mediating some aspect of OxPhos. This manuscript may give you some ideas about how to use eigenvector based decomp for your own data.

Hope that helps in some measure.