How to combine confidence intervals for a variance component of a mixed-effects model when using multiple imputation The logic of multiple imputation (MI) is to impute the missing values not once but several (typically M=5) times, resulting in M completed datasets. The M completed datasets are then analyzed with complete-data methods upon which the M estimates and their standard errors are combined using Rubin's formulas to obtain the "overall" estimate and its standard error.
Great so far, but i'm not sure how to apply this recipe when variance components of a mixed-effects model are concerned. The sampling distribution of a variance component is asymmetrical - therefore the corresponding confidence interval can't be given in the typical "estimate ± 1.96*se(estimate)" form. For this reason the R packages lme4 and nlme don't even provide the standard errors of the variance components, but only provide confidence intervals.
We can therefore perform MI on a dataset and then get M confidence intervals per variance component after fitting the same mixed-effect model on the M completed datasets. The question is how to combine these M intervals into one "overall" confidence interval.
I guess this should be possible - the authors of an article (yucel & demirtas (2010) Impact of non-normal random effects on inference by MI) seem to have done it, but they don't explain exactly how.
Any tips would be much obliged!
Cheers, Rok
 A: This is a great question! Not sure this is a full answer, however, I drop these few lines in case it helps.
It seems that Yucel and Demirtas (2010) refer to an older paper published in the JCGS, Computational strategies for multivariate linear mixed-effects models with missing values, which uses an hybrid EM/Fisher scoring approach for producing likelihood-based estimates of the VCs. It has been implemented in the R package mlmmm. I don't know, however, if it produces CIs.
Otherwise, I would definitely check the WinBUGS program, which is largely used for multilevel models, including those with missing data. I seem to remember it will only works if your MV are in the response variable, not in the covariates because we generally have to specify the full conditional distributions (if MV are present in the independent variables, it means that we must give a prior to the missing Xs, and that will be considered as a parameter to be estimated by WinBUGS...). It seems to apply to R as well, if I refer to the following thread on r-sig-mixed, missing data in lme, lmer, PROC MIXED. Also, it may be worth looking at the MLwiN software.
A: Repeated comment from above:
i'm not sure that a proper analytical solution to this problem even exists. I've looked at some additional literature, but this problem is elegantly overlooked everywhere. I've also noticed that Yucel & Demirtas (in the article i mentioned, page 798) write:

These multiply imputed datasets were
  used to estimate the model […] using
  the R package lme4 leading to 10
  sets of (beta, se(beta)), (sigma_b,
  se(sigma_b)) which were then combined
  using the MI combining rules defined
  by Rubin.

It seems they used some kind of shortcut to estimate the SE of the variance component (which is, of course, inappropriate, since the CI is asymmetrical) and then applied the classic formula.
A: Disclaimer: This idea might be foolish &  I'm not going to pretend to understand the theoretical implications of what I'm proposing.
"Suggestion"  :  Why don't you simply impute 100 (I know you normally do 5) datasets, run the lme4 or nmle, get the confidence intervals (you have 100 of them) and then:
Using a small interval width (say range / 1000 or something), test over the range of possible values of each parameter and include only those small intervals which appear in at least 95 of the 100 CIs.  You would then have a Monte Carlo "average" of your confidence intervals.
I'm sure there are issues (or perhaps theoretical problems) with this approach.  For instance, you could end up with a set of disjoint intervals.  This may or may not be a bad thing depending on your field. Note that this is only possible if you have at least two completely non-overlapping confidence intervals which are separated by a region with less than 95% coverage.
You might also consider something closer to the Bayesian treatment of missing data to get a posterior credible region which would certainly be better formed & more theoretically support than my ad-hoc suggestion.
