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

Wed May 27 03:51:24 CEST 2015

-----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-----