How can I get confidence intervals for fixed effects using the rlmer function (robustlmm package)?

Question

I need to fit a linear mixed model but my dependent variable is right-skewed with some big outliers. Thus I used the rlmer function of the robustlmm package. it works quite nicely, however what I am missing now is confidence intervals for the fixed effects. Is anyone working with this package and has some tips for me? Otherwise I would also be interested in other packages to fit robust linear mixed models.

Answer 1

• Wald confidence intervals: these assume that the sampling distribution of the parameters is multivariate Normal (a much weaker assumption than that the conditional distribution of the residuals is Normal). They are relatively easily to compute (for the fixed-effects parameters) by extracting the parameter values (fixef()) and the standard errors (sqrt(diag(vcov()))) and computing $\beta \pm z \cdot \sigma$ ...

Example:

library(robustlmm)
r <- rlmer(Reaction ~ Days + (1|Subject), sleepstudy)

confint.rlmerMod <- function(object,parm,level=0.95) {
     beta <- fixef(object)
     if (missing(parm)) parm <- names(beta)
     se <- sqrt(diag(vcov(object)))
     z <- qnorm((1+level)/2)
     ctab <- cbind(beta-z*se,beta+z*se)
     colnames(ctab) <- stats:::format.perc(c((1-level)/2,(1+level)/2),
                                           digits=3)
     return(ctab[parm,])
 }
 confint(r)
 ##                  2.5 %    97.5 %
 ## (Intercept) 235.575485 269.27845
 ## Days          9.211197  12.04305

• if you think the Wald intervals are likely to be nonsymmetric on the original scale, it might be possible to compute confidence intervals on another scale, e.g. the log scale, and back-transform.
• I don't think that profile likelihood confidence intervals are an option; at the very least you'd have to follow through the theory for robust linear models and see if there was a robust-likelihood analogue that followed the same asymptotic theory.
• the same problem applies for parametric bootstrapping.
• nonparametric bootstrapping is a possibility. See e.g. Confidence intervals on predictions for a non-linear mixed model (nlme) , Non-linear mixed model (nlme) with nested random effect, do not know how to include nested random effect when bootstrapping

Answer 2

It is also now possible to obtain confidence intervals of rlmerMod objects using the effects package.

Example:

library(robustlmm)
library(effects)

r <- rlmer(Reaction ~ Days + (1|Subject), sleepstudy)
as.data.frame(effect("Days",r))

###   Days      fit       se    lower    upper
###      0 252.4270 8.597854 235.4601 269.3938
###      2 273.6812 8.161895 257.5747 289.7877
###      4 294.9355 7.967757 279.2120 310.6589
###      7 326.8168 8.161895 310.7103 342.9234
###      9 348.0711 8.597854 331.1042 365.0379