extract level 2 residuals of merMod from lme() and test for spatial autocorrelation.
Thanks very much. Yes, that was a typo. I'm working with lme4 and the lmer function to make merMod objects. Sorry for any confusion. Here is a reproducible example of my attempt to extract residuals at the population level (ie random effects) by subtracting predictions from observations. library(lme4) fm1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy) # load library, fit merMod object test <- sleepstudy$Reaction - (predict(fm1, re.form= ~ (1| Subject))) # subtract predictions from observations summary(fm1) # examine residuals of random effects, which provides: # "Residual 960.5 30.99 " for variance and std dev, respectively var(test); sd(test) # which provides 869.8171 and 29.49266 # test <- sleepstudy$Reaction - (predict(fm1, re.form= NULL)) # unsurprisingly has same results as above 960.5 <> 869.8171 and 30.99 <> 29.49266 although they are the same direction and order of magnitude. I'm finding similar results with my actual data and model fits. Any advice is appreciated, thanks again! Dexter
On Sat, May 7, 2016 at 2:03 PM, Ben Bolker <bbolker at gmail.com> wrote:
Dexter Locke <dexter.locke at ...> writes:
Apologies for cross-posting, another version of this was posted on r-sig-geo at ... I want to extract the level 2 residuals from a merMod object created with lme() and am struggling.
Not sure what you mean here -- this is a little bit inconsistent. lmer (from lme4) makes merMod objects, lme (from nlme) makes lme objects. I'm going to assume you're talking about lmer (since that's a more likely typo).
How can I pull out a vector that is the residuals at level 2? Everything I see on stack exchange and the mixed models list serve points towards methods like VarCorr(). I'm not interested in a summary of the L2 residuals so much as a vector of the same length as the number of groups I have. Ultimately I want to test them for spatial autocorrelation. Apologies if this is documented somewhere, but I have
not
been able to find help online.
I always get confused about what levels mean in the context of mixed models (i.e. whether one counts down from the population level or up from the individual level), but the basic answer here is to get the *predictions* at the desired level and subtract them from the observed values, see ?predict.merMod -- the re.form argument is what you need. re.form=NA or re.form=~0 gives you predictions at the population level (i.e. neglecting all random effects), re.form=NULL gives you predictions at the individual level (i.e. including all random effects), intermediate cases can be gotten by using a formula. The results will be ordered the same as the original observation vector.
While I have found methods of incorporating spatial effects into a mixed model using corStruct, I am interested in first evaluating if that is an appropriate model for a given dataset by examining the level 2 residuals' spatial patterning - or lack thereof.
Seems reasonable.
Which slot contains the residuals for level 2 and how are they ordered?
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models