Interval estimates for multilevel models

Zhongyi Yuan writes,

I have a problem about linear mixed effect model.

Generally, when we want to do prediction in a linear mixed effect model, we obtain BLUPs of random effects first, then a point prediction. But how can we get an interval prediction , like in an ordinary linear model? It seems lme or lmer can’t handle that directly (right? maybe xtmixed can? It seems Stata is more powerful than R in dealing with mixed effect models. However, I am not sure. I seldom use Stata. ). While considering constructing prediction intervals like Bonferroni intervals( seems feasible if we could obtain the distribution of bias ), we need the estimated Correlation Matrix, which is given by R only in balanced models. Unfortunately, my model is unbalanced.

Besides, cuz I am coping with an econometric problem, I have to specify the within-group correlation structure. The default correlation structure in lme() in NULL, which is not the case. (By the way, it seems we can’t specify that in lmer(). What is the advantage of package lme4 over package nlme then). But the StdDev grows rapidly after a correlation structure specification. Prediction becomes inaccurate. I wonder if it makes any sense for us to do prediction in these circumstances.

My quick response: First off, lmer() does give standard errors (se.fixef() and se.ranef()). More generally you can get intervals using simulation, for example using mcsamp() (from the “arm” package). If you want arbitrary covariance matrices, you might have to program it yourself, although it’s not too hard if you go inside lm() and augment the data matrix. We discuss the algebra of it in the regression chapter of Bayesian Data Analysis. Certainly there’s no statistical reason why this can’t or shouldn’t be done.

This entry was posted in Uncategorized by Andrew. Bookmark the permalink.

2 thoughts on “Interval estimates for multilevel models

  1. I have a problem that may be easy to some of you, but I just can't find an answer to. The problem is about Bonferroni bounds on the probability of the union of events. These bounds are build on the marginals (let's call them Pi) and the pairwise joint probabilities Pij. In some contexts Bonferroni bounds can be too loose, but if the events are clusterized, meaning that each event ei is strongly correlated to only a subset of events, then they are pretty tight (at least that's how they come out in my application).

    Now, assume that we can compute the marginals Pi exactly, and that we have a die that generates outcomes with the right mutual correlations amongst the vents.

    The question is, is there an algorithm to estimate all Pij, starting from the marginals (as priors), with increasing accuracy as the number of realizations grow?

    For example, one could start with the estimates Pij ≈ Pi*Pj, then through the die and update those, through the die again, etc. I'd like to know the confidence of the estimates of Pij as we go along.

    As a final piece of information, I may add that in typical applications (I'm looking at), Pij can be very small, so just rolling the die and estimating Pij from averages may take a very long time do accomplish with acceptable accuravy.

  2. Estimating Pairwise Joint Probabilities from Marginals and Experiments —

    I have a problem that may be easy to some of you, but I just can't find an answer to. The problem is about Bonferroni bounds on the probability of the union of events. These bounds are build on the marginals (let's call them Pi) and the pairwise joint probabilities Pij. In some contexts Bonferroni bounds can be too loose, but if the events are clusterized, meaning that each event ei is strongly correlated to only a subset of events, then they are pretty tight (at least that's how they come out in my application).

    Now, assume that we can compute the marginals Pi exactly, and that we have a die that generates outcomes with the right mutual correlations amongst the vents.

    The question is, is there an algorithm to estimate all Pij, starting from the marginals (as priors), with increasing accuracy as the number of realizations grow?

    For example, one could start with the estimates Pij ≈ Pi*Pj, then through the die and update those, through the die again, etc. I'd like to know the confidence of the estimates of Pij as we go along.

    As a final piece of information, I may add that in typical applications (I'm looking at), Pij can be very small, so just rolling the die and estimating Pij from averages may take a very long time do accomplish with acceptable accuravy.

Comments are closed.