Measuring the association of a variable with 2 different survival outcomes

I have survival data for a cohort of patients and the outcome is either cardiovascular or non-CV death. I wish to assess whether some variables make it so that the patient is more likely do die from one cause or the other.

This question is much like my own (Comparing Hazard Rate Ratios for different outcomes on the same sample), but 11 years ago with few answers.

From what I cold gather, it seems the most appropriate method for my question would be the Lunn-McNeil one. Though I have seem some published papers using it, I have found surprisingly little information about it (be it on this website or otherwise).

Much of the competing risks literature regards CIC and Fine-Gray methods, which does not seem appropriate for this particular question.

So my questions are:

• Is the Lunn-McNeil the correct method for this situation?
• If so, is it available from statistical packages (R, preferably)? Or is it usually done "by hand", so to speak?

• @EduardoMoreira you can simply run tests on functions combining the $\beta_iX$ and $\beta_iY$ estimates for predictor $i$ and outcomes $X$ and $Y$. As the $\beta$ coefficients (in log-hazard scale) are (asymptotically) multivariate normal, after you have done the modeling you can use the coefficient point estimates and their variance-covariance matrix to do such comparisons. The simplest would be to evaluate $\beta_iX - \beta_iY$ to see if they are significantly different. Tests on the ratio $\beta_iX / \beta_iY$ are likely to be problematic, particularly if $\beta_iY$ is close to 0.