[DRAFT] adding support for the {mmrm} package

[kkmann]

Hi,

this is an attempt to implement support for the {mmrm} package https://github.com/openpharma/mmrm following the steps outlined in the {marginaleffects} documentation https://marginaleffects.com/vignettes/extensions.html#marginaleffects-extension.

I am still unclear about the role of the option setting.

options("marginaleffects_model_classes" = "lm_manual")

Is it expected that the user triggers this everytime they want to use {marginaleffects} with {mmrm}? This does not seem to be the case for other extensions.

My test case is

library(mmrm)
fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)
predictions(fit, newdata = insight::get_data(fit))

Also tagging @lang-benjamin @danielinteractive

[kkmann]

Currently my test case fails irrespective of me setting

options("marginaleffects_model_classes" = "mmrm")

or not.

Error: Unable to compute predicted values with this model. You can try to supply a
  different dataset to the `newdata` argument.
  
  Bug Tracker: https://github.com/vincentarelbundock/marginaleffects/issues

I did not open a separate issue in the bug tracker since this is likely a bug caused by my faulty implementation of the {mmrm} interface.

[vincentarelbundock]

Thanks for taking the time to do this, I really appreciate it!

The options() call is only necessary for temporary user-defined models. When we merge support in marginaleffects itself, we can do the following:

1. Add the package to the permanent list of supported models, to avoid calling options() in every session: R/sanitize_model.R
2. Add the package to the vignette of supported models by editing: data-raw/supported_models.csv
3. Add a bullet point to the news file: NEWS.md

Does your example work when you just supply data frame directly: newdata=fev_data?

[lang-benjamin]

Thanks, just adding some thoughts here.

Does your example work when you just supply data frame directly: newdata=fev_data?
Yes, that should work (fev_data has missings for the response and there will then be predictions for it).

The question to me is what the general default strategy regarding missing data should be. The element model$data will not omit missing data, whereas insight:get_data(model) will omit it (i.e. it uses the actual data used to fit the model). If we are interested in the actual data that was used to fit the model, instead of insight::get_data one could just use model.frame(model). This would not rely on the insight package.

[vincentarelbundock]

Ideally, on my end, I would like to go through insight::get_data(), and use listwise deletion by default. This would be consistent with how all the other models are treated in the package. In any case, it's trivial for users to supply another dataset to the newdata argument, if they want predictions for observations with missing response.

Also, looking at the code quickly, I suspect that the reason predictions don't work in the current PR is that get_predict() just calls predict(). But that function should return a data frame with two columns: rowid and estimate. If you only return a vector, it will like return this kind of error.

It might also be worth checking if you really need all those methods. Sometimes, the defaults work perfectly well. Just try something like:

library(marginaleffects)
get_predict.default(fit)

[kkmann]

Thanks for the feeback @vincentarelbundock, made some adjustments and starting to wrap my head around this. If I see it correctly we still have two things to address here.

• we still need to add the information into the csv file @lang-benjamin
• we need to look into the omitting missing data either via insight::get_data() or custom

[vincentarelbundock]

I was looking at this PR today and noticed that the current version does not return standard errors. The underlying reason for this is that the Jacobian is full of zeros:

library(marginaleffects)
library(mmrm)
fit <- mmrm(FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
            data = fev_data,
            weights = fev_data$WEIGHTS,
            method = "Kenward-Roger")
p <- predictions(fit)
attr(p, "jacobian")

The reason for that, is that get_predict() does not appear to make different predictions after we change the coefficient values stored in the model object using set_coef(). Under the hood, the way we compute standard errors is by slightly tweaking the coefficient estimates in a model, and making different predictions with the modified model object. This doesn't seem to work at the moment.

One possibility is that this is related to the problem with glmmTMB: #1064

I don't actually know if this is related, but I noticed that both glmmTMB and mmrm use Template Model Builder (TMB). If that's the source of the issue, it might be bad news... See the thread on glmmTMB for links to different attempts and problem definition.

[danielinteractive]

Thanks @vincentarelbundock and @kkmann and all, just looking at this thread as a bystander - yes, structurally the two packages will work similarly. However I also wonder if in this case, because we use TMB, there would be much better automatic differentiation ways to generate the Jacobian you need.

[vincentarelbundock]

@danielinteractive Thanks for chiming in!

Yes, that would certainly be more efficient. I think it makes sense for very common quantities like predictions. But the issue is that marginaleffects can compute literally dozens of quantities of interest (slopes, risk differences, marginal means, etc.), and we have to build distinct jacobians for all of these quantities.

My guess is that this would require a ton of work on the backend. But then again, I don't know TMB at all, so...

[clarkliming]

PR created at kkmann#1 to fix the issue that the jacobian is 0

reason is that predict in mmrm is always returning a conditional prediction (conditional on observed values)

[vincentarelbundock]

Sorry I left this PR to linger for so long. I looked into this today, and things mostly seem to work, but there are two lingering issues. The first is probably easy to solve. Not sure about the second.

library("mmrm")

fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)

newdata=NULL needs to include the response variable by default

By default, when the user leaves newdata=NULL (the default), marginaleffects uses sanitize_newdata() to extract the original data from the model object. Typically, this does not include the response, but it appears that the predict() method supplied by mmrm requires the response to make predictions.

I’m not sure if that’s required by the method. If it’s not actually required, it would be best if this were fixed upstream in mmrm. If this is a strong requirement, then we can implement a hack in marginaleffects to include the response.

Default breaks:

avg_predictions(fit)
> Error: Unable to compute predicted values with this model. You can try to supply a
>   different dataset to the `newdata` argument. This error was also raised:
>   
>   object 'FEV1' not found
>   
>   Bug Tracker: https://github.com/vincentarelbundock/marginaleffects/issues

Supplying the dataset explicitly is a workaround:

avg_predictions(fit, newdata = fev_data)
> 
>  Estimate Std. Error   z Pr(>|z|)   S 2.5 % 97.5 %
>      42.4      0.329 129   <0.001 Inf  41.8   43.1
> 
> Type:  response 

Duplicate indices

Internally, marginaleffects computes comparisons/slopes for multiple terms or comparisons in one shot, by stacking datasets. That way, we can compute the slope with respect to X and the slope with respect to Z in one shot, which is much more computationally efficient, since it doesn’t require us calling predict() too many times. Unfortunately, this breaks with mmrm, because the predict() method supplied by that package requires unique indices, and “stacking” duplicates the indices.

This means that we can compute the “effect” of a focal variable by specifying that varible explicitly, it

avg_comparisons(fit, 
    variables = "SEX",
    newdata = fev_data)
> 
>  Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
>     0.326      0.532 0.613     0.54 0.9 -0.717   1.37
> 
> Term: SEX
> Type:  response 
> Comparison: mean(Female) - mean(Male)

It also works if we specify a comparison between only two levels for a single focal variable:

avg_comparisons(fit, 
    variables = list("RACE" = c("Asian", "White")),
    newdata = fev_data)
> 
>  Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
>      5.64      0.666 8.48   <0.001 55.3  4.34   6.95
> 
> Term: RACE
> Type:  response 
> Comparison: mean(White) - mean(Asian)

But it doesn’t work if we try to compute the two possible contrasts at the same time: (1) Asian vs White and (2) Black vs. White.

avg_comparisons(fit, 
    variables = "RACE",
    newdata = fev_data)
>   Error in h_mmrm_tmb_data(object$formula_parts, newdata, weights = rep(1, : time
>   points have to be unique for each subject, detected following duplicates in data:

This error is related to the lines of code after this one: https://github.com/openpharma/mmrm/blob/424810d88b0507e9e88242b61d343852358b2df2/R/tmb.R#L122

I don’t know mmrm to know if this is essential. What happens if we bypass that check? Will the hardcoded ordering yield incorrect estimates?

[danielinteractive]

Thanks @vincentarelbundock , looking into this now on the mmrm end - thanks for the detailed pointers!

[danielinteractive]

Thanks again @vincentarelbundock! While we will be able to fix the first problem quickly (indeed this is a bug in mmrm), the second problem cannot easily be fixed on the mmrm side I think. For non-spatial covariance structures we need unique IDs/visit variable combinations here, otherwise downstream computations are not well defined. Is there maybe a point in the marginaleffects stack to define a method for mmrm objects such that in the stacking some kind of dummy renaming of the IDs is performed, ensuring unique IDs/visits for calling mmrm afterwards?

[vincentarelbundock]

Cool. Thanks for taking a look and finding a solution for issue 1.

I just spend a bunch of time trying to find a place to insert a hack that would append some random string to the subject ID column. Unfortunately, it's very hard. The problem is that this stacking logic is very "deep" in the marginaleffects. There isn't just one stacking. We stack when computing contrasts for different variables, but we also stack when computing several contrasts for a single variable (ex: different factor levels). And the stacking code is separate for different data types: numeric, logical, factor, etc.

So in order to support this, I would have to insert complex hacks in a dozen places in the code base, and don't think that's sensible design for software that wants to support 100+ model classes.

To be sure, that's less convenient, but maybe it's not the end of the world that users have to define a single contrast they are interested in, rather than get all of the them at once. This can probably be documented.

[danielinteractive]

Hi @vincentarelbundock , I discussed the issue openpharma/mmrm#461 further with @clarkliming , and it seems that a more principled solution to that will likely also fix the issue 2 here. So let us first sort that out and maybe both issues 1 and 2 will be solved in mmrm.

[vincentarelbundock]

Oh wow, OK. I really appreciate your help in working on this. It'll be really cool.if we can push this through the fi ish line

[adrianolszewski]

Hello @vincentarelbundock , @danielinteractive,

I just checked the mmrm repository and it seems that the main issue has been resolved:
openpharma/mmrm#463
openpharma/mmrm#461

Does it mean that the integration with marginaleffects is a step closer?
Such achievement would be an great news for the clinical trials industry interested in using R!

[danielinteractive]

Thanks @adrianolszewski for circling back on this, I am not sure if we solved the second issue, @clarkliming do you remember? Maybe @vincentarelbundock worth to try the code again with the current {mmrm} version?

[vincentarelbundock]

I merged the main branch into this feature Pull Request and tested a simple mmrm model. Things seem to be working OK.

I need your help! If you have experience with mmrm models or mixed models for repeated measures, please test this code and let me know if the results make substantive sense. While the functions appear to be working technically, I lack domain expertise to validate whether the slopes, predictions, and comparisons are mathematically and substantively correct for these model types.

Installation

First, install the development version with mmrm support:

remotes::install_github("kkmann/marginaleffects@mmrm-extension")

Load packages and fit model

library(marginaleffects)
library(mmrm)

# Fit an mmrm model using the built-in fev_data
# This models FEV1 (lung function) with an unstructured covariance matrix
fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)

Test core marginaleffects functions

Predictions

Get unit-level predictions and average predictions:

# Unit-level predictions (first few rows)
predictions(fit) |>
  head() |>
  data.frame()
#>   rowid estimate std.error statistic       p.value  s.value conf.low conf.high  df     FEV1                      RACE    SEX ARMCD AVISIT USUBJID
#> 1     1 41.20593 0.7662678  53.77484  0.000000e+00      Inf 39.70407  42.70779 Inf 39.97105 Black or African American Female   TRT   VIS2     PT1
#> 2     2 52.08639 1.2751130  40.84845  0.000000e+00      Inf 49.58722  54.58557 Inf 20.48379 Black or African American Female   TRT   VIS4     PT1
#> 3     3 35.61706 0.7636975  46.63766  0.000000e+00      Inf 34.12024  37.11388 Inf 31.45522                     Asian   Male   PBO   VIS2     PT2
#> 4     4 41.11959 0.6454274  63.70909  0.000000e+00      Inf 39.85457  42.38460 Inf 36.87889                     Asian   Male   PBO   VIS3     PT2
#> 5     5 45.83137 1.2760045  35.91788 1.606758e-282 936.0996 43.33045  48.33230 Inf 48.80809                     Asian   Male   PBO   VIS4     PT2
#> 6     6 37.47363 0.7408669  50.58078  0.000000e+00      Inf 36.02155  38.92570 Inf 35.98699 Black or African American Female   PBO   VIS2     PT3

# Unit-level predictions for an individual with average characteristics
predictions(fit, newdata = "mean")
#> 
#>  Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %  Df
#>      35.9      0.736 48.8   <0.001 Inf  34.5   37.4 Inf
#> 
#> Type: response

# Average predictions across the dataset
avg_predictions(fit)
#> 
#>  Estimate Std. Error   z Pr(>|z|)   S 2.5 % 97.5 %  Df
#>      42.3       0.33 128   <0.001 Inf  41.7     43 Inf
#> 
#> Type: response

Comparisons

Calculate average marginal comparisons for the SEX variable:

avg_comparisons(fit,
  variables = "SEX",
  newdata = fev_data
)
#> 
#>  Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
#>     0.326      0.532 0.613     0.54 0.9 -0.717   1.37
#> 
#> Term: SEX
#> Type: response
#> Comparison: Female - Male

Slopes (Marginal Effects)

Calculate average marginal effects:

avg_slopes(fit)
#> 
#>    Term                          Contrast Estimate Std. Error      z Pr(>|z|)     S  2.5 % 97.5 %
#>  ARMCD  TRT - PBO                            3.752      0.668  5.621   <0.001  25.6  2.444   5.06
#>  AVISIT VIS2 - VIS1                          4.819      0.565  8.534   <0.001  56.0  3.712   5.93
#>  AVISIT VIS3 - VIS1                         10.004      0.594 16.850   <0.001 209.2  8.840  11.17
#>  AVISIT VIS4 - VIS1                         15.358      0.926 16.590   <0.001 202.9 13.544  17.17
#>  RACE   Black or African American - Asian    1.530      0.624  2.451   0.0143   6.1  0.307   2.75
#>  RACE   White - Asian                        5.644      0.666  8.479   <0.001  55.3  4.339   6.95
#>  SEX    Female - Male                        0.326      0.532  0.613   0.5399   0.9 -0.717   1.37
#> 
#> Type: response

[danielinteractive]

Thanks a lot @vincentarelbundock !
@kkmann do you think you have any time to test this?

[lang-benjamin]

Thank you so much @danielinteractive @vincentarelbundock.
Just as a note so that we (at least I) don't forget about this: I think we should have another look at the default value for newdata in get_predict.mmrm as well. It may not reflect the data that was actually used in fitting the MMRM (which is the intention of the newdata argument). I can also hopefully have a closer look into that in the near future.