5
$\begingroup$

I got completely different results from lmer() and lme()! Just look at the coefficients' std.errors. Completely different in both cases. Why is that and which model is correct?

> mix1c = lmer(logInd ~ 0 + crit_i + Year:crit_i + (1 + Year|Taxon), data = datai)
> mix1d = lme(logInd ~ 0 + crit_i + Year:crit_i, random = ~ 1 + Year|Taxon, data = datai)
> 
> summary(mix1d)
Linear mixed-effects model fit by REML
 Data: datai 
       AIC      BIC    logLik
  4727.606 4799.598 -2351.803

Random effects:
 Formula: ~1 + Year | Taxon
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev       Corr  
(Intercept) 9.829727e-08 (Intr)
Year        3.248182e-04 0.619 
Residual    4.933979e-01       

Fixed effects: logInd ~ 0 + crit_i + Year:crit_i 
                 Value Std.Error   DF   t-value p-value
crit_iA      29.053940  4.660176   99  6.234515  0.0000
crit_iF       0.184840  3.188341   99  0.057974  0.9539
crit_iU      12.340580  5.464541   99  2.258301  0.0261
crit_iW       5.324854  5.152019   99  1.033547  0.3039
crit_iA:Year -0.012272  0.002336 2881 -5.253846  0.0000
crit_iF:Year  0.002237  0.001598 2881  1.399542  0.1618
crit_iU:Year -0.003870  0.002739 2881 -1.412988  0.1578
crit_iW:Year -0.000305  0.002582 2881 -0.118278  0.9059
 Correlation: 
             crit_A crit_F crit_U crit_W cr_A:Y cr_F:Y cr_U:Y
crit_iF       0                                              
crit_iU       0      0                                       
crit_iW       0      0      0                                
crit_iA:Year -1      0      0      0                         
crit_iF:Year  0     -1      0      0      0                  
crit_iU:Year  0      0     -1      0      0      0           
crit_iW:Year  0      0      0     -1      0      0      0    

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-6.98370498 -0.39653580  0.02349353  0.43356564  5.15742550 

Number of Observations: 2987
Number of Groups: 103 
> summary(mix1c)
Linear mixed model fit by REML 
Formula: logInd ~ 0 + crit_i + Year:crit_i + (1 + Year | Taxon) 
   Data: datai 
  AIC  BIC logLik deviance REMLdev
 2961 3033  -1469     2893    2937
Random effects:
 Groups   Name        Variance   Std.Dev. Corr   
 Taxon    (Intercept) 6.9112e+03 83.13360        
          Year        1.7582e-03  0.04193 -1.000 
 Residual             1.2239e-01  0.34985        
Number of obs: 2987, groups: Taxon, 103

Fixed effects:
               Estimate Std. Error t value
crit_iA      29.0539403 18.0295239   1.611
crit_iF       0.1848404 12.3352135   0.015
crit_iU      12.3405800 21.1414908   0.584
crit_iW       5.3248537 19.9323887   0.267
crit_iA:Year -0.0122717  0.0090916  -1.350
crit_iF:Year  0.0022365  0.0062202   0.360
crit_iU:Year -0.0038701  0.0106608  -0.363
crit_iW:Year -0.0003054  0.0100511  -0.030

Correlation of Fixed Effects:
            crit_A crit_F crit_U crit_W cr_A:Y cr_F:Y cr_U:Y
crit_iF      0.000                                          
crit_iU      0.000  0.000                                   
crit_iW      0.000  0.000  0.000                            
crit_iA:Yer -1.000  0.000  0.000  0.000                     
crit_iF:Yer  0.000 -1.000  0.000  0.000  0.000              
crit_iU:Yer  0.000  0.000 -1.000  0.000  0.000  0.000       
crit_iW:Yer  0.000  0.000  0.000 -1.000  0.000  0.000  0.000
>
Improve this question
$\endgroup$
4
  • $\begingroup$ The random effects are also completely different. How does the output compare when you don't include the interaction terms, or simplfy it even further ? $\endgroup$
    – Robert Long
    Oct 4, 2012 at 6:48
  • $\begingroup$ @longrob - the interaction terms are necessary - these are parameters I must get (= group specific slope). The only way to simplify is to remove the random effect for slope and use only (1|Taxon), in which case lmer and lme yield the same results. But, this is slightly different model. $\endgroup$
    – Tomas
    Oct 4, 2012 at 7:31
  • $\begingroup$ Can you post data/reproducible example? $\endgroup$
    – Ben Bolker
    Oct 4, 2012 at 23:14
  • $\begingroup$ I am still interested in a reproducible example for this question ... $\endgroup$
    – Ben Bolker
    Jun 17, 2016 at 2:44

2 Answers 2

Reset to default
7
$\begingroup$

lmer uses Laplace approximation, when the whole normal distribution of the random effect is approximated at its mode. This approximation is known to produce the estimates of the variance components that are biased down. lme uses a more thorough approximation via Gaussian quadrature approximation, but I neither know the default number of integration points nor the way to manipulate this number.

Note also that lmer produced a correlation of the random effects (the intercept vs. year) of -1. That's very bad, and is indicative of numeric problems. It hit some sort of a ridge in its approximation of the likelihood that it could not overcome (no wonder, given that this approximation is quite poor). lme did a somewhat better job, and came up with a correlation of 0.619. If you really have a year, as in, 2010, 2011, 2012, that's a very bad idea for mixed models, where high multicollinearity between factors may mean poor numeric stability and long convergence. You would have been much better off converting it to time centered near zero, rather than near 2010.

Improve this answer
$\endgroup$
4
  • $\begingroup$ But the lmer results have higher standard errors than lme - opposite of what you wrote? Anyway, the differences are so big that one of the results must be wrong, OK? $\endgroup$
    – Tomas
    Oct 3, 2012 at 19:59
  • $\begingroup$ I looked at it again, and found another problem -- see update. $\endgroup$
    – StasK
    Oct 5, 2012 at 2:14
  • 3
    $\begingroup$ (1) The first paragraph of this answer is false ... lmer does not use a Laplace approximation (this, or Gaussian quadrature, is/are only necessary for GLMMs/glmer). The J Stat Software paper on lmer (also included as a vignette with the package) explains this in great detail. (2) I agree that correlations of -1 (singular fits) are bad, but that doesn't by itself indicate that lmer is right and lme is wrong - singular fits are sometimes the "best" solution to the particular statistical question. (3) I agree that centering is a good idea. $\endgroup$
    – Ben Bolker
    Jun 17, 2016 at 2:55
  • $\begingroup$ Thanks, @BenBolker -- would you care the edit my answer to rectify? $\endgroup$
    – StasK
    Jun 21, 2016 at 19:14
0
$\begingroup$

The log-likelihood is substantially higher for the lmer() result, which would suggest it's closer to the true optimum.

Since the lmer() fit is on the boundary of the parameter space, which lme() reparametrises off to infinity (the 'Log-Cholesky' parametrisation), it's plausible that lme just didn't find the better solution.

That doesn't entirely answer the question of which answer is more useful, but it does argue against the lme one.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.