Correct interpretation of Cox proportional hazard model with single / multiple variates I have some data on survival from a particular cancer, and also data on gender and gene expression of a single gene of interest. We have good reason to think gene expression matters, and less reason to think gender matters but we still want to analyse whether it does. My understanding is that it is best to run a multivariate analysis simultaneously rather than individually. If I run a CoxPh in R for both, I would write something like
coxall <- coxph(Surv(time, event) ~ gend*exps, data = cancer)

This yields output suggesting gender is highly significant, as shown on the table and KM plot below (shown together for clarity)

Yet I'm not entirely convinced by this, as this seems unlike other cancers I've looked at. So I ran just a test for expression and not gender, and found it highly significant (P = 6.75e-07), log hazard plot shown below:

My suspicion however is that gene expression itself is markedly different between genders. Out of curiosity, I produced the following violin plot, which shows markedly different gene expression profiles for men and women...

My question is how should I interpret this? In other cancers I've looked at, expression level matters, and gender is not significant, though in these cancers the expression levels for men and women and broadly the same and there's not clear separation like here. My hunch is that in this case gender is a proxy for gene expression, and its showing up in the initial CoxPH model for that reason. But is there anyway to test this, and how should I interpret these results?

With a high correlation between gender and gene expression, it will be hard to distinguish an association of the gene with outcome beyond what is accounted for by gender. More flexible modeling of the form of the association between outcome and gene expression might help.
Say, for example, that some unknown characteristic of females gives them better survival after diagnosis of this type of cancer. Then any gene differentially expressed between males and females will appear to be associated with outcome, even if the gene itself has no direct association with outcome. That could be what's going on in your second plot.
The joint modeling of gender and gene expression is a good way to go. With this size of a data set you might do a bit better by modeling the gene expression with something beyond a linear function of its (log) expression. For example, fitting the form of its association with outcome more flexibly with a regression spline might help to unearth an association with outcome beyond what's expected based on gender alone.
It's also good practice to include in the model other clinical characteristics expected to be associated with outcome, like age and disease stage and therapy received. In survival analysis, omitting any outcome-associated predictor from a model can bias the coefficient estimates for the included predictors.
