Diferencias entre PROC Mixed y lme / lmer en R - grados de libertad

12

Nota: esta pregunta es una nueva publicación, ya que mi pregunta anterior tuvo que ser eliminada por razones legales.


Al comparar PROC MIXED de SAS con la función lmedel nlmepaquete en R, me topé con algunas diferencias bastante confusas. Más específicamente, los grados de libertad en las diferentes pruebas difieren entre PROC MIXEDy lme, y me preguntaba por qué.

Comience desde el siguiente conjunto de datos (el código R se proporciona a continuación):

  • ind: factor que indica al individuo donde se toma la medición
  • fac: órgano donde se toma la medida
  • trt: factor que indica el tratamiento
  • y: alguna variable de respuesta continua

La idea es construir los siguientes modelos simples:

y ~ trt + (ind): indcomo factor aleatorio y ~ trt + (fac(ind)): facanidado indcomo factor aleatorio

Tenga en cuenta que el último modelo debería causar singularidades, ya que solo hay 1 valor de ypor cada combinación de indy fac.

Primer modelo

En SAS, construyo el siguiente modelo:

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s;
    RANDOM ind /s;
run;

Según los tutoriales, el mismo modelo en R usando nlmedebería ser:

> require(nlme)
> options(contrasts=c(factor="contr.SAS",ordered="contr.poly"))
> m2<-lme(y~trt,random=~1|ind,data=Data)

Ambos modelos dan las mismas estimaciones para los coeficientes y su SE, pero cuando realizan una prueba F para el efecto de trt, utilizan una cantidad diferente de grados de libertad:

SAS : 
Type 3 Tests of Fixed Effects 
Effect Num DF Den DF     F  Value Pr > F 
trt         1      8  0.89        0.3724 

R : 
> anova(m2)
            numDF denDF  F-value p-value
(Intercept)     1     8 70.96836  <.0001
trt             1     6  0.89272  0.3812

Pregunta 1: ¿Cuál es la diferencia entre ambas pruebas? Ambos están equipados con REML y usan los mismos contrastes.

NOTA: Probé diferentes valores para la opción DDFM = (incluido BETWITHIN, que teóricamente debería dar los mismos resultados que lme)

Segundo modelo

En SAS:

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s;
    RANDOM fac(ind) /s;
run;

El modelo equivalente en R debería ser:

> m4<-lme(y~trt,random=~1|ind/fac,data=Data)

En este caso, hay algunas diferencias muy extrañas:

  • R encaja sin quejarse, mientras que SAS señala que el hessian final no es definitivo positivo (lo cual no me sorprende un poco, ver arriba)
  • El SE en los coeficientes difiere (es menor en SAS)
  • Nuevamente, la prueba F usó una cantidad diferente de DF (de hecho, en SAS esa cantidad = 0)

Salida SAS:

Effect     trt Estimate Std Error  DF t Value Pr > |t| 
Intercept        0.8863    0.1192  14    7.43 <.0001 
trt       Cont  -0.1788    0.1686   0   -1.06 . 

R salida:

> summary(m4)
...
Fixed effects: y ~ trt 
               Value Std.Error DF   t-value p-value
(Intercept)  0.88625 0.1337743  8  6.624963  0.0002
trtCont     -0.17875 0.1891855  6 -0.944840  0.3812
...

(Tenga en cuenta que en este caso, las pruebas F y T son equivalentes y usan el mismo DF).

Curiosamente, cuando se usa lme4en R, el modelo ni siquiera encaja:

> require(lme4)
> m4r <- lmer(y~trt+(1|ind/fac),data=Data)
Error in function (fr, FL, start, REML, verbose)  : 
  Number of levels of a grouping factor for the random effects
must be less than the number of observations

Pregunta 2 : ¿Cuál es la diferencia entre estos modelos con factores anidados? ¿Se especifican correctamente y, de ser así, cómo es que los resultados son tan diferentes?


Datos simulados en R:

Data <- structure(list(y = c(1.05, 0.86, 1.02, 1.14, 0.68, 1.05, 0.22, 
1.07, 0.46, 0.65, 0.41, 0.82, 0.6, 0.49, 0.68, 1.55), ind = structure(c(1L, 
2L, 3L, 1L, 3L, 4L, 4L, 2L, 5L, 6L, 7L, 8L, 6L, 5L, 7L, 8L), .Label = c("1", 
"2", "3", "4", "5", "6", "7", "8"), class = "factor"), fac = structure(c(1L, 
1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L), .Label = c("l", 
"r"), class = "factor"), trt = structure(c(2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Cont", 
"Treat"), class = "factor")), .Names = c("y", "ind", "fac", "trt"
), row.names = c(NA, -16L), class = "data.frame")

Datos simulados:

   y ind fac   trt
1.05   1   l Treat
0.86   2   l Treat
1.02   3   l Treat
1.14   1   r Treat
0.68   3   r Treat
1.05   4   l Treat
0.22   4   r Treat
1.07   2   r Treat
0.46   5   r  Cont
0.65   6   l  Cont
0.41   7   l  Cont
0.82   8   l  Cont
0.60   6   r  Cont
0.49   5   l  Cont
0.68   7   r  Cont
1.55   8   r  Cont
Joris Meys
fuente
@ Aaron: Encuentra tu respuesta incluida en esta publicación. Si pudieras copiar y pegar eso como respuesta, te doy el representante. Ha sido muy útil, por lo que realmente quiero mantenerlo aquí en validación cruzada. Después de que hayas hecho eso, borro tu respuesta de la pregunta.
Joris Meys
Estoy tratando de hacer que el equipo reviva tu Q original con esta desafortunada revisión eliminada para siempre, por lo tanto, hay una gran oportunidad de restaurar las respuestas originales y fusionarlas aquí.
@mbq: Eso sería bueno, aunque simulé algunos datos (que uso aquí) y edité la respuesta de Aaron en consecuencia. Para la otra respuesta, eso será un poco más complicado, pero también puedo intentarlo.
Joris Meys
La respuesta de Aaron es increíblemente buena. Espero que lo vean. Desafortunadamente, tu @Aaron no lo contactará a menos que haya participado en este hilo.
Wayne
1
Sí, esta fue una buena respuesta. Aquí le di un enlace a la publicación eliminada: stats.stackexchange.com/questions/26556/… Voy a agregar el enlace a la publicación actual.
Stéphane Laurent

Respuestas:

11

Para la primera pregunta, el método predeterminado en SAS para encontrar el df no es muy inteligente; busca términos en el efecto aleatorio que incluyan sintácticamente el efecto fijo, y lo usa. En este caso, como trtno se encuentra en ind, no está haciendo lo correcto. Nunca lo he intentado BETWITHINy no conozco los detalles, pero la opción Satterthwaite ( satterth) o el uso ind*trtcomo efecto aleatorio dan resultados correctos.

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s ddfm=satterth;
    RANDOM ind /s;
run;

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s;
    RANDOM ind*trt /s;
run;

En cuanto a la segunda pregunta, su código SAS no coincide con su código R; solo tiene un término para fac*ind, mientras que el código R tiene un término para ambos indy fac*ind. (Consulte la salida de Componentes de varianza para ver esto.) Agregar esto proporciona el mismo SE para trttodos los modelos en Q1 y Q2 (0.1892).

Como observa, este es un modelo extraño para ajustarse ya que el fac*indtérmino tiene una observación para cada nivel, por lo que es equivalente al término de error. Esto se refleja en la salida de SAS, donde el fac*indtérmino tiene variación cero. Esto también es lo que te dice el mensaje de error de lme4; La razón del error es que probablemente haya especificado algo incorrectamente al incluir el término de error en el modelo de dos maneras diferentes. Curiosamente, hay una ligera diferencia en el modelo nlme; es de alguna manera encontrar un término de varianza para el fac*indtérmino además del término de error, pero notará que la suma de estas dos variaciones es igual al término de error de SAS y nlme sin el fac*indtérmino. Sin embargo, el SE para trtpermanece igual (0.1892) como trtestá anidado enind, por lo que estos términos de menor varianza no lo afectan.

Finalmente, una nota general sobre los grados de libertad en estos modelos: se calculan después de que el modelo se ajusta, por lo que las diferencias en los grados de libertad entre diferentes programas u opciones de un programa no necesariamente significan que el modelo se ajuste de manera diferente. Para eso, uno debe mirar las estimaciones de los parámetros, tanto los parámetros de efectos fijos como los parámetros de covarianza.

Además, usar las aproximaciones t y F con un número dado de grados de libertad es bastante controvertido. No solo hay varias formas de aproximar el df, algunos creen que la práctica de hacerlo no es una buena idea de todos modos. Un par de palabras de consejo:

  1. Si todo está equilibrado, compare los resultados con el método tradicional de mínimos cuadrados, como deberían coincidir. Si está cerca del equilibrio, calcule usted mismo (asumiendo el equilibrio) para asegurarse de que los que está usando estén en el estadio correcto.

  2. Si tiene un tamaño de muestra grande, los grados de libertad no importan mucho ya que las distribuciones se acercan a un valor normal y chi-cuadrado.

  3. Vea los métodos de inferencia de Doug Bates. Su método más antiguo se basa en la simulación MCMC; Su método más reciente se basa en el perfil de la probabilidad.

Aaron dejó Stack Overflow
fuente
De hecho, es una buena respuesta, aunque creo que perfilar la probabilidad resuelve una pregunta diferente (IC apropiados en los parámetros de varianza donde el perfil no es cuadrático) que hacer la simulación MCMC (que maneja tanto la corrección de tamaño finito como la no cuadrática). Creo que bootMer (bootstrap paramétrico) está más cerca del equivalente para mcmcsamp que confint (profile (...)) ...
Ben Bolker
@BenBolker: Claro que podría ser. Doug Bates dio una charla aquí el mes pasado y habló sobre sus ideas sobre el perfil de la probabilidad. Eso es todo lo que sé hasta ahora.
Aaron dejó Stack Overflow el