Modelos mixtos / resultados REML: ¿Tenemos documentación sobre cómo se calculará el error estándar?Al usar SLS con REML, también conocido como efectos mixtos, un cliente obtiene diferentes resultados en diferentes paquetes (R / Python / JMP) para el error estándar para el intervalo de confianza. (R / Python / JMP: 0.05 / 0.08 / 0.002) este es un factor de 40, que es bastante alto. ¿Tenemos alguna documentación sobre cómo se calcula exactamente (fórmula)? El BLUPS en sí, así como las estimaciones tanto para FE como para ME, son el sema en los diferentes paquetes, pero el intervalo de predicción y el intervalo de confianza.
Pensé que podría tener algo que ver con una estructura de covarianza repetida subyacente diferente como predeterminada o la corrección de Kacker-Harville.
Añadiendo un ejemplo: https://www.rdocumentation.org/packages/doBy/versions/4.6-3/topics/data-dietox
- modelo de construcción Peso = Tiempo + Alimento + Cerdo (efecto aleatorio)
@jiancao @laural_jmp
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?Hola martin,
¿Puede compartir los datos y el script del modelo?Eso sería muy útil si puedes.Además, ¿qué errores estándar e intervalos de confianza son diferentes?También estoy intentando agregar a Elizabeth a esta pista.
Gracias.
Laura
@eclaassen
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?Martín,
Si se puede proporcionar un código R que produzca los intervalos de confianza y los errores estándar, sería útil. Me estoy quedando atrás en mi conocimiento de R, y solo he podido construir con éxito el modelo (usando los datos dietox y lmer ()) y obtener los valores predichos marginales y condicionales. Esos valores coinciden con JMP.
Utilizando los datos de dietox, ejecuté el modelo especificado en JMP y SAS (PROC MIXED), y los valores predichos marginales y condicionales y los intervalos de confianza coinciden (cuando el método DDFM correspondiente se especifica en SAS - KR2 para marginal y Satterthwaite para condicional) . SAS no tiene CI individuales. Hubo un problema con los grados de libertad para los CI individuales que se ha corregido en 16.0.
Sé que la matriz de covarianza corregida de Kackar-Harville se usa para pruebas / intervalos para estimaciones de efectos fijos solamente. El método Satterthwaite se utiliza si se trata de un efecto aleatorio.
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?El error estándar de 0,002 lo calcula el usuario a partir de las columnas de CI media condicional guardadas.
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?Originalmente había agregado las capturas de pantalla, pero no aparecieron en la publicación. Aquí lo intento de nuevo:
2.125 como Stderr para los componentes de varianza
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?Pido el código R.
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?Con respecto a los grados de libertad, acabamos de tener este pensamiento, y si usan JMP 15, ¿podría ser esta la razón del resultado diferente? Cuando ejecuto este ejemplo en 16 y 15, veo exactamente los mismos resultados.
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?@jiancao ,
¿No veo esta columna de fórmula? Me interesaría mucho verlo, porque no sé cómo lo calcularías sin los grados de libertad, que no son sencillos.
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?Hola Elisabeth,
adjunto ahora también está el código R. Usaron
checkintervals_dietox.R
* R 3.6.3
¿Cualquier otra cosa que necesites? Parece que también usaron lmer.
Re: Según Martin estaban comparando el error estándar ...Ugh. Acababa de escribir una respuesta larga y, cuando llegué a publicar, recibí un error de autenticación y se perdió la respuesta. Intentaré recrearlo ... y espero que no obtengamos una publicación duplicada.
¿Esa cantidad no es la SE? La diferencia = obs + crit * se - (obs - crit * se) = 2 * crit * se. Tomar la varianza de eso no será el error estándar que nos interesa. Sería la variación de la amplitud de los intervalos de confianza, que hasta donde yo sé no es nada "interesante". Si obtiene el error estándar de los valores predichos o los límites de confianza, obtiene un SE ~ 0.834, que se parece más a un SE promedio porque cada observación tiene su propio SE.
La discrepancia entre SAS y JMP es menor. Las probabilidades de (-2 log) son las mismas, por lo que nos encontramos con un poco de óptimos locales o algo de confusión numérica debido a que los valores en la tabla de datos no están bien redondeados. Incluso los valores predichos varían después del cuarto decimal a pesar de que los efectos fijos son los mismos.
/* outpm gives the marginal prediction information (values, se, CI) that corresponds to JMPs Mean Confidence Interval */
proc mixed data=sasuser.dietox nobound cl covtest;
class pig;
model weight= time feed / ddfm=kr2 s outpm=krpredmean;
random pig ;
run;
/* outp gives the conditional prediction information that corresponds to JMPs Conditional Mean CI */
proc mixed data=sasuser.dietox nobound;
class pig;
model weight= time feed / ddfm=satterth s outp=spred;
random pig;
run;
Realmente, cualquier discrepancia en el valor predicho de SE entre JMP y Python / R se debe a la diferencia en el SE para la varianza del grupo / cerdo, como se indica en las capturas de pantalla de Martin. La variación de la variación entra en juego. No sé cómo lo está calculando Python, y parece extraño que todos los demás SE estén de acuerdo y que uno no. Nuevamente, SAS y JMP esencialmente están de acuerdo y JMPs SLS-REML y Mixed están de acuerdo mientras usan rutinas ligeramente diferentes.
Re: Según Martin estaban comparando el error estándar ...Ah, la información perdida en el repost:
Los SE de predicción condicional están disponibles en Mixto de JMPPro. No sé por qué no están disponibles en SLS-REML cuando ambas plataformas ofrecen el SE marginal (Std Error of Predicted).
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?Parece que cualquier discrepancia en los intervalos de confianza proviene de los errores estándar calculados de los componentes de la varianza.
La salida de Python tiene la varianza de grupo / cerdo como 31.836 con SE de 2.125. (De la captura de pantalla de Martin)
La salida R es:
Random effects:
Groups Name Variance Std.Dev.
Pig (Intercept) 31.836 5.642
Residual 7.364 2.714
Number of obs: 789, groups: Pig, 72
JMP tiene:
SAS tiene:
| Estimaciones de parámetros de covarianza |
|---|
| Cov Parm |
Estimar |
Estándar
Error |
Valor Z |
Pr Z |
Alfa |
Inferior |
Superior |
|---|
| Cerdo |
31.8276 |
5.4936 |
5.79 |
<.0001 |
0,05 |
21.0604 |
42.5949 |
|---|
| Residual |
7.3645 |
0.3898 |
18,89 |
<.0001 |
0,05 |
6.6565 |
8.1924 |
|---|
Todos los paquetes obtienen la misma probabilidad de probabilidad.
Las estimaciones de los parámetros son esencialmente las mismas en todos los paquetes.
El SE para el componente de varianza de cerdo / grupo es comparable para SAS, JMP y R, pero muy diferente en Python.
El SE para la varianza residual es comparable para SAS y JMP, pero R es muy diferente y no veo ningún SE para esa estimación en la captura de pantalla de salida de Python.
Debido a que esos SE de los componentes de la varianza son necesarios en los cálculos posteriores, obtenemos diferencias en los intervalos de confianza de la predicción. No veo una manera de obtener los grados de libertad y / o SE que R está usando en los intervalos de confianza para compararlos directamente. Solo vemos las discrepancias resultantes en los intervalos.
Para ver los grados de libertad para varios intervalos, se pueden agregar pruebas personalizadas en JMP. Este script incluye pruebas personalizadas equivalentes al intervalo de confianza medio y al IC medio condicional. Puede ver el valor predicho y SE coincide con lo que se obtiene al guardar esas columnas. La información adicional es el DFDen utilizado para la prueba / CI. Ese DFDen se puede utilizar para encontrar el valor crítico a utilizar con el SE para el CI.
Fit Model(
Y( :Weight ),
Effects( :Feed, :Time ),
Random Effects( :Pig ),
NoBounds( 1 ),
Personality( "Standard Least Squares" ),
Method( "REML" ),
Emphasis( "Minimal Report" ),
Run(
:Weight << {Summary of Fit( 1 ), Analysis of Variance( 0 ),
Parameter Estimates( 1 ), Scaled Estimates( 0 ),
Plot Actual by Predicted( 0 ), Plot Regression( 0 ),
Plot Residual by Predicted( 0 ), Plot Studentized Residuals( 0 ),
Plot Effect Leverage( 1 ), Plot Residual by Normal Quantiles( 0 ),
Custom Test(
[1 5.200005 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0],
Label( "obs2 marginal" )
), Custom Test(
[1 5.200005 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0],
Label( "obs2 conditional" )
)}
)
);
Supongo que esto 'explica' de dónde viene el problema, pero no explica por qué. No sé cómo R o Python obtienen los SE de los componentes de la varianza, pero claramente hay una diferencia allí. Los procedimientos de JMP y SAS están completamente documentados (creo que el documento de JMP apunta al documento de SAS para obtener toda la información de inferencia estadística profunda).Documento relevante para el procedimiento Viya LMIXED (igual que PROC MIXED regular)
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?Noté que JMP y SAS lo etiquetan como Error estándar, mientras que el paquete R llama Std.Dev. ¿Hay una diferencia? En su informe de efectos fijos, R se etiqueta como Std. Error.
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?Probablemente no. Aunque técnicamente no es lo mismo, por lo general, los términos se usan indistintamente.
Re: modelos mixtos / resultados de REML: ¿tenemos documentación sobre cómo se calculará el error estándar?¡Gracias Elisabeth!
En cuanto al algoritmo de optimización utilizado, en el documento se menciona a Newton-Raphson, pero con una opción a cambiar en LMIXED. Supongo que en JMP solo usamos NR y no tenemos un cambio para cambiar, ¿correcto?
Y tampoco habrá ninguna opción para escribir la matriz de Hesse. Teniendo esta salida, podríamos comparar si ha habido una diferencia de iteración a iteración que difiera del algoritmo R, aunque no debería haber una matemáticamente. El trabajo de comparación lo dejaría hacer al cliente).
¡Gracias nuevamente por sus investigaciones! Eso fue muy útil.
