Problem in `predict.stan_nma_surv`.

[mattknowlesOHG]

Hi @dmphillippo,

When running through the newly diagnosed multiple myeloma example, I get the following

predict(ndmm_fit, type = "survival", level = "aggregate")
#> Error in sys.call(-3) : not that many frames on the stack

Similarly,

plot(predict(ndmm_fit, type = "survival", level = "aggregate"))
#> Error in (function (cond)  : error in evaluating the argument 'x' in selecting a method for function 'plot': i In argument: `.dup_id = 
#> 1:dplyr::n()`.
#> Caused by error: ! `.dup_id` must be size 0 or 1, not 2."

[dmphillippo]

Hi @mattknowlesOHG, thanks for raising this!

I have found the cause of the first issue; this is a bug introduced with the new marginal_effects() features in version 0.7.0. I have a fix https://github.com/dmphillippo/multinma/commit/98b1317f6ed728619cb72ad75211613af0934473. Until this is available (hopefully next week) a workaround is to specify times_seq when you call predict() for survival/hazard/cumulative hazard predictions. The times_seq argument is ignored if you specify times, but this will get you past the faulty code.

As for the second, I haven't been able to reproduce this yet... Can you let me know which version of dplyr you are using, packageVersion("dplyr")? Are you literally just running through the vignette code? If not, can you give me a small reproducible example?

[mattknowlesOHG]

Hi @dmphillippo, thanks for the quick response! I will give your workaround a try. I am using dplyr version 1.1.4, and was running through the vignette code. For what it's worth, I'm using R v4.3.1.

[dmphillippo]

Hi @mattknowlesOHG, the first issue should be fixed with the latest update (v0.7.1) which is now on CRAN. Mac binaries are available now; windows binaries are still in the queue to be built so in the meantime I suggest installing from r-universe:

install.packages("multinma", repos = c("https://dmphillippo.r-universe.dev", getOption("repos")))

As for the second issue - I can see where the error must be being produced, and this would happen if the prediction data were somehow empty for certain arms. But I can't figure out how or why this is occurring for you! I've made the predict code more robust to empty prediction data, so maybe that will fix it.

[mattknowlesOHG]

Thanks @dmphillippo, the good news is that the plots are now generating using this version! The bad news is that I still can't quite reproduce plots in the Vignette. I get this plot when running plot(predict(ndmm_fit_nph, type = "survival")).

What's interesting is if I run just predict(ndmm_fit, type = "survival"), I can't find any NAs, so I'm not sure if the issue is coming from the plotting function as opposed to the print function? I'm working on the non-PH object if that's useful to know, so I have aux_by = c(.study, .trt) in my nma() call. If I don't include that, then the prediction plot works as expected.

[dmphillippo]

Excellent, I'm glad that's working now!

Your plot output is correct. The non-PH model with baseline hazards stratified by treatment arm as well as study cannot produce predictions for the unobserved treatments in each population. The shapes of the baseline hazards in each arm are treated as separate and unrelated.

If you want a non-PH model that can produce estimates for unobserved treatment arms (i.e. can actually be used for decision-making, rather than just checking the PH assumption!) then you could fit a model that puts treatment effects (or other covariate effects) on the shape of the baseline hazard. For the M-spline model these are treatment effects on the spline coefficients. You can do this with aux_regression = ~.trt.

We talk about this in a little more detail in the final paragraph of the Assessing the proportional hazards assumption section in the vignette.

I think the .dup_id must be size 0 or 1, not 2 error was probably coming from calling predict() on the non-PH model too? The code previously didn't account for arms where no predictions could be made. The updated code handles this properly now.

[mattknowlesOHG]

Ahh this is really helpful to know, thank you so much for your help with this!

[dmphillippo]

My pleasure, thank you for reporting the bugs!

[dkarletsos]

Hi @dmphillippo,

Above you say that:

If you want a non-PH model that can produce estimates for unobserved treatment arms (i.e. can actually be used for decision-making, rather than just checking the PH assumption!) then you could fit a model that puts treatment effects (or other covariate effects) on the shape of the baseline hazard. For the M-spline model these are treatment effects on the spline coefficients. You can do this with aux_regression = ~.trt.

Would aux_regression = ~.trt stand also when a parametric distribution is used (e.g. Weibull) instead of M-spline?

[dmphillippo]

Hi @dkarletsos yes, the aux_regression argument works for any survival distribution in the package. For the Weibull for example this will put a model on the log shape parameters.