I've fit an lmer model with the following (albeit made up output):

Random effects:
 Groups        Name        Std.Dev.
 day:sample (Intercept)    0.09
 sample        (Intercept) 0.42
 Residual                  0.023 

I'd really like to build a confidence interval for each effect using the following formula:

$ \frac{(n-1)s^2}{\chi^2_{\alpha/2, n-1}},\frac{(n-1)s^2}{\chi^2_{1-\alpha/2,n-1}} $

Is there a way to conveniently get out the degrees of freedom?

  • 1
    I think you want to check lmerTest. There are a number of approximations to approximate the d.f. in a mixed-effects model for the fixed effects (eg. Satterthwaite, Kenward-Roger, etc.) For your case it seems to me that you overcomplicating your life. You do assume each effect to be Gaussian. Just use the standard deviation to get the confidence interval of your choice. – usεr11852 Apr 18 '15 at 0:00
  • 3
    @usεr11852 In a mixed-effects model you assume that each effect is Gaussian, but the parameter is the variance of the Gaussian distribution, not the mean. So the distribution of its estimator will be very skewed, and the normal ± ~2 standard deviations confidence interval will not be appropriate. – Karl Ove Hufthammer Apr 18 '15 at 6:58
  • 1
    @KarlOveHufthammer: You are right to point this out; I see what you (and probably the OP) mean. I thought he was concerned about the means and/or the realizations of the random effects as he mentioned degrees of freedom. – usεr11852 Apr 18 '15 at 9:11
  • degrees of freedom are "problematic" for mixed-models, see: stat.ethz.ch/pipermail/r-help/2006-May/094765.html and stats.stackexchange.com/questions/84268/… – Tim Apr 18 '15 at 9:56
up vote 13 down vote accepted

I would instead just create profile likelihood confidence intervals. They’re reliable, and very easy to calculate using the ‘lme4’ package. Example:

> library(lme4)
> fm = lmer(Reaction ~ Days + (Days | Subject),
            data=sleepstudy)
> summary(fm)
[…]
Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 612.09   24.740       
          Days         35.07    5.922   0.07
 Residual             654.94   25.592       

You can now calculate the profile likelihood confidence intervals with the confint() function:

> confint(fm, oldNames=FALSE)
Computing profile confidence intervals ...
                               2.5 %  97.5 %
sd_(Intercept)|Subject        14.381  37.716
cor_Days.(Intercept)|Subject  -0.482   0.685
sd_Days|Subject                3.801   8.753
sigma                         22.898  28.858
(Intercept)                  237.681 265.130
Days                           7.359  13.576

You can also use the parametric bootstrap to calculate confidence intervals. Here’s the R syntax (using the parm argument to restrict which parameters we want confidence intervals for):

> confint(fm, method="boot", nsim=1000, parm=1:3)
Computing bootstrap confidence intervals ...
                              2.5 % 97.5 %
sd_(Intercept)|Subject       11.886 35.390
cor_Days.(Intercept)|Subject -0.504  0.929
sd_Days|Subject               3.347  8.283

The results will naturally vary somewhat for each run. You can increase nsim to decrease this variation, but this will also increase the time it takes to estimate the confidence intervals.

  • Good answer (+1). I would also mention the fact that one can also get CIs from parametric bootstrap in this case. This thread contains a very interesting discussion on the matter. – usεr11852 Apr 18 '15 at 9:14
  • @usεr11852 Thanks for the suggestion. I have now added an example using the parametric bootstrap. – Karl Ove Hufthammer Apr 18 '15 at 9:32

Degrees for freedom for mixed-models are "problematic". For reading more on it you can check the lmer, p-values and all that post by Douglas Bates. Also r-sig-mixed-models FAQ summarizes the reasons why it is bothersome:

  • In general it is not clear that the null distribution of the computed ratio of sums of squares is really an F distribution, for any choice of denominator degrees of freedom. While this is true for special cases that correspond to classical experimental designs (nested, split-plot, randomized block, etc.), it is apparently not true for more complex designs (unbalanced, GLMMs, temporal or spatial correlation, etc.).
  • For each simple degrees-of-freedom recipe that has been suggested (trace of the hat matrix, etc.) there seems to be at least one fairly simple counterexample where the recipe fails badly.
  • Other df approximation schemes that have been suggested (Satterthwaite, Kenward-Roger, etc.) would apparently be fairly hard to implement in lme4/nlme,
    (...)
  • Because the primary authors of lme4 are not convinced of the utility of the general approach of testing with reference to an approximate null distribution, and because of the overhead of anyone else digging into the code to enable the relevant functionality (as a patch or an add-on), this situation is unlikely to change in the future.

The FAQ gives also some alternatives

  • use MASS::glmmPQL (uses old nlme rules approximately equivalent to SAS 'inner-outer' rules) for GLMMs, or (n)lme for LMMs
  • Guess the denominator df from standard rules (for standard designs) and apply them to t or F tests
  • Run the model in lme (if possible) and use the denominator df reported there (which follow a simple 'inner-outer' rule which should correspond to the canonical answer for simple/orthogonal designs), applied to t or F tests. For the explicit specification of the rules that lme uses, see page 91 of Pinheiro and Bates — this page is available on Google Books
  • use SAS, Genstat (AS-REML), Stata?
  • Assume infinite denominator df (i.e. Z/chi-squared test rather than t/F) if number of groups is large (>45? Various rules of thumb for how large is "approximately infinite" have been posed, including [in Angrist and Pischke's ''Mostly Harmless Econometrics''], 42 (in homage to Douglas Adams)

But if you are interested in confidence intervals there are better approaches, e.g. based on bootstrap as suggested by Karl Ove Hufthammer in his answer, or the ones proposed in the FAQ.

  • "Guess the denominator df from standard rules (for standard designs) and apply them to t or F tests"; I'd really like if someone could elaborate on that. For instance, for a nested design (f. ex. patients vs controls, several samples per subject; with subject ID being the random effect), how do we get the degrees of freedom for such a design? – Arnaud A Jan 24 at 0:05

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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