[R-sig-ME] Extracting Standard Errors of Uncorrelated Random Effects?

[Reinhold Kliegl]

May I add one consideration that guides my approach in model building
and ask a question? I usually delete non-significant correlation
parameters and variance components with the argument that my current
set of data is most likely not rich enough to support stable estimates
of these model parameters. From this perspective, can I expect that
the simple model yields more stable estimates of the parameters than
the full model? The fact that I drop non-significant parameters does
not imply that I accept the null hypothesis about them. In other
words, I consider it very likely that with a larger, more reliable set
of data than the present one I would be able to estimate these
parameters (and keep them in the model.) Therefore, in a way, I expect
model complexity (and theoretical impact) to grow as I improve the
reliability of my measures or increase the data base.

Reinhold Kliegl

On Wed, Dec 14, 2011 at 10:34 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
> On Wed, Dec 14, 2011 at 2:42 PM, Derek Dunfield <dunfield at mit.edu> wrote:
>> Dr. Bates:
>> Many thanks for the quick response - this was exactly what I was looking for.
>
>> An interesting aside - in Andrew Gelman's "Data Analysis Using
>> Regression and Multilevel/Hierarchical Models", which does a great job
>> outlining multilevels models using R, he doesn't seem to mention
>> uncorrelated models at all (other than a special case where all groups
>> share the same intercept). In your opinion, should I always aim to use
>> an uncorrelated fit if there is no significant difference between it
>> and the correlated model?
>
> There are differences of opinion on this.  I have spent the last 30+
> years in the Department of Statistics at the University of Wisconsin -
> Madison, a department that was founded by George Box.  To me it is
> natural to pursue parsimony in a model, which means that I delete
> terms that do not appear to be contributing significantly to the model
> fit.  Like the quote attributed to Einstein, "Make things as simple as
> possible, but no simpler."
>
> In fact, in this particular case it is of interest to test the
> hypothesis of no correlation of the random effects because the
> experimenters want to know if they can predict the extent of sleep
> deprivation's effect on response time from the initial response time.
> (I.e., are those who have faster response times initially less
> affected by sleep deprivation?)
>
> Others (feel free to chime in here, Ben) believe that "model building"
> by deleting apparently insignificant terms results in overfitting of
> the model and I don't dispute that.  In some ways it depends on what
> the objective in fitting the model is.  For the purposes of prediction
> I want to pay attention to the bias-variance trade-off and aim for a
> simple, adequate model.  For the purposes of establishing the
> significance of a fixed-effects term, preliminary simplification of
> the model may bias the effect of the model.
>
> The part of the "stay with the initial model regardless" approach that
> I don't like is that I am not convinced that the initial model is
> necessarily a good model.
>
>
>> Many thanks again for your excellent work on lme4.
>> Derek
>> --
>> Derek Dunfield, PhD
>> Postdoctoral Fellow, MIT Intelligence Initiative
>> Sloan School of Management
>> MIT Center for Neuroeconomics, Prelec Lab
>> 77 Massachusetts Ave E62-585
>> Cambridge MA 02139
>>
>>
>>
>> On Wed, Dec 14, 2011 at 3:16 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
>>> On Wed, Dec 14, 2011 at 1:05 PM, Derek Dunfield <dunfield at mit.edu> wrote:
>>>> Hello all:
>>>> I was hoping to extract SE for my random effects to get a better idea
>>>> of their precision for each grouping factor.
>>>>
>>>> Andrew Gelman has developed a nice package called "Arm" which does
>>>> this seamlessly for models with correlated random effects such as
>>>> Reaction ~ 1 + Days + (1 + Days | Subject)
>>>>
>>>> When I try to run this on an uncorrelated model such as
>>>> Reaction ~ 1 + Days + (1 | Subject) + (0 + Days | Subject)
>>>>
>>>> I get the error "Error in .local(object, ...) : Code not written yet"
>>>
>>> Unfortunately, it is complicated.  Strangely enough, the case of
>>> correlated random effects is easier than the case of uncorrelated
>>> effects.  The values that you want to get are the solutions to a
>>> system like L x = cbind(e_i, e_j) where the indices i and j are the
>>> random effects corresponding to the k'th level of the grouping factor
>>> and L is the sparse Cholesky factor.  For correlated random effects
>>> the i and j are going to be adjacent.  For uncorrelated random effects
>>> they may end up being adjacent after the fill-reducing permutation but
>>> that is not guaranteed.
>>>
>>> In this particular case things work out fine.  If I fit that model as
>>> fm1 then I can check that the random effects for each subject are
>>> adjacent by examining
>>>
>>> image(fm1 at L)
>>>
>>> where we can see that each of the pairs of blocks starting with 1:2
>>> are correlated within the block but orthogonal to the other blocks.
>>> So, setting
>>>
>>> L <- fm1 at L
>>> Id <- Diagonal(x = rep.int(1, ncol(L))   # Identity matrix of the same size as L
>>>
>>>  the conditional variance-covariance matrix of the orthogonal random
>>> effects for subject 1 is something like
>>>
>>> crossprod(solve(L, Id[,1:2], sys="L")) * lme4:::sigma(fm1)^2
>>>
>>> Before using such an expression I would check it against the results
>>> from the "Arm" package for the model with correlated random effects.
>>> It is entirely possible that I have left something out of the
>>> calculation.
>>>
>>> If you wonder why something that looks simple like this is not part of
>>> the lme4 package, it is because a lot of the simplicity comes from the
>>> fact that both random-effects terms have the same grouping factor.
>>> The general case is much messier.
>>>
>>>> I believe this may be due to some unfinished code in the lme4 package
>>>> for ranef, though I'm not sure.
>>>>
>>>> Does anyone have an idea how I could go about calculating these standard errors?
>>>>
>>>> Many thanks in advance for you help!
>>>> :Derek
>>>>
>>>> --
>>>> Derek Dunfield, PhD
>>>> Postdoctoral Fellow, MIT Intelligence Initiative
>>>> Sloan School of Management
>>>> MIT Center for Neuroeconomics, Prelec Lab
>>>> 77 Massachusetts Ave E62-585
>>>> Cambridge MA 02139
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models