Verificación de supuestos del modelo de datos.
La validez de los resultados obtenidos en cualquier análisis de varianza ANOVA queda supeditado a que los supuestos del modelo se cumplan. Estos supuestos son:
- Normalidad
- Varianza constante (igual varianza de los tratamientos) u homocedasticidad.
- Independencia.
Esto es, la respuesta \(y_{ij}\) se debe distribuir de manera normal, con la misma varianza en cada tratamiento y las mediciones deben ser independientes.
\[y_{ij}\sim N(\mu+\tau_i~,~\sigma^2)~, \\ Igual~varianza~\sigma^2~para~todos~tratamientos~i~, \\ Mediciones~independientes\]
Estos supuestos sobre \(y_{ij}\) se traducen en supuestos sobre el término del error aleatorio \(\epsilon_{ij}\) en el modelo \(y_{ij}=\mu+\tau_i+\epsilon_{ij}\), en la práçtica se realizan la verificación de los supuestos sobre los residuales
Definición de residuales \(e_{ij}\)
Es una práctica común utilizar la muestra de ** residuales o residuos** para comprobar los supuestos del modelo, ya que si los supuestos se cumplen, los residuos o residuales se pueden ver como una muestra aleatoria de una distribución normal con media cero y varianza constante, esto es:
\[e_{ij} \sim N(0,\sigma^2)\]
Los residuales \(e_{ij}\) se definen como la diferencia entre la observación \(y_{ij}\) y la respuesta predicha por el modelo de datos \(\hat{y_{ij}}\) por lo tanto:
\[e_{ij} = y_{ij}-\hat{y_{ij}}\]
El cálculo de residuales permite hacer un diagnóstico más directo de la calidad del modelo, ya que su magnitud señala qué tan bien describe a los datos el modelo estadístico.
Recordemos que el modelo que describe los datos en un diseño de experimentos de un factor es:
\[y_{ij}=\mu+\tau_i+\epsilon_{ij}\]
Cuando se realiza el ANOVA, y sólo cuando éste resulta significativo (rechazo \(H_o\)), entonces se procede a estimar el modelo ajustado o modelo de trabajo dado por:
\[\hat{y_{ij}}=\hat{\mu}+\hat{\tau_i}\]
El término del error \(\epsilon_{ij}\) desaparece de la expresión del modelo de datos puesto que su valor esperado \(E(\epsilon_{ij})=0\)
También es importante recordar que:
\[\hat{\mu}=\bar{y_{..}} \\ \tau_i = \bar{y_{i.}} - \bar{y_{..}}\]
Por lo que el modelo se puede escribir como:
\[\hat{y_{ij}}=\hat{\mu}+\hat{\tau_i} = \bar{y_{..}}-(\bar{y_{i.}}+\bar{y_{..}}) \rightarrow\]
\[\hat{y_{ij}} = \bar{y_{i.}}\]
Dado lo anterior, el residual \(e_{ij}\) asociado a la observación \(y_{ij}\) está dado por:
\[e_{ij}=y_{ij}-\hat{y_{ij}}= y_{ij} - \bar{y_{i.}}\]
Los supuestos del modelo, traducidos a los residuales \(e_{ij}\) son:
- Los \(e_{ij}\) siguen una distribución normal con media cero.
\[e_{ij}\sim N(0,\sigma^2)\]
Los \(e_{ij}\) de cada tratamiento tienen la misma varianza \(\sigma^2\). Homocedasticidad de la varianza de los residuales.
Los \(e_{ij}\) son independientes entre sí.
Para comprobar cada supuesto existen pruebas analíticas y gráficas que veremos a continuación.
Verificación de supuesto de normalidad de residuales \(e_{ij}\)
Verificación gráfica
Un procedimiento gráfico para verificar el supuesto de normalidad de los residuales consiste en graficar los residuos versus una probabilidad acumulada o versus la normal inversa de la probabilidad acumulada. Para ello se siguen los siguientes pasos:
Es necesario poseer los cálculos de los \(N\) residuales \(e_{ij}\).
Ordenar los \(N\) valores de los residuales en orden ascendente y asignarles en ese orden los valores de \(1~a~N\).Sean \(r_i\) con \(i=1,2,...,N\) los datos de los residuales \(e_{ij}\) en orden creciente.
Calcular una posición de graficación para cada dato en función de su rango y del total de observaciones como:
\[p_i=\frac{i-0.5}{N}~,~i=1,2,...,N\]
- Calcular la normal inversa \(Z_i\) de \(\frac{i-0.5}{N}~,~i=1,2,...,N\), es decir:
\[Z_i=\phi^{-1}\left(\frac{i-0.5}{N}\right)\]
- Graficar en el eje de las \(x\) los \(r_i\) y en el eje de las \(y\) la posición \(p_i\) o \(Z_i\)
Verificación gráfica ejemplo
Tomemos el ejemplo inicial de la evaluación de la resistencia a la tensión de una fibra sintética para distintos porcentajes de algodón en la fibra. Los residuales calculados son los siguiente:
| Residuales |
|---|
| -2.8 |
| -2.8 |
| 5.2 |
| 1.2 |
| -0.8 |
| -3.4 |
| 1.6 |
| -3.4 |
| 2.6 |
| 2.6 |
| -3.6 |
| 0.4 |
| 0.4 |
| 1.4 |
| 1.4 |
| -2.6 |
| 3.4 |
| 0.4 |
| -2.6 |
| 1.4 |
| -3.8 |
| -0.8 |
| 0.2 |
| 4.2 |
| 0.2 |
Organizamos los residuales de manera ascendente, calculamos la posición \(\frac{i-0.5}{N}\) y \(Z_i\)y obtenemos la siguiente información
| Posición \(i\) | Residuales \(r_i\) | \(\frac{i-0,5}{N}\) | \(Z_i\) |
|---|---|---|---|
| 1 | -3.8 | 0.02 | -2.0537489 |
| 2 | -3.6 | 0.06 | -1.5547736 |
| 3 | -3.4 | 0.10 | -1.2815516 |
| 4 | -3.4 | 0.14 | -1.0803193 |
| 5 | -2.8 | 0.18 | -0.9153651 |
| 6 | -2.8 | 0.22 | -0.7721932 |
| 7 | -2.6 | 0.26 | -0.6433454 |
| 8 | -2.6 | 0.30 | -0.5244005 |
| 9 | -0.8 | 0.34 | -0.4124631 |
| 10 | -0.8 | 0.38 | -0.3054808 |
| 11 | 0.2 | 0.42 | -0.2018935 |
| 12 | 0.2 | 0.46 | -0.1004337 |
| 13 | 0.4 | 0.50 | 0.0000000 |
| 14 | 0.4 | 0.54 | 0.1004337 |
| 15 | 0.4 | 0.58 | 0.2018935 |
| 16 | 1.2 | 0.62 | 0.3054808 |
| 17 | 1.4 | 0.66 | 0.4124631 |
| 18 | 1.4 | 0.70 | 0.5244005 |
| 19 | 1.4 | 0.74 | 0.6433454 |
| 20 | 1.6 | 0.78 | 0.7721932 |
| 21 | 2.6 | 0.82 | 0.9153651 |
| 22 | 2.6 | 0.86 | 1.0803193 |
| 23 | 3.4 | 0.90 | 1.2815516 |
| 24 | 4.2 | 0.94 | 1.5547736 |
| 25 | 5.2 | 0.98 | 2.0537489 |
Ahora podemos graficar los residuales \(e_{ij}\) organizados de manera ascendente \(r_i\) versus \(\frac{i-0.5}{N}\)
Podemos graficar los residuales \(e_{ij}\) organizados de manera ascendente \(r_i\) versus \(Z_i\)
Esta gráfica tiene las escalas de tal manera que si los residuos siguen una distribución normal, al graficarlos tienden a quedar alineados en una línea recta; por lo tanto, si claramente no se alinean se concluye que el supuesto de normalidad no es correcto.
Cabe enfatizar el hecho de que el ajuste de los puntos a una recta no tiene que ser perfecto, dado que el análisis de varianza resiste pequeñas y moderadas desviaciones al supuesto de normalidad.
Verificación analítica.
Prueba de Shapiro - Wilks
Consideremos una muestra aleatoria de datos \(x_1,x_2,...x_n\) que procede de cierta función desconocida denotada \(F(x)\). Se quiere verificar si dichos datos fueron generados por un proceso normal, mediante las hipótesis estadísticas.
- Planteamos las hipótesis
\[H_o: e_{ij}\sim N(\mu,\sigma^2)\\H_1: e_{ij} \nsim N(\mu,\sigma^2)\]
Organizamos los residuales \(e_{ij}\) de manera ascendente, denotemos los datos ordenados como: \(r_i= x_{(1)},x_{(2)},...,x_{(n)}\)
Calculamos estadístico de prueba \(W\) definido por:
\[W=\frac{1}{(n-1)S^2}\left[\sum_{i=1}^{h}{a_i(x_{n-i+1}-x_i)}\right]^2~ con~ h= \left\lbrace\begin{array}{c} \frac{n}{2}~si~n~es~par \\ \frac{n-1}{2}~si~n~es~impar\end{array}\right.\]
- Los coeficientes \(a_1,a_2,...,a_n\) se obtienen de tablas precalculadas.
- \(S^2\) es la varianza muestral de los \(x_i\).
- Comparamos estadístico teórico de tabla.
Si
\[W>W_{1-\alpha}~\rightarrow ~ Rechazo~H_o\]
\(\alpha\) es el nivel de significancia
Solución de situación inicial en R
La situación inicial es la siguiente:
Un ingeniero de desarrollo de productos tiene interés en investigar la resistencia a la tensión de una fibra sintética nueva que se usará para hacer tela de camisas para caballero. El ingeniero sabe por experiencia previa que la resistencia a la tensión se afecta por el peso porcentual del algodón utilizado en la mezcla de materiales de la fibra. Además, sospecha que al aumentar el contenido de algodón se incrementará la resistencia, al menos en un principio. Sabe asimismo que el contenido de algodón deberá variar entre \(10\%\) y \(40\%\) para que el producto final tenga otras características de calidad que se desean (como la capacidad de ser sometido a un tratamiento de planchado permanente). El ingeniero decide probar ejemplares en cinco niveles del peso porcentual del algodón: \(15\%\), \(20\%\), \(25\%\), \(30\%\) Y \(35\%\) por ciento. También decide probar cinco ejemplares en cada nivel del contenido de algodón. La siguiente tabla muestra los datos:
| Peso porcentual algodón | 1 | 2 | 3 | 4 | 5 | Total | Promedio |
|---|---|---|---|---|---|---|---|
| 15% | 7 | 7 | 15 | 11 | 9 | 49 | 9.80 |
| 20% | 12 | 17 | 12 | 18 | 18 | 77 | 15.40 |
| 25% | 14 | 18 | 18 | 19 | 19 | 88 | 17.60 |
| 30% | 19 | 25 | 22 | 19 | 23 | 108 | 21.60 |
| 35% | 7 | 10 | 11 | 15 | 11 | 54 | 10.80 |
| 376 | 15.04 |
La solución en R para verificar la normalidad sería la siguiente:
# Análisis de varianza
library(readxl)
datos <- read_excel("Datos_base_1.xlsx")
datos$Porc_Algodon<- as.factor(datos$Porc_Algodon)
modelo <- lm(datos$Resistencia~datos$Porc_Algodon)
anova <- aov(modelo)
# Residuales
residuales <- anova$residuals
# Test Shapiro Wilk
shapiro.test(residuales)##
## Shapiro-Wilk normality test
##
## data: residuales
## W = 0.94387, p-value = 0.1818De R obtenemos los siguientes resultados.
Estadístico de prueba: \(W= 0.94387\)
P-value: \(0.1818\)
De las tablas de cuantiles teóricos para el estadístico W, obtenemos el siguiente dato, si \(\alpha = 0.05\) y la cantidad de residuales es \(n=25\)
- \(W_{1-\alpha} = W_{1-0.05} = W_{0.95} = 0.985\)
Por lo que:
\[W= 0.9438 \ngtr W_{0.95} = 0.985 \]
Por lo que no se rechaza la hipótesis nula \(H_0\) y los residuales \(e_{ij}\) proceden de una distribución normal.
Solución problema de clase en R
Para el problema usado en clases el procedimiento de solución para el Test de Shapiro Wilk es el siguiente:
# Análisis de varianza
cuero <- c("A", "A","A","A","A","A", "B", "B","B","B","B","B","C","C","C","C","C","C","D","D","D","D","D","D")
desgaste <- c(264,260,258,241,262,255,208,220,216,200,213,206,220,263,219,225,230,228,217,226,215,227,220,222)
cuero <- as.factor(cuero)
modelo <- lm(desgaste~cuero)
anova <- aov(modelo)
# Residuales
residuales <- anova$residuals
# Test de Shapiro - Wilk
shapiro.test(residuales)##
## Shapiro-Wilk normality test
##
## data: residuales
## W = 0.88326, p-value = 0.00967De R obtenemos los siguientes resultados.
Estadístico de prueba: \(W= 0.88326\)
P-value: \(0.00967\)
De las tablas de cuantiles teóricos para el estadístico W, obtenemos el siguiente dato, si \(\alpha = 0.05\) y la cantidad de residuales es \(n=24\)
- \(W_{1-\alpha} = W_{1-0.05} = W_{0.95} = 0.984\)
Por lo que:
\[W= 0.88326 \ngtr W_{0.95} = 0.984 \]
Por lo que no se rechaza la hipótesis nula \(H_0\) y los residuales \(e_{ij}\) proceden de una distribución normal.
Prueba de Kolmogorov - Smirnov (K-S)
Consideremos una muestra aleatoria de datos \(x_1,x_2,...x_n\) que procede de cierta función desconocida denotada \(F(x)\). Se quiere verificar si dichos datos fueron generados por un proceso normal, mediante las hipótesis estadísticas.
- Planteamos las hipótesis
\[H_o: e_{ij}\sim N(\mu,\sigma^2)\\H_1: e_{ij} \nsim N(\mu,\sigma^2)\]
Organizamos los residuales \(e_{ij}\) de manera ascendente, denotemos los datos ordenados como: \(x_i=x_{(1)},x_{(2)},...,x_{(n)}~con~i=1,2,3...n\).
Calculamos probabilidad teórica (percentil) \(P_i\) de la forma:
\[P_i=i/n\]
- Estandarizamos los residuales organizados \(x_i\) de la siguiente manera:
\[Z_i=\frac{x_i-\bar{x}}{S}\]
- \(\bar{x}\) corresponde a la media muestra.
- \(S\) es la desviación estándar de la muestra.
Calculamos la probabilidad \(P(Z_i)\) para una distribución normal estándar.
Calculamos las distancias:
- Distancia 1.
\[D_1=|P(Z_i)-P_i|\]
- Distancia 2.
\[D_2=|P(Z_i)-P_{i-1}|\]
- Escojo el valor para el estadístico de prueba \(D\)
\[D=max\left\lbrace D_1,D_2 \right\rbrace\]
- Comparamos el estadístico de prueba \(D\) con el estadístico teórico \(KS\)
\[KS=\frac{C_\alpha}{k(n)}\]
\(C_\alpha\) viene dado por tablas de la siguiente manera
| \(C_\alpha\) por distribución | \(\alpha=0.1\) | \(\alpha=0.05\) | \(\alpha=0.01\) |
|---|---|---|---|
| Normal | 0.819 | 0.895 | 1.035 |
| Exponencial | 0.990 | 1.094 | 1.308 |
| Weibull \(n=10\) | 0.760 | 0.819 | 0.944 |
El valor de \(k(n)\) también se encuentra tabulado, y lo encontramos de la siguiente forma:
| \(K(n)\) por distribución | $K(n) |
|---|---|
| Normal | \(k(n)=\sqrt{n}-0,01+\frac{0,85}{\sqrt{n}}\) |
| Exponencial | \(k(n)=\sqrt{n}+0,12+\frac{0,11}{\sqrt{n}}\) |
| Weibull \(n=10\) | \(k(n)=\sqrt{n}\) |
Si
\[D>KS \rightarrow Rechazo~H_o\]
Solución de situación inicial en R
# Análisis de varianza
library(readxl)
datos <- read_excel("Datos_base_1.xlsx")
datos$Porc_Algodon<- as.factor(datos$Porc_Algodon)
modelo <- lm(datos$Resistencia~datos$Porc_Algodon)
anova <- aov(modelo)
# Residuales
residuales <- anova$residuals
# Test Kolmogorov Smirnov
library(nortest)
lillie.test(residuales)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuales
## D = 0.16212, p-value = 0.08889De los resultados de R obtenemos:
\(D = 0.16212\)
\(P-value = 0.08889\)
El cálculo del estadístico de referencia resulta en:
- \(KS =0.17345\)
Como
\[D = 0.16212 \ngtr KS= 0.17345 \]
No existe evidencia estadística para rechazar \(H_0\) por lo que los residuales provienen de una distribución normal.