robustlmm package and emtrends()

[Ann1Ko]

Hi, Is there any chance to implement rlmerMod objects (from robustlmm package) as supported models to the emmeans functions? I calculated a model with the rlmer function and try to do a post hoc slope analysis using emtrends. I actually did not find any alternative solution to do these slope analysis with my robust model.

Best regards, Annika

[rvlenth]

It's probably possible. I wonder if the developer @kollerma of robustlmm is interested in adding such support. What it requires is a recover_data method (usually pretty trivial) and an emm_basis method (which is a lot like the inside of a predict method for new data). There is a vignette that provides details and examples.

I like to encourage developers to create these, as they know the nuances of their code and their model objects, and what options might be of interest. I'd be happy to help the developer get it working. Also, if changes are made or features added in robustlmm, these methods can also be updated accordingly.

[rvlenth]

In the meantime, the qdrg() function may be a way to get interim support. Since you want slopes of trends, you will have to provide two different values (most easily 1 unit apart) of the trend variable and manually compute difference quotients using update(pairs(..., simple = "x", reverse = TRUE), by = NULL) where "x" is the trend variable; this will compute the "x" difference for each combination of the other factors.

[Ann1Ko]

Ok! Thanks a lot.

As I'm relatively new to R Code, could you shortly elaborate how you would use the update/pairs function for different levels of the trend-variable?

Best regards Annika

[rvlenth]

You can't use the 'object' argument here. You need to provide all the other required documented arguments instead.

[rvlenth]

Here's an example that I got to work, based on the generated dataset ubds in my package:

> library(robustlmm)
> library(emmeans)
> rmod = rlmer(y ~ A * x + (1|B), data = ubds)
> summary(rmod)
Robust linear mixed model fit by DAStau 
Formula: y ~ A * x + (1 | B) 
   Data: ubds 

... (some output omitted) ...

Fixed effects:
            Estimate Std. Error t value
(Intercept)  15.1508     8.0242   1.888
A2           11.0407    11.2985   0.977
A3          -14.8074    13.3265  -1.111
x             0.1602     0.6708   0.239
A2:x         -0.8688     0.9419  -0.922
A3:x          1.3916     1.1041   1.260

... (rest of output omitted) ...

To use qdrg, we need to manually provide what it needs, because the rlmer results do not follow what is expected if we use object. Also, the covariance matrix from the object is of a fancy type and needs to be converted to an ordinary matrix. I'm going to specify two x values so we can compute trends.

> rrg = qdrg(formula = ~ A * x, data = ubds, coef = fixef(rmod), vcov = as.matrix(vcov(rmod)), 
+            at = list(x = c(1,2)))
> # visualize the fit:
> emmip(rrg, A ~ x)

Now let's estimate the slopes by taking differences for x:

> trends = update(pairs(rrg, reverse = TRUE, simple = "x"), by = NULL)
> trends
 contrast A estimate    SE  df z.ratio p.value
 2 - 1    1    0.160 0.671 Inf  0.239  0.8112 
 2 - 1    2   -0.709 0.757 Inf -0.937  0.3490 
 2 - 1    3    1.552 0.939 Inf  1.653  0.0984 

The slopes are shown as contrast between x = 2 and x = 1. Two things to observe:

• The slopes are consistent with what we see in the plot -- slightly positive at A=1, negative for A=2, positive at A=3.
• The statistics for A=1 are identical to those from summary(rmod) for the coefficient of x. This is correct because A=1 is the reference level for A. So this serves as a sanity check on what we did.

I hope this helps.

[kollerma]

I have just submitted a new version of robustlmm to CRAN that includes basic support for emmeans.

[rvlenth]

Terrific! I'll check it out when it's available

[Ann1Ko]

Awesome! Thanks a lot

[rvlenth]

I tried it with a model for nlme::Oats (with defaults for all but the formula and data) and it works without a hitch. So it looks luike all the needed pieces are in place. Yay!

Compared with the lmer fit of the same model, the EMMs look similar but not identical, as I'd expect. The SEs of the rlmer fit are somewhat smaller than those for the lmer fit. I was a little bit surprised by this, as I'd have expected them to possibly be a bit larger due to down-weighting some observations and such. But it all depends on the scale estimates used, etc. so am not sure if that indicates anything.

I can imagine that at some future time, you may provide more options in the vcov() method; and if so, you'll probably want to make those available to emmeans, via extra argument(s) to emm_basis.rlmerMod().

[rvlenth]

Closing this issue, as new methods in robustlmm package seem to work. I'll add mention to the models vignette