[R-sig-ME] extract level 2 residuals of merMod from lme() and test for spatial autocorrelation.

[Ben Bolker]

  Briefly (without spending too much time digging into the details), I
suspect the discrepancy between the variance of the residuals and the
reported residual variance may be that the residual is calculated
based on the *penalized* residual sum of squares ... ?  May be worth
taking a look at the formulas in the lme4 JSS paper (which is included
as a vignette in the package) to clarify/confirm.


On Wed, May 11, 2016 at 9:32 AM, Dexter Locke <dexter.locke at gmail.com> wrote:
> 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
>
>