[R] Very small random effect estimation in lmer but not in stata xtmixed

Joerg Luedicke joerg.luedicke at gmail.com
Fri May 4 02:37:42 CEST 2012


I just realized that I have sent this only to Mohammed. So here it is
to the list:

These two model fits should yield the same results. In fact, if we use
some simulated data generated by the model

yik=2+0.5*xik1+0.25*xik2+uk+eik , with uk~N(0, 1) and eik~N(0, 1)

and compare the results between Stata and R:

. xtmixed y xik1 xik2 || id:, reml

Mixed-effects REML regression                   Number of obs      =      5000
Group variable: id                              Number of groups   =      1000

                                               Obs per group: min =         5
                                                              avg =       5.0
                                                              max =         5


                                               Wald chi2(2)       =   2083.89
Log restricted-likelihood = -8013.0388          Prob > chi2        =    0.0000

------------------------------------------------------------------------------
          y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       xik1 |   .4708663    .025753    18.28   0.000     .4203914    .5213411
       xik2 |   .2726184   .0256147    10.64   0.000     .2224145    .3228222
      _cons |   2.007218   .0354583    56.61   0.000     1.937721    2.076715
------------------------------------------------------------------------------

------------------------------------------------------------------------------
 Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  sd(_cons) |   1.028595    .027431      .9762119    1.083789
-----------------------------+------------------------------------------------
               sd(Residual) |   .9980663   .0111614      .9764285    1.020184
------------------------------------------------------------------------------

Linear mixed model fit by REML
Formula: y ~ xik1 + xik2 + (1 | id)
  Data: data3
  AIC   BIC logLik deviance REMLdev
 16036 16069  -8013    16009   16026
Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 1.05801  1.02859
 Residual             0.99614  0.99807
Number of obs: 5000, groups: id, 1000

Fixed effects:
           Estimate Std. Error t value
(Intercept)  2.00722    0.03546   56.61
xik1         0.47087    0.02575   18.28
xik2         0.27262    0.02561   10.64

Correlation of Fixed Effects:
    (Intr) xik1
xik1 -0.007
xik2  0.005 -0.798

we can see that the results are _exactly_ the same, both for fixed and
random effects. The correlation of fixed effects, of which results are
coming along when using the -summary- function in R are irrelevant for
the model fit. In the above data, variables xik1 and xik2 are drawn
from a multivariate normal distribution with r=0.8.

Thomas, OPs model contains only varying intercepts so there are no
correlations of random effects in his model.

I don't know what went wrong. Mohammed, how are your variables
measured? Maybe huge differences in scale or something like that are
causing a problem in one of the packages. Perhaps try to center or
standardize your independent variable and see if that helps.

Also note that in general, postings that are concerning
multilevel/mixed-effects models might have a better home over at the R
mixed model forum
(https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models).

Joerg

On Thu, May 3, 2012 at 4:50 AM, Mohammed Mohammed
<M.A.MOHAMMED at bham.ac.uk> wrote:
> Hi folks
>
> I am using the lmer function (in the lme4 library) to analyse some data where individuals are clustered into sets (using the SetID variable) with a single fixed effect (cc - 0 or 1). The lmer model and output is shown below.
> Whilst the fixed effects are consistent with stata (using xtmixed, see below), the std dev of the random effect for SetID is very very small (3.5803e-05)compared to stata's (see below 1.002). Any ideas why this should be happening please....?
>
>
> LMER MODEL
>
> summary(lmer(AnxietyScore ~ cc + (1|SetID), data=mydf))
> Linear mixed model fit by REML
> Formula: AnxietyScore ~ cc + (1 | SetID)
>   Data: mydf
>   AIC   BIC logLik deviance REMLdev
> 493.4 503.4 -242.7    486.6   485.4
> Random effects:
> Groups   Name        Variance   Std.Dev.
>  SetID    (Intercept) 1.2819e-09 3.5803e-05
> Residual             1.3352e+01 3.6540e+00
> Number of obs: 90, groups: SetID, 33
>
> Fixed effects:
>            Estimate Std. Error t value
> (Intercept)   3.1064     0.5330   5.828
> cc            2.3122     0.7711   2.999
>
> Correlation of Fixed Effects:
>   (Intr)
> cc -0.691
>
>
> STATA XTMIXED
>
> xtmixed anxietyscore cc || setid:, reml
>
> Mixed-effects REML regression                   Number of obs      =        90
> Group variable: setid                           Number of groups   =        33
> Log restricted-likelihood = -242.48259          Prob > chi2        =    0.0023
> ------------------------------------------------------------------------------
> anxietyscore |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
> -------------+----------------------------------------------------------------
>          cc |   2.289007   .7492766     3.05   0.002     .8204519    3.757562
>       _cons |   3.116074   .5464282     5.70   0.000     2.045094    4.187053
> ------------------------------------------------------------------------------
> ------------------------------------------------------------------------------
>  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
> -----------------------------+------------------------------------------------
> setid: Identity              |
>                   sd(_cons) |   1.002484    .797775      .2107137    4.769382
> -----------------------------+------------------------------------------------
>                sd(Residual) |   3.515888   .3281988      2.928045     4.22175
> ------------------------------------------------------------------------------
>
>
> with thanks & best wishes
> M
>
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



More information about the R-help mailing list