¿Por qué los valores p cambian de importancia cuando se cambia el orden de las covariables en el modelo aov?

10

Tengo un conjunto de datos de 482 observaciones.

data=Populationfull

Voy a hacer un análisis de asociación de genotipo para 3 SNP. Estoy tratando de construir un modelo para mi análisis y estoy usando el aov (y ~ x, data = ...). Para un rasgo tengo varios efectos fijos y covariables que he incluido en el modelo, así:

Starts <- aov(Starts~Sex+DMRT3+Birthyear+Country+Earnings+Voltsec+Autosec, data=Populationfull)

summary(Starts)
                Df Sum Sq Mean Sq F value   Pr(>F)    
Sex              3  17.90    5.97  42.844  < 2e-16 ***
DMRT3            2   1.14    0.57   4.110    0.017 *  
Birthyear        9   5.59    0.62   4.461 1.26e-05 ***
Country          1  11.28   11.28  81.005  < 2e-16 ***
Earnings         1 109.01  109.01 782.838  < 2e-16 ***
Voltsec          1  12.27   12.27  88.086  < 2e-16 ***
Autosec          1   8.97    8.97  64.443 8.27e-15 ***
Residuals      463  64.48    0.14                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Descubrí que si cambiaba el orden de las variables en el modelo, obtenía diferentes valores de p, por favor vea a continuación.

Starts2 <- aov(Starts~Voltsec+Autosec+Sex+DMRT3+Birthyear+Country+Earnings, data=Populationfull)

summary(Starts2)

                Df Sum Sq Mean Sq F value   Pr(>F)    
Voltsec   1   2.18    2.18  15.627 8.92e-05 ***
Autosec   1 100.60  100.60 722.443  < 2e-16 ***
Sex              3  10.43    3.48  24.962 5.50e-15 ***
DMRT3            2   0.82    0.41   2.957  0.05294 .  
Birthyear        9   3.25    0.36   2.591  0.00638 ** 
Country          1   2.25    2.25  16.183 6.72e-05 ***
Earnings      1  46.64   46.64 334.903  < 2e-16 ***
Residuals      463  64.48    0.14                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

¿Por qué obtengo diferentes valores de p según el orden en que se codifican las variables / factores / covariables / efectos fijos (?)? ¿Hay alguna forma de "corregirlo"? ¿Puede ser que estoy usando el modelo equivocado? Todavía soy bastante nuevo en R, así que si puedes ayudarme con esto, hazlo realmente simple para que pueda entender la respuesta jeje ... ¡Gracias, espero que alguien pueda ayudarme a entender esto!

r anova Rbeginner
fuente
3
Proporcione algunos datos de muestra para Populationfullque su problema sea reproducible . Esto no sucede con el ejemplo de la aov()página de ayuda. summary(aov(yield ~ block + N + P + K, npk)); summary(aov(yield ~ K + P + block + N , npk))
MrFlick
Los valores de p están cambiando porque todo el campo de valores ha cambiado. tu primera carrera Earnings 1 109.01 109.01 782.838 < 2e-16 ***tu segunda carrera Earnings 1 46.64 46.64 334.903 < 2e-16 ***. Tus resultados no son los mismos. Comience verificando que no haya hecho más que reordenar variables.
1
TAMBIÉN: en el segundo modelo, usa Ganar, no Ganancias ... si hay dos variables de nombres diferentes, podría ser un problema si eso no es un error tipográfico al copiar en el espacio de preguntas SO.
Sí, los valores cambian pero ¿por qué? He usado exactamente las mismas columnas del mismo marco de datos exacto en ambos modelos (lo de Ganancias vs Ganancias en el segundo modelo es solo que lo escribí mal, lo he corregido ahora).
Rbeginner
1
Esto está sucediendo porque tienes un diseño desequilibrado. Encontrará mucha ayuda sobre esto si busca en este foro o simplemente "ANOVA desequilibrado en Google". Recomiendo buscar en el carpaquete: implementa ANOVA Tipo II y Tipo III, que no dependen del orden de las variables, mientras que lo aovhace ANOVA Tipo I.
Slow loris

Respuestas:

15

El problema proviene de la forma en que aov()realiza sus pruebas de importancia predeterminadas. Utiliza lo que se llama análisis ANOVA "Tipo I", en el que las pruebas se realizan en el orden en que las variables se especifican en su modelo. Entonces, en el primer ejemplo, determina cuánta varianza se explica sexy prueba su importancia, luego qué porción de la varianza restante se explica DMRT3y prueba su importancia en términos de esa varianza restante , y así sucesivamente. En el segundo ejemplo, DMRT3solo se evalúa después Voltsec, Autosecy sex, en ese orden, por lo que queda menos varianza para DMRT3explicar.

Si dos variables predictoras están correlacionadas, entonces la primera ingresada en el modelo obtendrá un "crédito" completo, dejando menos variación para ser "explicada por" la segunda, que por lo tanto puede parecer menos "estadísticamente significativa" que la primera, incluso si es no, funcionalmente. Esta pregunta y su respuesta explican los diferentes tipos de análisis ANOVA.

Una forma de evitar esto es extraerse de las restricciones del ANOVA clásico y utilizar un modelo lineal simple, con lm()R, en lugar de aov(). Esto analiza efectivamente todos los predictores en paralelo, "corrigiendo" todos los predictores a la vez. En ese caso, dos predictores correlacionados pueden terminar teniendo grandes errores estándar de sus coeficientes de regresión estimados, y sus coeficientes pueden diferir entre diferentes muestras de la población, pero al menos el orden en que ingrese las variables en la especificación del modelo no importará.

Si su variable de respuesta es algún tipo de variable de conteo, como su nombre lo Startsindica, entonces probablemente no debería usar ANOVA de todos modos ya que es poco probable que los residuos se distribuyan normalmente, como lo requiere la interpretación del valor p . Las variables de conteo se manejan mejor con modelos lineales generalizados (p. Ej., glm()En R), que pueden considerarse como una generalización de lm()otros tipos de estructuras de error residual.

EdM
fuente