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

Pantelis Z. Hadjipantelis kalakouentin at gmail.com
Wed May 27 05:19:16 CEST 2015 

Hmm... interesting rationale.
Thank you for your insights on this computational matter; they are 
invaluable.

Pantelis

On 26/05/15 18:51, Ben Bolker 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.
>
>    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

More information about the R-sig-mixed-models mailing list