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

Question

Nota : esta pregunta es un repost, ya que mi anterior pregunta tenía que ser eliminado por razones legales.

---

Si bien la comparación de PROC MIXED de SAS con la función lme de la nlme paquete en R, me topé con algunos bastante confuso diferencias. Más específicamente, los grados de libertad en las diferentes pruebas difieren entre PROC MIXED y lme, y me preguntaba por qué.

Inicio de la siguiente conjunto de datos (código R indica a continuación) :

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

La idea es construir los siguientes modelos :

y ~ trt + (ind) : ind como un factor aleatorio y ~ trt + (fac(ind)) : fac anidado en ind como un factor aleatorio

Tenga en cuenta que el último modelo debe causar singularidades, como sólo hay 1 valor de y para cada combinación de ind y fac.

Primer Modelo

En SAS, voy a construir el siguiente modelo :

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

De acuerdo a las tutorías, el mismo modelo en R usando nlme debe 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 de los coeficientes y sus SE, pero cuando la realización de una prueba F para el efecto de la 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

Pregunta1: ¿Cuál es la diferencia entre ambas pruebas? Ambos están equipados uso de REML, y el uso de los mismos contrastes.

NOTA: he intentado con diferentes valores para el DDFM= opción (incluyendo BETWITHIN, que teóricamente deberían dar los mismos resultados que lme)

El 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 debe ser :

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

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

• R encaja sin quejarse, mientras que SAS se observa que el final de hess no es positiva definida (que no me sorprende un poco, ver arriba)
• La SE en los coeficientes difieren (es menor en SAS)
• De nuevo, la prueba F se utiliza una cantidad diferente de DF (de hecho, en SAS que cantidad = 0)

SAS de salida :

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 . 

Salida R :

> 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, la F y la prueba de la T son equivalentes y utilizar el mismo DF.)

Curiosamente, cuando se utiliza lme4 en R el modelo ni siquiera caben :

> 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? Son ellos especificado correctamente y si es así, ¿cómo viene los resultados son tan diferentes?

---

Los 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")

Los 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

Answer 1

Para la primera pregunta, el método por defecto en SAS para encontrar el df no es muy inteligente, se ve términos en el efecto aleatorio que sintácticamente incluir el efecto fijo, y los usos que. En este caso, ya trt no se encuentra en ind, no está haciendo lo correcto. Nunca he intentado BETWITHIN y no conozco los detalles, pero el Satterthwaite opción (satterth) o el uso de ind*trt como el de efectos aleatorios dar 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, tu SAS código no muy a la altura de su código R; sólo tiene un plazo para fac*ind, mientras que el código R tiene un plazo para ambas ind y fac*ind. (Ver los Componentes de Varianza de salida para ver esto). La adición de este da el mismo SE de trt en todos los modelos, tanto en la Q1 y Q2 (0.1892).

Como nota, este es un extraño modelo de ajuste como el fac*ind plazo tiene una observación para cada nivel, por lo que es equivalente a que el término de error. Esto se refleja en el SAS de salida, donde el fac*ind plazo tiene cero de la varianza. Esto es también lo que el mensaje de error de lme4 está diciendo; el motivo del error es que lo más probable es que mal especificada algo como lo que estamos incluyendo el término de error en el modelo de dos maneras diferentes. Curiosamente, hay una ligera diferencia en la nlme modelo; es de alguna manera encontrar una varianza plazo para la fac*ind plazo, además del término de error, pero te darás cuenta de que la suma de estos dos varianzas iguales el término de error de SAS y nlme sin la fac*ind plazo. Sin embargo, la SE de trt sigue siendo el mismo (0.1892) como trt está anidado en ind, por lo que la reducción de la varianza de los términos no afectan.

Por último, una nota general sobre los grados de libertad en estos modelos: son calculadas después de que el modelo es adecuado, y así las diferencias en los grados de libertad entre los diferentes programas o las opciones de un programa no significa necesariamente que el modelo se ajuste de manera diferente. Para eso, uno debe mirar a las estimaciones de los parámetros, tanto los parámetros de efecto fijo y la covarianza de los parámetros.

Además, el uso de la t y la F aproximaciones con un determinado número de grados de libertad es bastante controvertido. No sólo hay varias maneras 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, comparar los resultados con la tradicional método de mínimos cuadrados, ya que deben estar de acuerdo. Si es cerca de equilibrado, calcular usted mismo (suponiendo que el equilibrio), de modo que usted puede asegurarse de que el que está utilizando en el derecho de béisbol.

2. Si usted tiene una muestra de gran tamaño, los grados de libertad no importa mucho, ya que las distribuciones de acercarse a una normal y chi-cuadrado.

3. Echa un vistazo a Doug Bates de los métodos de inferencia. Su método anterior es basado en simulación MCMC; su más reciente método se basa en el perfil de la la probabilidad.