Computing standard errors (SEs) and confidence intervals for BLUPs in lmer/lme4

Useful email exchange over here. Key bits—something like this will work:

rfg <- ranef(my.lmer, postVar=TRUE)
my.se = sqrt(as.numeric(attributes(rfg$group)$postVar))

But note Dimitris Rizopoulos’s remark:

“… maybe one thing that I think needs to be kept in mind is that the posterior variances that ranef(…, postVar = TRUE) returns condition on the MLEs and do not take their variability into account.”

Douglas Bates notes over here:

“… I prefer to call the quantities that are returned by the ranef extractor ‘the conditional modes of the random effects’. If you want to be precise, these are the conditional modes (for a linear mixed model they are also the conditional means) of the random effects B given Y = y, evaluated at the parameter estimates. One can also evaluate the condtional variance-covariance of B given Y = y and hence obtain a prediction interval. These are returned by ranef when the optional argument ‘postVar’ is TRUE…

“BTW, the reason that I say ‘conditional modes’, rather than ‘conditional means’, is so the term can apply to generalized linear mixed models, nonlinear mixed models and generalized nonlinear mixed models. The conditional distribution of B given Y is a continuous multivariate distribution that can be evaluated explicitly for a linear mixed model but is known only up to a constant for the generalized forms of the model. We can still optimize the conditional density of B given Y = y for particular values of the parameters but doing so provides the conditional modes, not the conditional means.

Goldstein (see e.g., his book Multilevel Statistical Models) calculates that multiplying by 1.39 rather than 1.96 allows pairwise comparison of these estimates.

Also the arm package has a function se.ranef. Code for that below:

function (object)
{
    se.bygroup <- ranef(object, postVar = TRUE)
    n.groupings <- length(se.bygroup)
    for (m in 1:n.groupings) {
        vars.m <- attr(se.bygroup[[m]], "postVar")
        K <- dim(vars.m)[1]
        J <- dim(vars.m)[3]
        se.bygroup[[m]] <- array(NA, c(J, K))
        for (j in 1:J) {
            se.bygroup[[m]][j, ] <- sqrt(diag(as.matrix(vars.m[,
                , j])))
        }
        names.full <- dimnames(se.bygroup)
        dimnames(se.bygroup[[m]]) <- list(names.full[[1]], names.full[[2]])
    }
    return(se.bygroup)
}
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s