Dear all,
I'm trying to understand how to extract BLUPs of the random effects from
LMMs, and GLMMs, fitted with MCMCglmm.
To use a simple example:
library(lme4)
library(MCMCglmm)
data(PlodiaPO)
m2 <- lmer(PO ~ 1 + (1 | FSfamily), PlodiaPO)
MC2 <- MCMCglmm(PO ~ 1, random = ~ FSfamily, data=PlodiaPO, verbose=F, pr=T)
The help file for "ranef" in the lme4 package says this function
"extract[s] the conditional modes of the random effects from a fitted model
object. For linear mixed models the conditional modes of the random effects
are also the conditional means."
Bracketing some conceptual fuzziness, threads like <
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2007q1/000128.html>
suggest "ranef" yields BLUPs.
But how do I get the equivalent from an MCMCglmm object (like MC2 here)?
Two possibilities:
colMeans(MC2$Sol[,-1])
or
posterior.mode(MC2$Sol[,-1])
In the thread <
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2012q1/017565.html>,
Jarrod says that, with some conditions, "The marginal posterior modes of
the random effects should coincide with BLUPs...". So does that mean
"posterior.mode" is the function to use here?
And given the statement about modes and means in the statement about ranef
above, alternatively, would colMeans work with MCMCglmm in the case of a
LMM, but not GLMM?
In this simple example, the three approaches yield very similar results:
cor(cbind(ranef(m2)$FSfamily, colMeans(MC2$Sol[,-1]),
posterior.mode(MC2$Sol[,-1])))
But maybe that doesn't always happen.
(I also came across the thread <
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/015947.html>, but
that raises some broader questions of principle which I'll park for now.)
Any assistance would be greatly appreciated.
Many thanks,
Malcolm
[[alternative HTML version deleted]]