Why does lmer (unlike SAS) output a non-zero variance of random effect if it is also included as fixed: (1|site) + site?

Question

(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

Answer 1

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

Answer 2

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.

Answer 3

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

Answer 4

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.