[R-sig-ME] lme4 vs. HLM (program) - differences in estimation?

[Douglas Bates]

On Tue, Nov 18, 2008 at 3:51 PM, Felix Schönbrodt <nicebread at gmx.net> wrote:
> Thanks very much for your thorough remarks!
>
> I suspect you are not using the same data between HLM and R and that may be
> the problem. That is, in R you create a variable called ses.gmc thinking
> this is a group mean centered variable. But, HLM "automagically" does group
> mean centering for you if you ask it to.
>
> When you work in HLM are you using the exact same data file that you created
> for your use in R? Or, are you relying on HLM to group mean center for you?
> If so, I suspect that is the issue. In fact, we can see that the difference
> in your estimates lies in this variable and this is the variable you
> manipulate in R, so I suspect the error may no be a function of estimation
> differences, but of data differences.
>
> I followed your suggestion about the possibly different centering approach
> in HLM and did a run only with raw variables (no centering). The small
> differences stay.
>
> I would check two things: does HLM estimate two variances and a
> covariance for the random effects and is it using the REML criterion
> or the ML criterion.
>
> HLM is using the REML criterion; concerning the variances for random effects
> I couldn't dig deep enough till now to answer this question ...
> From a lot of comparisons I have run now, I can conclude the following:
> - fixed effects usually are equal up to 3 digits after the decimal point
> - random variance components can show small deviations if they are very
> close to zero (and are not siginificant anyway - following the p-value
> reported by HLM, which is usually > 0.50 in these cases). However, I am
> aware of the discussion concerning the appropriateness of p-values reported
> by HLM (see
> also http://markmail.org/search/?q=lme4+p-value#query:lme4%20p-value+page:1+mid:3t4hxxrlh3uvb7kh+state:results).

That's an interesting point about the biggest apparent differences
being in the small variance estimates.  It has been a while since I
looked at descriptions of the computational methods in HLM but I
believe that they work in terms of the precision of the random effects
(the inverse of the variance) instead of the variance.  In the way
that the model was typically expressed when that software was written
it is easier to work with the precision matrix than with the variance
matrix of the random effects.

When you look at the situation in more detail, however, you find that
it is possible to have estimates of zero for the variances of the
random effects whereas you cannot have an estimate of zero for the
precision, corresponding to infinite variance.  Optimizing the
log-likelihood with respect to the variance parameters (or, better,
the standard deviations) is more stable than optimizing with respect
to the precision.  The instability of optimizing with respect to the
precision parameters is most noticeable for small variance components
(large precision estimates).  It is particularly dramatic when the
variance goes to zero (precision to infinity).  Some of the many,
incompatible changes to the lme4 package - changes that have
inconvenienced many of its users whom I thank for their patience -
have been specifically for the purpose of stabilizing the estimation
situation for very small variance components.

The place where the choice of estimation algorithm matters the most is
in the estimation of small variance components so it doesn't surprise
me that differences crop up there.


> Maybe I was too nitpicking about these small differences - they only seem to
> appear in coefficients/ variance components that are insignificant anyway.
> But the discussion was definitely productive for me, and I think I now can
> convince my colleagues to use lme4 instead ;-)
> Maybe someoner wants to dig deeper into this issue; the data set from the
> HLM program (HS&B data) can be downloaded here for
> example: http://www.hlm-online.com/datasets/HSBDataset/HSBdata.zip.
>
> Thanks,
> Felix
>