add getSurvival

[0xMuluh]

ping #93

[0xMuluh]

This wrapper currently defaults to using pairwise log-ratio (PLR) transformations via coda4microbiome::coda_coxnet() , consistent with the approach described in Calle et al., 2024.

plr was selected by the authors over clr, ilr, or alr for reasons related to interpretability, sparsity, and model perfomance in high-dimensional microbiome survival analysis. While this default appears justified and performs well empirically, it may not be optimal in all cases.

If plr proves sufficient across expected use cases, we can add the unit tests around this assumption.

If user needs or new evidence suggest flexibility is needed, we should refactor to allow other transformation via a parameter (e.g., transform.method = "plr", "clr", "ilr", etc.). or mia::transformAssay()

Currently the method returns a list of objects including ggplots and dfs/matrices. We can consider addSurvival() and store in metadata?

Open to discussion on whether this needs to be prioritized now or flagged for future maintenance..

[TuomasBorman]

Thanks for this clear PR!

One thing to discuss is how much functionality we want to utilize from coda4microbiome package. If there are plans for broader interlinking, this PR might be fine. However, under the hood, coda_coxnet() is utilizing survival::coxph() and glmnet::cv.glmnet() and survminer for plotting.

The implementation in coda_coxnet() is not complex but includes some overhead that -- in my opinion -- does not necessarily fit well to our current design choices (transformation, printing, plotting).

Instead, this functionality for cox model could be implemented directly in miaTime package. We could take inspiration from the coda4microbiome package, but it seems that there is licence mismatch so this should be taken into account.

[antagomir]

Thanks!

Indeed, I would not wrap coda4microbiome. We may like to test and run standard survival models (Cox) also with other transformation.

So let's keep separate:

• data preparation (e.g. plr -> transformAssay)
• statistical analysis / testing (e.g. Cox model, and glmnet if needed for feature selection)
• summaries (e.g. tables)
• visualizations (-> miaViz?)

This can take inspiration from coda4microbiome and cite it but Cox & GLM are more general and previously published methods. Moreover, the coda4microbiome paper did not include many comparisons to alternative transformations & methods. Our implementation should be modular and support different combinations of choices (e.g. transformation -> testing -> presentation).

The miaTime already includes data from the paper as example data, that could be used to verify the workflow (should yield the same result).

[0xMuluh]

Thoughts. @antagomir @TuomasBorman

[0xMuluh]

Consider adding pairwise logratios plr to transformAssay()

[antagomir]

I think that mia PR #740 already provides this.

[TuomasBorman]

Thanks! Can you also add unit tests?

[TuomasBorman]

Check argument names that we are already using. I think we are using col.var which would suit also here

[TuomasBorman]

Could also be defined as a formula

[0xMuluh]

Thanks! Can you also add unit tests?

unit tests added after #107

[TuomasBorman]

Thanks!

[TuomasBorman]

It is not evident what this method does.

Can you add Detail section?

[TuomasBorman]

Coefficients coul dbe added to rowData and risk_scores to colData in addSurvival

[TuomasBorman]

What is the intention of multiplying the values with 2?

[0xMuluh]

to have interpretable contrasts that live in space [-1, 1]. after the operation L1 norm is 2

[antagomir]

In manpage description clearly cite the publication where the data is from (note that there is a special way to add citations in roxygen). Now you only refer to R package. But there is also formal publication to cite.

Cite all references in the manpage text so that it is clear why they are added on the reference list.

Running the manpage examples gives error

The assay contains already only strictly positive values. Pseudocount is not added.
Error in list(coefficients = coefs, risk_scores = risk_scores, c_index = c_index, :
argument 6 is empty

Either the function name should be getCox or something explicit, or if it is getSurvival then there should be a method argument that specific what survival method one wants to use.

[antagomir]

These are time-to-event models. How about renaming status.col into event.col?

[antagomir]

getSurvival has

A list with model summaries:

    • ‘coefficients’: estimated model coefficients

    • ‘selected_features’: nonzero features in penalized model

    • ‘risk_scores’: predicted risk scores

    • ‘cindex.apparent’: apparent concordance index

    • ‘cv.cindex.mean’: mean cross-validated C-index (if penalized)

    • ‘cv.cindex.sd’: SD of cross-validated C-index (if penalized)

    • ‘fit’: fitted model object (coxph or cv.glmnet)

-> Harmonize to use "_" as the separator for all (do not mix)

-> cindex.apparent could be C.apparent, or just C ?

-> For clarity I would rename cv.cindex.mean to C.mean.cv and same for the SD ting

-> should selected_features be just nonzero_features if these are the same thing

-> coefficients is usually just coef in many R methods -> harmonize?

[0xMuluh]

In manpage description clearly cite the publication where the data is from (note that there is a special way to add citations in roxygen). Now you only refer to R package. But there is also formal publication to cite.

The data has not been used. Will use this data instead: b752170

[0xMuluh]

Running the manpage examples gives error

The assay contains already only strictly positive values. Pseudocount is not added.
Error in list(coefficients = coefs, risk_scores = risk_scores, c_index = c_index, :
argument 6 is empty

Thanks, there was a typo , in the return list.