4
$\begingroup$

(I am completely editing my post. The root question is the same, but I am now posting simulated data to illustrate my issue).

I ran (I beleive) the same model in both proc mixed and R (lmer) but have gotten slightly different results.

I have simplified my problem from a previous post. In this example, I have created a dataset that came from two sites (10 observations from each). I would like to control for site as a random effect.

Originally, I had additional fixed effects that I wanted to control for, but what was not obvious to me at the time was that the fixed effects I entered into the model uniquely identified each site. In essence, I was fitting a model that had site as both a fixed and random effect. To simplify for this post, I am simply going to fit a model that has site as both a fixed and random effect....you would never do this in practice, but I am simply illustrating how each program handles this situation.

When I run this model in R, no errors are printed and R calculates a positive/large variance for site. When I fit the same (i beleive) model in SAS, a 0 variance estimate is calculated for site. The 0 variance from SAS makes sense, since there should be no variability left over after completely controlling for site as a fixed effect.

Here is the data:

site    y
1   6.67949
1   4.667
1   5.91541
1   6.72867
1   5.52195
1   5.98493
1   5.75202
1   6.99084
1   6.19884
1   7.26799
2   1.78078
2   3.68979
2   2.63699
2   3.37598
2   4.07128
2   2.90417
2   3.0458
2   4.94125
2   4.28889
2   2.88859

Here is the R code and output:

l <- lmer(y ~ (1|site) + site, data=fake_data)
summary(l)

Output:

Linear mixed model fit by REML ['lmerMod']
Formula: y ~ (1 | site) + site
   Data: fake_data

REML criterion at convergence: 49.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8662 -0.5452 -0.1016  0.7029  1.8630 

Random effects:
 Groups   Name        Variance Std.Dev.
 site     (Intercept) 3.7025   1.9242  
 Residual             0.7182   0.8475  
Number of obs: 20, groups:  site, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept)    6.171      1.943   3.176
site2         -2.808      2.747  -1.022

Correlation of Fixed Effects:
      (Intr)
site2 -0.707

Here is the SAS code and output:

proc mixed data=dat;
 class site (ref=first);
 model y = site/solution;
 random intercept / subject=site;
run;

Output:

Covariance Parameter Estimates 
Cov Parm Subject Estimate 
Intercept site   0 
Residual         0.7182 

Fit Statistics 
-2 Res Log Likelihood 49.7 
AIC (Smaller is Better) 51.7 
AICC (Smaller is Better) 52.0 
BIC (Smaller is Better) 50.4 

Solution for Fixed Effects 
Effect    site   Estimate Standard Error DF t Value Pr > |t| 
Intercept         6.1707  0.2680         0  23.03   . 
          site 2 -2.8084  0.3790         0  -7.41   . 
          site 1  0        .             .    .     . 

What is even more odd is that all other estimates are identical.

Again, this is a generalization from a real dataset. The above model has an obvious issue (you would never want to fit a term as both fixed and random)...but using R, you would not know there is a problem. SAS highlights that there is 0 variance from this model.

Update:

Many are pointing out an error message in their R output. I am not getting this error message. That is the issue! I had been running these models without any indication something was wrong. I am using: R version 3.4.3 (2017-11-30) RStudio Version 1.0.143 lme4_1.1-15

The model is obviously silly, but SAS gives two indications something is wrong....a 0 variance estimate and an Note about the Hessian. Currently I am getting no such error in R.

UPDATE:

Here is the end of the output from str(l)

  ..@ optinfo:List of 7
  .. ..$ optimizer: chr "bobyqa"
  .. ..$ control  :List of 1
  .. .. ..$ iprint: int 0
  .. ..$ derivs   :List of 2
  .. .. ..$ gradient: num 3.2e-10
  .. .. ..$ Hessian : num [1, 1] 2.86e-06
  .. ..$ conv     :List of 2
  .. .. ..$ opt : int 0
  .. .. ..$ lme4: list()
  .. ..$ feval    : int 20
  .. ..$ warnings : list()
  .. ..$ val      : num 2.27

Another update:

First, I want to say thanks for all the help and input below.

As Ben Bolker pointed out below, the number one issue is for some reason I am not getting an error message. I would like to continue working in R. Is there anything else I can run/check "manually" to be sure everything is running correctly? I was thinking maybe I could output the Hessian matrix and/or it's determinant to verify it is positive. I haven't found anything online on how to do that. All I could find is outputting the variance/covariance matrix for the fixed effects...or only the random effects.

Also, I went ahead and calculated the REML "by hand" to better wrap my head around what was the issue. The code is below, and it appears the REML function is somehow invariant to the site random effect. That is why both SAS and R produce identical output for all other covariates. Technically both sets of results are equivalent, although, to me, SAS's output is more intuitive. The 0 variance is also, as others have pointed out, the result from the ML estimation.

R code for REML function:

y <- as.matrix(fake_data['y'])
X <- cbind(rep(1, each=20),
           c(rep(0, each=10), rep(1, each=10)))

calc_REML <- function(int,b,s,res,y,X){
  #int = intercept, b = beta for site 2, 
  #s = site variance, res = residual variance,
  #y = output vector, X = fixed effects design matrix

  #Create Sigma Matrix
  smat <- matrix(s, nrow = 10, ncol = 10)
  zmat <- matrix(0, nrow = 10, ncol = 10)
  sig1 <- rbind(cbind(smat, zmat), cbind(zmat, smat))
  sig2 <- diag(res, nrow=20, ncol=20)
  sigma <- sig1 + sig2

  #Create Fix effects matrix
  beta <- as.matrix(c(int, b))
  mu <- X %*% beta

  #Calculate REML
  reml <- -((20 - rankMatrix(X))/2)*log(2*pi) -   
    .5*log(det(sigma)) -
    .5*log(det(t(X) %*% solve(sigma) %*% X)) -
    .5*(t(y - mu) %*% solve(sigma) %*% (y - mu))

  neg2reml <- as.numeric(-2*reml)
  neg2reml
}

Output for various values of site variance:

> calc_REML(6.171,-2.808,3.7025,0.7182,y,X) #R's solution
[1] 49.72954
> calc_REML(6.171,-2.808,0,     0.7182,y,X) #SAS solution
[1] 49.72954
> calc_REML(6.171,-2.808,10000, 0.7182,y,X) #Extreme solution
[1] 49.72954
$\endgroup$

migrated from stackoverflow.com Feb 15 '18 at 1:39

This question came from our site for professional and enthusiast programmers.

  • $\begingroup$ Posting the data somewhere or having code that simulates data that matches your design would be helpful. $\endgroup$ – JimB Feb 14 '18 at 23:44
  • $\begingroup$ In the second set of models why would you have NewSiteID (i.e., the factory) as both a fixed effect and a random effect? And (to me) your coding of the fixed effects in both SAS and R seems a bit odd in that NewSiteID is already a class variable in SAS and a factor (i.e., categorical variable) in R and you don't need to add in 14 dummy variables. Why not run the following in R: lmer(ln_i ~ NewSiteID + (1 | NewSiteID:ID), data=d) ? Then the results match with SAS. And R is giving you a warning that something is off. $\endgroup$ – JimB Feb 15 '18 at 3:13
  • $\begingroup$ I agree that the second set of models I am specifying is wrong/silly. You would never do that in practice (enter Site as both a fixed and random effect). However, if you do, the variance for Site I would expect to be 0....since there should be no leftover variance after fully controlling for it with the fixed effects. $\endgroup$ – itch987 Feb 15 '18 at 14:01
  • $\begingroup$ @amoeba /solution gives the estimates of the fixed effects in the "Solution for Fixed Effects" table. $\endgroup$ – JimB Feb 15 '18 at 15:19
6
$\begingroup$

The default behavior for lme4::lmer is to use restricted maximum likelihood estimation. If you switch to 'pure' maximum likelihood estimation, the result of 0 variance replicates:

l2 <- lmer(y ~ (1|site) + site, data=fake_data, REML = FALSE)
summary(l2)
# Linear mixed model fit by maximum likelihood  ['lmerMod']
# Formula: y ~ (1 | site) + site
#    Data: fake_data
# 
#      AIC      BIC   logLik deviance df.resid 
#       56       60      -24       48       16 
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -1.9671 -0.5747 -0.1071  0.7409  1.9638 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev.
#  site     (Intercept) 0.0000   0.000   
#  Residual             0.6464   0.804   
# Number of obs: 20, groups:  site, 2
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)   8.9791     0.5685  15.794
# site         -2.8084     0.3596  -7.811
# 
# Correlation of Fixed Effects:
#      (Intr)
# site -0.949

I think you can find out more about the specification of the restricted maximum likelihood equation in: Bates, D., & DebRoy, S. (2004). Linear mixed models and penalized least squares. Journal of Multivariate Analysis, 91(1), 1–17. https://doi.org/10.1016/j.jmva.2004.04.013

$\endgroup$
  • $\begingroup$ can the OP confirm that PROC MIXED is using ML by default, and that we get similar results to R's REML output if we switch to REML in SAS? $\endgroup$ – Ben Bolker Feb 15 '18 at 19:22
  • $\begingroup$ @Ben According to the answer by Chris that was just posted here, SAS is using REML but OP did not correctly specify the random intercept model. But anyway, can one intuitively understand why lmer produces non-zero site variance with REML? $\endgroup$ – amoeba Feb 15 '18 at 20:00
  • $\begingroup$ Yes, SAS by default runs REML. The post by Chris was not completely correct. His specfication and mine should be the same since site is a class variable with 2 levels. $\endgroup$ – itch987 Feb 15 '18 at 20:33
  • 4
    $\begingroup$ There are a few issues here. The main one is why lme4 is not reliably producing a warning message! The second one is the intuition behind lmer producing a non-zero variance under REML. The third (which I don't claim to understand) is what the correct specification in SAS is, or at least what specifications are equivalent in R vs SAS. $\endgroup$ – Ben Bolker Feb 15 '18 at 22:07
1
$\begingroup$

Check your SAS log. I ran your code and get

NOTE: Convergence criteria met but final Hessian is not positive definite.
NOTE: Estimated G matrix is not positive definite.
NOTE: PROCEDURE MIXED used (Total process time):
      real time           0.03 seconds
      cpu time            0.03 seconds

So SAS is failing the same way R is failing, but the results are reported differently.

$\endgroup$
1
$\begingroup$

The SAS model you specify is a random slope model. If you run

proc mixed data=fake_data;
 class site (ref=first);
 model y = site/solution;
 random intercept / subject=site;
run;

(aside, REML is the default)

you will get a non-zero estimate in the output

Covariance Parameter Estimates 
Cov Parm Subject Estimate 
Intercept site 0.08978 
Residual       0.7182 

In the log is the warning about the Hessian:

NOTE: Convergence criteria met but final Hessian is not positive definite.
NOTE: PROCEDURE MIXED used (Total process time):
      real time           0.08 seconds
      cpu time            0.04 seconds
$\endgroup$
  • $\begingroup$ that is not what i get when i run your code. Both my specification and yours are equivalent since site is a class variable with only 2 levels. What version of SAS are you running? $\endgroup$ – itch987 Feb 15 '18 at 20:29
  • $\begingroup$ NOTE: Copyright (c) 2016 by SAS Institute Inc., Cary, NC, USA. NOTE: SAS (r) Proprietary Software 9.4 (TS1M5) NOTE: This session is executing on the X64_7PRO platform. NOTE: Updated analytical products: SAS/STAT 14.3 SAS/ETS 14.3 SAS/OR 14.3 SAS/IML 14.3 SAS/QC 14.3 NOTE: Additional host information: X64_7PRO WIN 6.1.7601 Service Pack 1 Workstation $\endgroup$ – Chris Andrews Feb 15 '18 at 21:11
  • $\begingroup$ Covariance Parameter Estimates Cov Parm Subject Estimate Intercept site 0 Residual 0.7182 Fit Statistics -2 Res Log Likelihood 49.7 AIC (Smaller is Better) 51.7 AICC (Smaller is Better) 52.0 BIC (Smaller is Better) 50.4 Solution for Fixed Effects Effect site Estimate Standard Error DF t Value Pr > |t| Intercept 6.1707 0.2680 0 23.03 . site 2 -2.8084 0.3790 0 -7.41 . site 1 0 . . . . $\endgroup$ – itch987 Feb 16 '18 at 14:30
  • $\begingroup$ Sorry, I tried to copy my output from SAS using your code to show I still get the exact same out as my post above. The formatting is lost. $\endgroup$ – itch987 Feb 16 '18 at 14:33
1
$\begingroup$

Your edits make the issue so much clearer. Thank you!

It might be nicer if lmer gave a value of zero for the site variance component but that does assume that the appropriate action is always to state that the site variance component is zero rather than warning that there is a logical conflict between the specification of the fixed effects and the random effects.

lmer does give two very strong warnings:

unable to evaluate scaled gradient
Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

The Model failed to converge... warning should be warning enough.

Update:

When I run str(l) I get the following at the end of that output:

..@ optinfo:List of 7
  .. ..$ optimizer: chr "bobyqa"
  .. ..$ control  :List of 1
  .. .. ..$ iprint: int 0
  .. ..$ derivs   :List of 2
  .. .. ..$ gradient: num -4.26e-10
  .. .. ..$ Hessian : num [1, 1] -2.86e-06
  .. ..$ conv     :List of 2
  .. .. ..$ opt : int 0
  .. .. ..$ lme4:List of 2
  .. .. .. ..$ code    : int -3
  .. .. .. ..$ messages: chr [1:2] "unable to evaluate scaled gradient" "Model failed to converge: degenerate  Hessian with 1 negative eigenvalues"
  .. ..$ feval    : int 20
  .. ..$ warnings : list()
  .. ..$ val      : num 2.51

I don't know why your output does not show that.

$\endgroup$
  • $\begingroup$ I guess i am not seeing those errors in R. Above is all the output from R when i run the summary(). Do I need to ask for something else? $\endgroup$ – itch987 Feb 15 '18 at 16:43
  • $\begingroup$ What version of R and lme4 are you using? Can you see the warnings when you execute str(l) ? $\endgroup$ – JimB Feb 15 '18 at 16:58
  • 1
    $\begingroup$ @JimB I also don't see any warnings with the latest release of R and fully updated CRAN packages. What version of R and lme4 are you using? $\endgroup$ – Ista Feb 15 '18 at 18:26
  • $\begingroup$ Touche! R 3.3.3 and lme4 1.1-12. (I'll load the newest versions and see what happens.) But I'm not using R-Studio. What happens when you run it outside of R-Studio? $\endgroup$ – JimB Feb 15 '18 at 18:51
  • $\begingroup$ In 1.1-15 get the warnings when using other optimizers (like "Nelder_Mead" and "nloptwrap"), but not with the default "bobyqa". $\endgroup$ – aosmith Feb 15 '18 at 18:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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