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

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

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  • Thanks for sharing!!! May 17 '20 at 9:09 • I have a doubt. When I run coefs$df I get NULL. I guess it should be coefs[1,3] instead of coefs$df to extract the t- value from coefs, right? Aug 10 '20 at 10:35 • This is odd. You can run coefs to see the columns names. In my case the third column is called "df" so, coefs$df is equivalent to coefs[1,3]. Which command you use doesn't really matter as long as you access the correct column including the dfs.
– Juju
Aug 10 '20 at 14:48
• You need lmerTest packages. Typical lmer packages does not output degrees of freedom. After library(lmerTest), use coef(summary(model)) again. Oct 16 '20 at 2:05