Obtaining p-values in a robustlmm mixed model via Satterthwaite-approximated DFs of the equivalent lme4 model?

Question

I've used lme4 to fit a mixed model and could obtain p-values by using the lmerTest or afex packages. However due to heteroskedasticity (Levene Test) I also fit a robust model (rlmer command in the robustlmm package). Unfortunately lmerTest and afex do not work with rlmerMod objects. So I was looking for other ways to obtain p-values.

Doing some online-research I found the recent publication by Geniole et al. (2019) and they used the Satterthwaite-approximated degrees of freedom generated by the lme4 model in combination with the output of the robustlmm model to obtain p-values for the robust model (see their Supplemental Material).

As the the summary of the rlmer model gives me a t-value I could calculate the p-values for the robust model with the help of the Dfs from the non-robust model via the following formula:

p.value.robust = 2*pt(abs(t.value.robust), Satterthwaithe.DF.non.robust, lower=FALSE)

Is that a correct way to do this? Do you know other ways to calculate p-values for the robust model?

Literature:

Geniole Shawn N., Proietti Valentina, Bird Brian M., Ortiz Triana L., Bonin Pierre L., Goldfarb Bernard, Watson Neil V. and Carré Justin M. 2019. Testosterone reduces the threat premium in competitive resource division. 286. Proceedings of the Royal Society B. http://doi.org/10.1098/rspb.2019.0720

Answer 1

I'm going to answer that myself in case that others may be interested at some point. In the end I did follow the approach mentioned above. The code would look like this:

library(lme4)
library(robustlmm)

data(Dyestuff, package = "lme4")

# fit a mixed model
model <- lmer(Yield ~ 1|Batch, data=Dyestuff)
# fit the robust equivalent
robust.model <- rlmer(Yield ~ 1|Batch, data=Dyestuff)

# get coefficients from non-robust model to extract Satterthwaite approximated DFs
coefs <- data.frame(coef(summary(model)))

# get coefficients from robust model to extract t-values
coefs.robust <- coef(summary(robust.model))

# calculate p-values based on robust t-values and non-robust approx. DFs
p.values <- 2*pt(abs(coefs.robust[,3]), coefs$df, lower=FALSE)
p.values

Answer 2

I would like to add for future reference, that the tab_model function of sjPlot package can also calculate p values based on Satterthwaite-approximated dfs.