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

[Malcolm Fairbrother]

Hi Dexter,
I may be wrong, but it sounds like you may simply be looking for ranef(fm1)
?
However, it seems you are also troubled by the fact that the standard
deviation of the residuals is not the same as the residual-level standard
deviation.
That is, sigma(fm1) and as.data.frame(VarCorr(fm1))[2,] are not equal to
sd(resid(fm1)).
What is resid(fm1)? Compare:
head(resid(fm1))
head(sleepstudy$Reaction - (predict(fm1, re.form= NULL)))
head(with(sleepstudy, Reaction - fixef(fm1)[1] - fixef(fm1)[2]*Days -
ranef(fm1)$Subject[Subject,]))
In other words, you are looking at the residuals *after* taking into
account the random effects (i.e., the group-level errors).
In that vein, compare also:
as.data.frame(VarCorr(fm1))[1,]
sd(ranef(fm1)$Subject[,1])
The standard deviation of the random effects (group-level errors) is,
therefore, not the same as the group-level standard deviation. At both
levels, the variance of the random errors you can extract is smaller than
the variance of the random errors you can't see. So what's going on?
Basically, the random effects are shrunk towards zero, while the estimate
provided of the group-level standard deviation is pre-shrinkage. And
because the residuals you can extract are calculated taking into account
the random effects, the residuals--and thus their SD--are also affected by
the shrinkage.
Does that help?
- Malcolm



> Date: Wed, 11 May 2016 10:30:12 -0400
> From: Ben Bolker <bbolker at gmail.com>
> To: Dexter Locke <dexter.locke at gmail.com>
> Cc: "r-sig-mixed-models at r-project.org"
>         <r-sig-mixed-models at r-project.org>
> Subject: Re: [R-sig-ME] extract level 2 residuals of merMod from lme()
>         and test for spatial autocorrelation.
>
>   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?
>

	[[alternative HTML version deleted]]