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?


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


2 Answers 2


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:


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)
  • $\begingroup$ Thanks for sharing!!! $\endgroup$
    – Dekike
    May 17, 2020 at 9:09
  • $\begingroup$ 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? $\endgroup$
    – Dekike
    Aug 10, 2020 at 10:35
  • $\begingroup$ 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. $\endgroup$
    – Julian
    Aug 10, 2020 at 14:48
  • 1
    $\begingroup$ You need lmerTest packages. Typical lmer packages does not output degrees of freedom. After library(lmerTest), use coef(summary(model)) again. $\endgroup$
    – ESKim
    Oct 16, 2020 at 2:05

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.

  • $\begingroup$ That's awesome, thanks for sharing this information! $\endgroup$
    – Julian
    Apr 14, 2022 at 7:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.