[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