¿Por qué lmer (a diferencia de SAS) muestra una varianza distinta de cero del efecto aleatorio si también se incluye como fijo: (1|sitio) + sitio?

Question

(Estoy editando completamente mi post. La cuestión de fondo es la misma, pero ahora estoy publicando datos simulados para ilustrar mi problema).

He ejecutado (creo) el mismo modelo tanto en proc mixed como en R (lmer) pero he obtenido resultados ligeramente diferentes.

He simplificado mi problema de un post anterior. En este ejemplo, he creado un conjunto de datos procedentes de dos sitios (10 observaciones de cada uno). Me gustaría controlar el centro como efecto aleatorio.

En un principio, quería controlar otros efectos fijos, pero lo que no me resultó obvio en aquel momento fue que los efectos fijos que introduje en el modelo identificaban de forma exclusiva cada lugar. En esencia, estaba ajustando un modelo que tenía el lugar como efecto fijo y aleatorio. Para simplificar en este artículo, simplemente voy a ajustar un modelo que tiene el lugar como efecto fijo y aleatorio.... nunca haría esto en la práctica, pero simplemente estoy ilustrando cómo cada programa maneja esta situación.

Cuando ejecuto este modelo en R, no se imprimen errores y R calcula una varianza positiva/grande para el sitio. Cuando ejecuto el mismo modelo (creo) en SAS, se calcula una varianza 0 para el sitio. La varianza 0 de SAS tiene sentido, ya que no debería quedar variabilidad después de controlar completamente el lugar como efecto fijo.

Estos son los datos:

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

Aquí está el código R y la salida:

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

Salida:

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

Aquí está el código SAS y la salida:

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

Salida:

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

Lo que es aún más impar es que todas las demás estimaciones son idénticas.

De nuevo, se trata de una generalización a partir de un conjunto de datos reales. El modelo anterior tiene un problema obvio (nunca querría ajustar un término como fijo y aleatorio a la vez)... pero usando R, no sabría que hay un problema. SAS destaca que hay 0 varianza de este modelo.

Actualización:

Muchos señalan un mensaje de error en la salida de R. A mí no me aparece este mensaje de error. Es decir el ¡problema! He estado utilizando estos modelos sin ningún indicio de que algo fuera mal. Estoy utilizando R versión 3.4.3 (2017-11-30) RStudio versión 1.0.143 lme4_1.1-15

El modelo es obviamente tonto, pero SAS da dos indicaciones de que algo está mal....a 0 estimación de la varianza y una Nota sobre el hessiano. Actualmente no estoy recibiendo tal error en R.

ACTUALIZACIÓN:

Aquí está el final de la salida de 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

Otra actualización:

En primer lugar, quiero dar las gracias por toda la ayuda y las aportaciones recibidas.

Como Ben Bolker señaló más adelante, el problema número uno es que por alguna razón no estoy recibiendo un mensaje de error. Me gustaría seguir trabajando en R. ¿Hay algo más que pueda ejecutar/comprobar "manualmente" para asegurarme de que todo funciona correctamente? Estaba pensando que tal vez podría salida de la matriz hessiana y / o su determinante para verificar que es positivo. No he encontrado nada en Internet sobre cómo hacerlo. Todo lo que pude encontrar es la salida de la matriz de varianza / covarianza de los efectos fijos ... o sólo los efectos aleatorios.

Además, he calculado el REML "a mano" para entender mejor cuál era el problema. El código está abajo, y parece que la función REML es de alguna manera invariante al efecto aleatorio del sitio. Por eso tanto SAS como R producen resultados idénticos para todas las demás covariables. Técnicamente, ambos conjuntos de resultados son equivalentes, aunque, para mí, el resultado de SAS es más intuitivo. La varianza 0 es también, como otros han señalado, el resultado de la estimación ML.

Código R para la función REML:

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
}

Salida para distintos valores de la varianza del emplazamiento:

> 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

El comportamiento por defecto para lme4::lmer es utilizar la estimación de máxima verosimilitud restringida. Si cambia a la estimación de máxima verosimilitud "pura", el resultado de las réplicas de varianza 0:

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

Creo que puedes encontrar más información sobre la especificación de la ecuación de máxima verosimilitud restringida en: Bates, D., & DebRoy, S. (2004). Linear mixed models and penalized least squares. Revista de Análisis Multivariante , 91(1), 1-17. https://doi.org/10.1016/j.jmva.2004.04.013

Answer 2

Tus ediciones aclaran mucho la cuestión. Muchas gracias.

Sería mejor si lmer dio un valor de cero para el componente de varianza del sitio, pero eso supone que la acción apropiada es siempre declarar que el componente de varianza del sitio es cero en lugar de advertir que hay un conflicto lógico entre la especificación de los efectos fijos y los efectos aleatorios.

lmer hace dos advertencias muy fuertes:

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

En Model failed to converge... advertencia debería ser advertencia suficiente.

Actualización:

Cuando corro str(l) Al final de esa salida obtengo lo siguiente:

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

No sé por qué su salida no muestra eso.

Answer 3

Comprueba tu registro SAS. He ejecutado su código y obtengo

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

Así que SAS está fallando de la misma manera que R, pero los resultados se reportan de manera diferente.

Answer 4

El modelo SAS que especifica es un modelo de pendiente aleatoria. Si ejecuta

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

(aparte, REML es el valor por defecto)

obtendrá una estimación distinta de cero en la salida

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

En el registro aparece la advertencia sobre el hessiano:

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