3
$\begingroup$

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.

$\endgroup$
6
$\begingroup$
  • 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
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for your answer. I already did that myself, but I am not so happy with Wald confidence intervals in this situation, as the normality assumption is obviously violated and my distribution is not symmetric... $\endgroup$ – Steffi Sep 8 '16 at 11:12
  • $\begingroup$ see edits. A non-symmetric response distribution doesn't imply that the Wald CIs are bad ... $\endgroup$ – Ben Bolker Sep 8 '16 at 14:27
  • $\begingroup$ Thank you for expanding your answer! I was also thinking about bootstrapping, but I was/am unsure how to do that correctly in this situation, incorporating the random effects etc. I will have a closer look at your links and see if I manage! It would be interesting also for comparison with the Wald confidence intervals. About your suggestion to calculate the CIs on another scale: There I see the problem with the interpretation. The model is for a surgery trial and the surgeons need to be able to interpret the results - many of them do not have an in-depth knowledge about statistics. $\endgroup$ – Steffi Sep 9 '16 at 6:32
  • $\begingroup$ Hi all! @BenBolker, I am using robustlmm and I will use your code to calculate confidence intervals. I was wondering if I have to give some specific reference if I use your method to calculate the confidence intervals since I am not using any package. By the way, I read a lot of comments from you everywhere, and I have to tell you THANK YOU VERY MUCH for sharing your knowledge. $\endgroup$ – Dekike Apr 25 at 13:49
  • 1
    $\begingroup$ I think that some generic reference to Wald confidence intervals would be sufficient ... $\endgroup$ – Ben Bolker Apr 25 at 15:56
3
$\begingroup$

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
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ one doubt, How should I specified in this package if I want confidence intervals for two fixed factors and their interaction? When I specified the interaction term I get multiple rows for each level-combination between the fixed factors, but I don't want that. I want just the C.I for each fixed factor and the interaction. Do you know how to do this? I checked help for the package but I can't get what I want. If I use the code proposed by Ben Bolker, I get what I want, but I wanted to check if with the package I could get the same. Thanks! $\endgroup$ – Dekike Apr 25 at 14:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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