Baseline risk meta-regression

[ndunnewind]

Hi @dmphillippo. As discussed some time ago, this is an initial implementation of baseline risk meta-regression. It's not ready to be merged yet, but I first wanted to ask if you agree with this implementation.

[dmphillippo]

Thanks @ndunnewind, this is looking really nice so far.

A few thoughts and ideas below. All of this is very much up for discussion! :)

1. You have gone down the route of having a baseline_risk = TRUE/FALSE flag, which works well and is simple. However, this does currently restrict the model to having a single regression term for all non-reference treatments. This may be the most common model, but it would be nice to allow more flexible specification, e.g. independent interactions by treatment, or by class.

   a. I had thought about possibly implementing this through the regression formula interface, with another "special" for baseline risk, e.g. .mu, so you would specify the baseline risk regression as regression = ~.mu:.trt. The usual class_interactions options would apply to this. But this may be far too complicated... I do think that having the baseline_risk flag is probably more straightforward.

   b. Instead, we could have another argument to specify the type of interactions with baseline risk, similar to the class_interactions argument. E.g. baseline_interactions = c("independent", "common"). Perhaps we should refactor/alias class_interactions to trt_interactions, to better align with this new argument. Or even change baseline_risk from a boolean to a string to specify the type of baseline risk model, e.g. baseline_risk = c("none", "treatment", "class").

2. The new prior spec prior_br is nice; I would probably have shared this with prior_reg, but having separate prior_br is a little more flexible.

3. Currently this only works for agd_arm only models. Implementing for agd_contrast data should be similar. Can we get this working for models with ipd too? Is that even a model that would ever be fitted? Perhaps...

4. calculate_baseline_risk() is great, this should maybe be an internal function rather than exported?

5. Really great to put the certolizumab data in. Docs will need fleshing out a little, but we can use this for the vignette.

6. Tests in the vignette are excellent, thanks. These will need hardening against stochasticity - they fail at random at the moment due to sampling error. The solution is to run a larger number of iterations when testing.

7. Plot functionality I imagine would be along the lines of the vaccine example, i.e. demonstrate the ggplot code to create the plots, rather than add a new function? Unless this new function could also be used to plot regression terms in general?

8. The relative_effects() and predict() functions will need updating.

Thank you!

[ndunnewind]

Thanks for the extensive feedback @dmphillippo!

I removed the baseline_risk and prior_br arguments from nma(). Baseline risk meta-regression can now be specified using regression = ~.mu:.trt, like you suggested. The current implementation probably slows down all "normal" models a bit, because I moved some Stan code from transformed data to transformed parameters to make this implementation work. I'll look into this later to make sure performance of other models is not affected.

For relative_effects(), I think the sampled study-specific intercepts should be used by default to calculate study-specific relative effect estimates. A data.frame with specific values of .mu should be passed as newdata to calculate non-study-specific relative effects (see for example the vignette). Do you agree?

The sampled study-specific intercepts should also be used by default for predict() to make study-specific absolute predictions. When passing a baseline distribution, the sampled values from this distribution should be used as .mu, and .mu should not be in newdata. What do you think?

[dmphillippo]

Thanks @ndunnewind, this is looking great 🙂

The current implementation probably slows down all "normal" models a bit, because I moved some Stan code from transformed data to transformed parameters to make this implementation work.

I see you've moved the design matrices to transformed parameters; I am slightly concerned about this slowing down models. Probably there is minimal difference for standard NMAs, but IPD NMA and ML-NMR with all the numerical integration may take a heavier hit. Is there a way to do this keeping the design matrices in transformed data? Perhaps you can split off the baseline risk columns of X in transformed data and then deal with just these bits in transformed parameters and/or model blocks? You were previously adding on the baseline risk terms to the linear predictor in the model block, I think, which could work?

For relative_effects(), I think the sampled study-specific intercepts should be used by default to calculate study-specific relative effect estimates. A data.frame with specific values of .mu should be passed as newdata to calculate non-study-specific relative effects (see for example the vignette). Do you agree?

The sampled study-specific intercepts should also be used by default for predict() to make study-specific absolute predictions. When passing a baseline distribution, the sampled values from this distribution should be used as .mu, and .mu should not be in newdata. What do you think?

Yes on both of these counts, I think that's the right approach.

[ndunnewind]

I modified relative_effects() and predict() as discussed, came up with a more efficient stan implementation, and improved the vignette. I believe this is now almost ready to be merged. @dmphillippo, what do you think?

I still have some issues making the tests pass on all operating systems. Even with 10K iterations there are differences of 2 in the RE dic() pd across systems.