Package support for `mmrm` implementing mixed models for repeated measures

[danielinteractive]

Hi @florianhartig ,

just got pointed to this nice package by @hugesingleton - we built mmrm (https://openpharma.github.io/mmrm) implementing mixed models for repeated measures and it seems there is some interest in adding DHARMa support for it. Do you have instructions for package support contributions that we could follow? Or would you have time and interest to take a stab at this yourself?

thanks, best regards Daniel

[florianhartig]

Hello @danielinteractive,

the steps to add a model to DHARMa are described in https://github.com/florianhartig/DHARMa/wiki/Adding-new-R-packages-to-DHARMA

I had a quick look at mmrm, and it seems to me that most things are unproblematic, but the functionality missing in mmrm to work with DHARMa is an option to simulate from the fitted model.

If you implement this, it would be highly preferable to also implement the re.form syntax of lme4 to condition on fitted REs when using predict() and simulate().

If these functions are be available, I could write the wrappers to include mmrm in DHARMa. Of Course, you could do that as well via a PR, but this is probably easier for me.

[danielinteractive]

Thanks @florianhartig for the fast response!

For clarification maybe worth highlighting that we don't have any random effects in the mmrm package. It is only using a parametrized covariance structure with fixed effects to model the repeated measurements. I guess in that case the re.form syntax you mention is not applicable?

[florianhartig]

Hi Daniel,

OK, on re-reading, you do state this quite clearly in the intro ;)

Excuse my ignorance though, but if you specify in https://openpharma.github.io/mmrm/main/ the term us(AVISIT | USUBJID), what does that really mean in terms of the likelihood that is estimated? Isn't this "unstructured covariance matrix for the timepoints (also called visits in the clinical trial context, here given by AVISIT) within the subjects (here USUBJID)" essentially a RE (1|subject), even it it not estimated?

It is clear that if you don't estimate the REs, you can't condition on them, so all my comments about re.form are void.

One additional comment / thought from me: while it still seems easy to technically couple mmrm to DHARMa, I am wondering now a bit about the sense of this, because I understand now that the specification of residual covariance structures is quite central to mmrm, but in this case, quantile residuals are not always easy to interpret, and usually conditioning on the covariance structure is useful.

To explain what I mean: say you fit an AR1 autoregressive model. Marginal residuals (quantile or not) will still be autocorrelated, even if you simulate. However, if you can condition on the estimated AR1 process (as you can in regression packages such as glmmTMB that fit the AR1 via REs), you can test if the remaining error is iid normal. When you can't condition, however, the residuals will carry less information.

I know to little about mmrm to understand what kind of problems your users would like to test for. I just want to warn that often, just simulating from the global model will be less informative than when you have the option to condition on your variance (RE) and covariance structures.

[danielinteractive]

Thanks Florian,

good questions :-) Actually I don't know because I did not come across the use case for DHARMa myself yet. Maybe @hugesingleton can reply your use case question?

Re: mmrm likelihood, please see https://openpharma.github.io/mmrm/main/articles/algorithm.html for the definition. Basically it is just a general linear model.

cheers Daniel

[florianhartig]

OK, I see, your likelihood is essentially a GLS and you can estimate with MLE or REML. Then

1. Indeed, there is no estimates to condition on
2. You can still simulate from the model, and that can be useful in some situations

However, you have to realise that if you have a structured error, and you can't condition on your estimates, this structured error will persist in the DHARMa residuals, so you cannot check, e.g., if you have a remaining structured error AFTER accounting for your specified covariance structure, which is something many users want to do.

I would still say that there is no harm coupling mmrm to DHARMa. Essentially, what is needed is a function that allows the user to simulate from your MLE under the model specified in https://openpharma.github.io/mmrm/main/articles/algorithm.html f

Note my comments on best practices for the signature of this function.

[danielinteractive]

@florianhartig would it be possible to standardize the residuals by multiplying with the inverse Cholesky root of the estimated covariance matrix? could that be a way to check for any "additional" structure in residuals?

[florianhartig]

Hello Daniel,

good point - Andrea-Havron contacted me about this a while ago, and I have implemented an experimental option for this in https://rdrr.io/cran/DHARMa/man/simulateResiduals.html. From the help:

Residual auto-correlation: a common problem is residual autocorrelation. Spatial, temporal and phylogenetic autocorrelation can be tested with testSpatialAutocorrelation and testTemporalAutocorrelation. If simulations are unconditional, residual correlations will be maintained, even if the autocorrelation is addressed by an appropriate CAR structure. This may be a problem, because autocorrelation may create apparently systematic patterns in plots or tests such as testUniformity. To reduce this problem, either simulate conditional on fitted correlated REs, or you could try to rotate residuals via the rotation parameter (the latter will likely only work in approximately linear models). See getQuantile for details on the rotation.

In the current implementation, you can either provide an estimated covariance structure (what you suggest) or the covariance can be estimated from simulations.

In terms of precision, the latter is clearly suboptional (because you need a lot of simulations to estimate the covariance), but I favour it because it keeps the DHARMa idea of simulation-based residuals.

Calculating a "rotation-corrected" Pearson residual is something that I would rather implement in mmrm, because this requires package-specific calculations / access to the estimated parameters.

[florianhartig]

This goes all back to the question what the users really want - I suppose what most users want is actually the DHARMa plots and not so much the specific way DHARMa calculates the residuals.

For a pure GLS, I don't really see an advantage of simulation-based residuals (the principle of DHARMa), because all necessary standardisations can be done analytically.

The only thing in DHARMa would be useful (because it requires simulations) is the simulated likelihood ratio test, which could be used to get p-values for the addition of a variance structure.

[florianhartig]

Hello @danielinteractive - I just wanted to check on the status of this - will you implement a simulate option? In this case, I'll add the option to couple to DHARMa

[danielinteractive]

Thanks @florianhartig yes we are on the way to add the required methods. Including simulate. But I guess it will take a few weeks

[florianhartig]

Awesome, and sure, no worries about the time, just tracking down open issues!

[danielinteractive]

Hi @florianhartig , just coming back to this - happy to inform that mmrm on GitHub (and hopefully soon CRAN) contains simulate and predict method (see https://github.com/openpharma/mmrm/releases/tag/v0.3.5 for all news) and I would guess that should be sufficient to get the DHARMa interface working? Please let me know if anything is missing on the mmrm side.

[florianhartig]

Hi Daniel,

sorry for replying so late - that's great news. I will check this as soon as possible!

[danielinteractive]

Thanks @florianhartig - let me know if you hit any problems with the current CRAN release of mmrm :)

[danielinteractive]

We have the new version on CRAN meanwhile :-) so the way should be free now

[florianhartig]

Hi Daniel,

thanks for the info - will put this on the list for the next DHARMa release.

[danielinteractive]

Super, thanks @florianhartig and let me know if any questions come up at any point 😄