~.trt in aux_regression not consistently included in predict?

[rogula]

Dr. Phillippo - Thank you so much for making this package available. We have been using the nma function for modeling survival data, and have benefited greatly from the code and documentation you have shared.

When specifying both the ‘regression’ and ‘aux_regression’ arguments in the nma function, we have encountered what we think might be a bug when calling the predict function: the ~.trt coefficient appears to be ignored when predicting e.g., type = "survival" and other 'types'.

By adapting your R code (simulated_parametric_survival.R), we can reproduce this error using the following changes:

First, update the simulated data to produce non-proportional hazard (PH) datasets: survdat_AB <- simsurv("weibull", lambdas = weib_sim_pars["scaleAB"], gammas = weib_sim_pars["shapeAB"], x = covdat_AB, betas = betas, maxt = 1, tde = c(trtB = 5)) %>% mutate(status = if_else(runif(nAB) <= cens_rate, 0L, status))

survdat_AC <- simsurv("weibull", lambdas = weib_sim_pars["scaleAC"], gammas = weib_sim_pars["shapeAC"], x = covdat_AC, betas = betas, maxt = 1, tde = c(trtC = 1)) %>% mutate(status = if_else(runif(nAC) <= cens_rate, 0L, status))

Then, when fitting the following model (using the IPD dataset), the model coefficients correctly capture the treatment effect on the auxiliary parameters (we have confirmed this separately), but the predict function does not appear to include it - S(t) plots do not match those we generate separately and do not reflect non-PH curves: weib_IPD_nph1 <- nma(sim_net_IPD, regression = ~(x1 + x2 + x3)*.trt, likelihood = "weibull", prior_intercept = normal(0, 100), prior_reg = normal(0, 100), prior_trt = normal(0, 100), prior_aux = half_normal(10), aux_regression = ~.trt, QR = TRUE)

Interestingly, the predict function appears to work as expected when the ‘regression’ argument is left as the default, as in the model specified below:

weib_IPD_nph2 <- nma(sim_net_IPD,

regression = ~(x1 + x2 + x3)*.trt,

                 likelihood = "weibull",
                 prior_intercept = normal(0, 100),
                 prior_reg = normal(0, 100),
                 prior_trt = normal(0, 100),
                 prior_aux = half_normal(10),
                 aux_regression = ~.trt,
                 QR = TRUE)

We have explored other parameterizations which appear to work as expected using the predict function, but require additional covariates in the ‘aux_regression’ such as this one: weib_IPD_nph4 <- nma(sim_net_IPD, regression = ~(x1 + x2 + x3)*.trt, likelihood = "weibull", prior_intercept = normal(0, 100), prior_reg = normal(0, 100), prior_trt = normal(0, 100), prior_aux = half_normal(10), aux_regression = ~.trt + x3:.trtclass, QR = TRUE)

We are happy to provide additional details if it would be helpful. Please also let us know if this error is due to our misunderstanding how to correctly parameterize non-PH models.

R version 4.4.0 (2024-04-24 ucrt) Platform: x86_64-w64-mingw32/x64 Running under: Windows 11 x64 (build 22631) multinma_0.7.1.9000

[dmphillippo]

Hi @rogula, thank you very much for reporting this and for the super clear reproducible example - I really appreciate it!

I believe that I have found the cause of the problem and fixed this. As you say, the model was estimated properly but then the aux_regression treatment effects were omitted from the predictions in some cases.

The latest development version now contains the fix, you can install this from r-universe with:

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

Hopefully this resolves the problem for you? Note that in my testing the weibull non-PH model still doesn't fit very well to the B arm of the AB trial due to the very strong time-varying effect. Using an M-spline non-PH model does much better.

[rogula]

We installed the new package and this appears to be resolved. This is amazing, thank you! We have also been using spline models for our data, and agree they offer better fit to the non-PH data (although we are encountering convergence issues); we used the Weibull as an example simply because it was easier to reproduce externally in order to confirm the error.

[dmphillippo]

Excellent!

Regarding the convergence issues, when you fit a non-PH spline model with aux_regression, the model must use the same knot locations across all studies. This restriction is not in place when aux_regression is not specified. See the documentation of the knots argument in nma() and the make_knots() function.

When each study has its own knots, the knot placement is straightforward and does not really affect convergence or fit. But when we have to share knots across all studies this can lead to issues when the studies have different lengths of follow-up. The reason is that when a knot is placed just before the last event time(s) in a study, there is then not much information to estimate the spline hazard in the following interval. Placing a knot just after the last event time has no such issue.

Calculating good knot locations when they must be shared across all studies is hard to automate - the package has a good go with the make_knots() function, but the knot locations generated are still not ideal in all scenarios.

I would suggest looking at the length of follow-up in all studies, and comparing those to the knot locations - look at the code in the make_knots() examples for ideas on how to do this. Then tweak the knot locations moving/adding/removing before fitting the model with these knots, making sure that there are no knots just before the end of follow-up in any studies.

This is not really documented anywhere yet (sorry!). I am currently working on a paper that presents the M-spline model and the non-PH regression which will have all of these details in 🙂

[rogula]

Thank you! We do appear to have sufficient events in each study between knots, so we're exploring to see if we can diagnose/resolve the cause of the divergent transitions. It's great to hear that an M-spline and non-PH paper is in the works - we look forward to it!