# 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 very interesting question. I look forward to hearing from your results, if you want to share them... – chl Sep 29 '10 at 19:25
@chl: I can send you the tables with the results when I'm done, but I won't be inventing anything new really. So far I'm just planning to compare MI under a two-level imputation model (R package pan) to MI under a simple normal model (ignoring the two-level structure, R package norm) and listwise deletion. Under different sample sizes, values of the variance component etc. This should be enough for the seminar (I'm a PhD student), but not exactly groundbreaking. If you have any ideas on how to "jazz up" the simulation study, i'd love to hear. – Rok Oct 6 '10 at 7:37
One other thing: 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 looked over 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.” – Rok Oct 6 '10 at 7:40
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. – Rok Oct 6 '10 at 7:47
Ok, thx for that. Can you put your comments into an answer so that it can be voted? – chl Oct 6 '10 at 9:25

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.

-
Thanks a lot for your reply! In principle i'm also interested in how to solve a concrete problem like the one i described (thus thanks for the WinBUGS tip). But at the moment i'm trying to do a simulation study for a seminar paper in which i'd examine the performance (coverage rates etc.) of MI under model misspecification. I suppose i'll just forget about the variance components if I can't find a solution and focus on the fixed effects, but it's frustrating to give up. – Rok Sep 13 '10 at 13:15
@Rok Great idea for the simulation! I'll look forward for this particular issue. I suppose you already search on the r-sig-mixed mailing and Gelman's book on multilevel regression... – chl Sep 13 '10 at 13:20
I looked now, tanks for the references! Unfortunately, there's nothing on MI in the r-sig-mixed archives; and Gelman only gives the basic formula on how to combine inferences from MI when we have variation within and between imputations given (§25.7). – Rok Sep 13 '10 at 15:14

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.

-
I appreciate you came back to share your experience with this problem. Unfortunately, I have no real solution but maybe other suggestions will come up. – chl Oct 6 '10 at 14:30
"Elegantly overlooked"... that's a useful phrase for reviewing the literature if ever I've heard one. – Matt Parker Sep 29 '11 at 16:59

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.

-