[R-sig-ME] [R] different results from lme and lmer function

Wed May 27 16:43:47 CEST 2015

On Tue, May 26, 2015 at 8:52 PM Ben Bolker <bbolker at gmail.com> wrote:

> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> On 15-05-26 09:45 PM, Pantelis Z. Hadjipantelis wrote:
> > Ben,
> >
> > I fully accept that the two procedures should converge at the same
> > solution but as the OP does not give the versions used. Therefore I
> > am not sure that the fitting is done using the same procedure
> > ("bobyqa" vs. "Nelder_Mead"). Different optim. algorithms with
> > different initial values might converge to different local minima.
> > The OP has only 39 subjects so for a model of this size
> > over-fitting will not be unheard off given the number of
> > parameters.
> >
> > I am curious though for your later comment: I understand that in
> > the case of an unstable model you might expect lme4 to be slightly
> > better, but wouldn't the singular fit in the correlation (-1)
> > suggest that lmer's fit is sub-optimal?
> >
> > All best, Pantelis
>
>   No, I suspect that the singular fit actually represents the best
> fit/minimal achievable REML criterion for this problem.  I can't say
> for sure (and of course none of us can without having the original
> data to play with), but based on my previous experiences I do suspect
> that this is just a very flat surface and that lme has stopped
> *slightly* too soon.
>

Exactly.  You may recall the discussion in "Fitting Linear Mixed-Effects
Models with lme4" (http://arxiv.org/abs/1406.5823, some day, I hope within
my lifetime, to appear in jstatsoft.org) regarding solving the penalized
linear least squares problem, that the formulation in lme uses the relative
precision (inverse of the variance) factor and lmer uses the relative
covariance factor.  An important difference between the two forms is that
the lmer formulation can represent a singular model but the lme formulation
can't.

The lmer fit is singular and lme must stop before the singular model.

>
>   Ben Bolker
>
>
> >
> >
> > On 05/26/2015 05:50 PM, Ben Bolker wrote:
> >> I agree that the difference is trivially small/practically
> >> unimportant.  The point here is that -- just for those of us who
> >> are interested in the details of the methods -- lme4 and lme are
> >> fitting *exactly* the same model, by similar methods, so in
> >> general they should converge to the same answer (to a somewhat
> >> closer tolerance than this).  Generally when they don't it's
> >> because the model is slightly unstable, and I have generally
> >> found that lme4 does slightly better (but I wouldn't rule out the
> >> opposite case).
> >>
> >> cheers Ben Bolker
> >>
> >>
> >> On Tue, May 26, 2015 at 8:47 PM, John Sorkin
> >> <jsorkin at grecc.umaryland.edu> wrote:
> >>> Ben, I doubt the very small difference in log likelihood gives
> >>> much, if any information about which model is a better fit.
> >>> Even if we overlook the limited precision of the estimate of
> >>> the REML criterion, the difference is so small as to me of
> >>> minimal importance. John
> >>>
> >>> John David Sorkin M.D., Ph.D.
> >>>
> >>> Professor of Medicine
> >>>
> >>> Chief, Biostatistics and Informatics
> >>>
> >>> University of Maryland School of Medicine Division of
> >>> Gerontology and Geriatric Medicine
> >>>
> >>> Baltimore VA Medical Center
> >>>
> >>> 10 North Greene Street
> >>>
> >>> GRECC (BT/18/GR)
> >>>
> >>> Baltimore, MD 21201-1524
> >>>
> >>> (Phone) 410-605-7119
> >>>
> >>> (Fax) 410-605-7913 (Please call phone number above prior to
> >>> faxing)
> >>>
> >>>
> >>> On May 26, 2015, at 8:03 PM, Ben Bolker <bbolker at gmail.com>
> >>> wrote:
> >>>
> >>> These actually aren't terribly different from each other.  I
> >>> suspect that lmer is slightly closer to the correct answer,
> >>> because lme reports a "log-likelihood" (really -1/2 times the
> >>> REML criterion) of 49.30021, while lmer reports a REML
> >>> criterion of  -98.8 -> slightly better fit at -R/2 = 49.4.  The
> >>> residual sds are 0.0447 (lme) vs. 0.0442 (lmer); the intercept
> >>> sd estimate is 0.016 vs 0.0089, admittedly a bit low, and both
> >>> month sds are very small.  lmer indicates a singular fit
> >>> (correlation of -1).    If you look at the confidence intervals
> >>> on these estimates (confint(fitted_model) in lme4;
> >>> intervals(fitted_model) in lme) I think you'll find that the
> >>> confidence intervals are much wider than these differences (you
> >>> may even find that lme reports that it can't give you the
> >>> intervals because the Hessian [curvature] matrix is not
> >>> positive definite).
> >>>
> >>> Both should be comparable to SAS PROC MIXED results, I think,
> >>> if you get the syntax right ...
> >>>
> >>> On Tue, May 26, 2015 at 7:09 PM, li li <hannah.hlx at gmail.com>
> >>> wrote:
> >>>
> >>> Hi all,
> >>>
> >>> I am fitting a random slope and random intercept model using R.
> >>> I used both lme and lmer funciton for the same model. However I
> >>> got different results as shown below (different variance
> >>> component estimates and so on). I think that is really
> >>> confusing. They should produce close results. Anyone has any
> >>> thoughts or suggestions. Also, which one should be comparable
> >>> to sas results?
> >>>
> >>> Thanks!
> >>>
> >>> Hanna
> >>>
> >>>
> >>> ## using lme function
> >>>
> >>> mod_lme <- lme(ti  ~ type*months, random=~ 1+months|lot,
> >>> na.action=na.omit,
> >>>
> >>> + data=one, control = lmeControl(opt = "optim"))
> >>>
> >>> summary(mod_lme)
> >>>
> >>> Linear mixed-effects model fit by REML
> >>>
> >>> Data: one
> >>>
> >>> AIC       BIC   logLik
> >>>
> >>> -82.60042 -70.15763 49.30021
> >>>
> >>>
> >>> Random effects:
> >>>
> >>> Formula: ~1 + months | lot
> >>>
> >>> Structure: General positive-definite, Log-Cholesky
> >>> parametrization
> >>>
> >>> StdDev       Corr
> >>>
> >>> (Intercept) 8.907584e-03 (Intr)
> >>>
> >>> months      6.039781e-05 -0.096
> >>>
> >>> Residual    4.471243e-02
> >>>
> >>>
> >>> Fixed effects: ti ~ type * months
> >>>
> >>> Value   Std.Error DF   t-value p-value
> >>>
> >>> (Intercept)     0.25831245 0.016891587 31 15.292373  0.0000
> >>>
> >>> type            0.13502089 0.026676101  4  5.061493  0.0072
> >>>
> >>> months          0.00804790 0.001218941 31  6.602368  0.0000
> >>>
> >>> type:months -0.00693679 0.002981859 31 -2.326329  0.0267
> >>>
> >>> Correlation:
> >>>
> >>> (Intr) typPPQ months
> >>>
> >>> type           -0.633
> >>>
> >>> months         -0.785  0.497
> >>>
> >>> type:months  0.321 -0.762 -0.409
> >>>
> >>>
> >>> Standardized Within-Group Residuals:
> >>>
> >>> Min            Q1           Med            Q3           Max
> >>>
> >>> -2.162856e+00 -1.962972e-01 -2.771184e-05  3.749035e-01
> >>> 2.088392e+00
> >>>
> >>>
> >>> Number of Observations: 39
> >>>
> >>> Number of Groups: 6
> >>>
> >>>
> >>>
> >>>
> >>>
> >>> ###Using lmer function
> >>>
> >>> mod_lmer <-lmer(ti  ~ type*months+(1+months|lot),
> >>> na.action=na.omit, data=one)
> >>>
> >>> summary(mod_lmer)
> >>>
> >>> Linear mixed model fit by REML t-tests use Satterthwaite
> >>> approximations to
> >>>
> >>> degrees of freedom [merModLmerTest]
> >>>
> >>> Formula: ti ~ type * months + (1 + months | lot)
> >>>
> >>> Data: one
> >>>
> >>>
> >>> REML criterion at convergence: -98.8
> >>>
> >>>
> >>> Scaled residuals:
> >>>
> >>> Min      1Q  Median      3Q     Max
> >>>
> >>> -2.1347 -0.2156 -0.0067  0.3615  2.0840
> >>>
> >>>
> >>> Random effects:
> >>>
> >>> Groups   Name        Variance  Std.Dev.  Corr
> >>>
> >>> lot      (Intercept) 2.870e-04 0.0169424
> >>>
> >>> months      4.135e-07 0.0006431 -1.00
> >>>
> >>> Residual             1.950e-03 0.0441644
> >>>
> >>> Number of obs: 39, groups:  lot, 6
> >>>
> >>>
> >>> Fixed effects:
> >>>
> >>> Estimate Std. Error        df t value Pr(>|t|)
> >>>
> >>> (Intercept)     0.258312   0.018661  4.820000  13.842 4.59e-05
> >>> ***
> >>>
> >>> type         0.135021   0.028880  6.802000   4.675  0.00245 **
> >>>
> >>> months          0.008048   0.001259 11.943000   6.390 3.53e-05
> >>> ***
> >>>
> >>> type:months -0.006937   0.002991 28.910000  -2.319  0.02767 *
> >>>
> >>> ---
> >>>
> >>> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> >>>
> >>>
> >>> Correlation of Fixed Effects:
> >>>
> >>> (Intr) typPPQ months
> >>>
> >>> type     -0.646
> >>>
> >>> months      -0.825  0.533
> >>>
> >>> type:month  0.347 -0.768 -0.421
> >>>
> >>>
> >>> _______________________________________________
> >>>
> >>> R-sig-mixed-models at r-project.org mailing list
> >>>
> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >>>
> >>>
> >>> ______________________________________________
> >>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more,
> >>> see 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.
> >>>
> >>>
> >>> Confidentiality Statement:
> >>>
> >>> This email message, including any attachments, is for
> >>> ...{{dropped:7}}
> >> _______________________________________________
> >> 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
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.11 (GNU/Linux)
>
> iQEcBAEBAgAGBQJVZSMcAAoJEOCV5YRblxUH+gMH/iqfJA2s56KxPMfON+liOtIy
> cvkOk+OO+xBjDfcz9axQ8i8t7RwihONcR02rCmzrpqlTfOjSQqHvwcTuancQVDEg
> /QpMfamxn2t3+ah6VWCDg+gB6uhaGnUOPc+06V7rLw920bYBLuN67s74n1xGCnC2
> f+wa1lzZ+YLbuZl73Q5KQwaKZmVX+BqDc0eyN2CGtjOZO1A0n7Wewc5IeOHPjh5k
> QH0tynG7/Amv3YLGWLA3BYsoiYOW9AxGJy1FOKfY0XuwWm3+VsOEn8RFw15yCU8b
> BbVRVbmrv6+Rjy3QOF83wWRHKTtAlb96+yxAt4Qi7Cazuplya5MchHNSmZyAFmg=
> =8Jtz
> -----END PGP SIGNATURE-----
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>

	[[alternative HTML version deleted]]