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

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

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.