Calculadora de Shapiro–Wilk (normalidad)

Explicación

El test de Shapiro–Wilk es una prueba estadística diseñada para evaluar si una muestra de datos puede considerarse proveniente de una distribución normal. Es uno de los tests de normalidad más potentes y recomendados para tamaños muestrales pequeños y moderados.

El test compara:

• Los datos ordenados de tu muestra
• l os valores esperados de una muestra normal ordenada

La hipótesis nula plantea que los datos siguen una normal; la alternativa plantea que no siguen una normal. El estadístico del test se denota por \(W\), y toma valores entre 0 y 1: cuanto más cerca de 1, mayor compatibilidad con normalidad.

El estadístico de Shapiro–Wilk se construye combinando los datos ordenados \(x_{(i)}\) con coeficientes \(a_i\) que dependen de \(n\):

\(W = \dfrac{\left(\sum_{i=1}^{m} a_i(x_{(n+1-i)}-x_{(i)})\right)^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}\) donde:

• \(x_{(i)}\) son los datos ordenados,
• \(a_i\) son coeficientes calculados a partir de los cuantiles normales teóricos,
• el denominador es la variabilidad total de la muestra.

Simulación de Monte Carlo para el p-valor

La distribución exacta del estadístico \(W\) bajo \(H_0\) no tiene una forma cerrada sencilla. Existen aproximaciones analíticas (la de Royston, que es la que usa R), pero esta calculadora estima el p-valor mediante simulación de Monte Carlo. El procedimiento es el siguiente:

1. Se genera una muestra aleatoria de tamaño \(n\) a partir de \(N(0,1)\). Puede usarse la normal estándar porque \(W\) es invariante a cambios de media y escala: solo depende de la forma de los datos ordenados, no de su nivel o dispersión.
2. Se calcula el estadístico \(W\) en esa muestra simulada.
3. Se repite \(B\) veces (por ejemplo, \(B = 5000\)). Con ello se obtiene una distribución empírica de \(W\) bajo \(H_0\): es decir, cómo se comportaría \(W\) si los datos realmente fueran normales.
4. Se cuenta el número \(b\) de simulaciones al menos tan extremas como la observada. En Shapiro–Wilk son los valores pequeños de \(W\) los que indican falta de normalidad, así que \(b = \#\{W_{\text{sim}} \leq W_{\text{obs}}\}\), y el p-valor Monte Carlo se calcula como \[\hat p = \dfrac{b+1}{B+1}\] Un \(W\) observado muy pequeño (lejos de 1) caerá en la cola izquierda de la distribución simulada, dando un p-valor pequeño y llevando a rechazar la normalidad.

¿Por qué \((b+1)/(B+1)\) y no \(b/B\)?

La fracción \(b/B\) es el estimador natural de la probabilidad de cola \(p = P(W \leq W_{\text{obs}} \mid H_0)\) (de hecho es insesgado para ella), pero cuando el p-valor se usa para decidir con un número finito de réplicas conviene la corrección \((b+1)/(B+1)\), recomendada en la literatura de pruebas Monte Carlo (Davison y Hinkley, 1997; North, Curtis y Sham, 2002; Phipson y Smyth, 2010). La justificación es la siguiente:

• Validez exacta del contraste. Bajo \(H_0\), el \(W\) observado y los \(B\) simulados son \(B+1\) valores intercambiables de la misma distribución, de modo que la posición del observado entre ellos es uniforme en \(\{1, \dots, B+1\}\). Por eso \(P\bigl(\hat p \leq \alpha\bigr) \leq \alpha\) para cualquier \(B\), con igualdad cuando \(\alpha(B+1)\) es entero: la tasa de error de tipo I nunca supera el \(\alpha\) nominal. Con \(b/B\), en cambio, \(P(b/B \leq \alpha) = (\lfloor \alpha B \rfloor + 1)/(B+1)\), que puede superar ligeramente \(\alpha\) (con \(B = 5000\) y \(\alpha = 0{,}05\) es \(251/5001 \approx 0{,}0502\)).
• Nunca da \(p = 0\). Si ninguna réplica es tan extrema como la muestra, \(b/B\) devolvería un p-valor nulo, que sugiere una certeza imposible con \(B\) réplicas. La versión corregida devuelve \(1/(B+1)\), que es la menor probabilidad que una simulación de tamaño \(B\) puede resolver.
• El coste es despreciable. La diferencia entre ambos estimadores es como mucho \(1/(B+1)\) (≈ 0,0002 con \(B = 5000\)), muy por debajo del error de simulación salvo en la cola extrema, así que solo afecta a casos límite, y en ellos lo hace en la dirección conservadora.

El histograma del gráfico muestra precisamente esa distribución empírica: las barras acumulan los \(B\) valores simulados de \(W\). El marcador rojo señala el \(W\) observado en tu muestra; si cae en la cola izquierda (zona de valores bajos), hay evidencia contra la normalidad. El marcador verde indica el valor crítico \(W_\alpha\): aproximadamente el percentil \(\alpha\) de la distribución simulada, elegido de forma que \(W_{\text{obs}} < W_\alpha\) equivale exactamente a \(\hat p < \alpha\), de modo que el gráfico y la decisión siempre coinciden.

Incertidumbre Monte Carlo

La semilla aleatoria permite reproducir exactamente el mismo resultado: fijar la misma semilla garantiza que las \(B\) muestras simuladas sean idénticas en distintas ejecuciones, lo que facilita compartir y verificar resultados. Pero el p-valor Monte Carlo sigue siendo una estimación del p-valor que se obtendría con infinitas réplicas, \(p = P(W \leq W_{\text{obs}} \mid H_0)\). Como el mismo cálculo de \(W\) se aplica a la muestra observada y a las simuladas, el error de simulación es la única diferencia entre ambos.

El número de réplicas tan extremas sigue una distribución binomial, \(b \sim \text{Bin}(B, p)\), así que el error estándar del p-valor estimado es aproximadamente \(\sqrt{p(1-p)/B}\): depende tanto del propio p-valor como del número de réplicas. Con \(B = 5000\), cerca de \(p = 0{,}05\) es \(\approx 0{,}0031\) (un margen aproximado del 95 % de \(\pm 0{,}006\)), y cerca de \(p = 0{,}5\) es \(\approx 0{,}007\). Además del error estándar, la calculadora muestra un intervalo de confianza del 95 % de Wilson para \(p\) a partir de \(b\) y \(B\), que sigue siendo informativo cuando \(b\) es pequeño o nulo (con \(b = 0\) su límite superior es \(\approx 3{,}84/B\)). Si \(\alpha\) cae dentro de ese intervalo, la simulación no basta para saber de qué lado de \(\alpha\) está el p-valor y la calculadora lo advierte.

Esa precisión es aceptable para muchos usos, pero no es universalmente suficiente. El error disminuye como \(1/\sqrt{B}\) (cuadruplicar \(B\) lo reduce a la mitad) y, para lograr un error estándar \(\delta\), hacen falta unas \(B \approx p(1-p)/\delta^2\) réplicas: por ejemplo, unas 47 500 para \(\delta = 0{,}001\) cerca de \(p = 0{,}05\). Si el p-valor queda cerca de \(\alpha\), o si se necesita estimar con precisión p-valores muy pequeños, conviene aumentar \(B\).

La decisión formal se realiza con el p-valor: si es menor que el nivel de significación \(\alpha\), se rechaza la normalidad.

Hipótesis y estadístico

\(H_0\): la muestra procede de una distribución normal

\(H_1\): la muestra no procede de una distribución normal

\(W = \dfrac{\left(\sum_{i=1}^{m} a_i(x_{(n+1-i)}-x_{(i)})\right)^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}\)

Contraste rápido

En Shapiro–Wilk, el p-valor es la base de la decisión estadística. El estadístico \(W\) resume cuánto se desvían los datos ordenados respecto al patrón esperado bajo normalidad.

\(\text{Si } p < \alpha\Rightarrow\) rechazar \(H_0\)

\(\text{Si } p \ge \alpha\Rightarrow\) no rechazar \(H_0\)

Un resultado no significativo no demuestra normalidad perfecta: indica que, con la evidencia disponible en esa muestra, no se detectan desviaciones claras. Conviene complementar con gráfico Q–Q e histograma.

Ejemplo resuelto

En un proceso de control de calidad se miden las dimensiones (en mm) de 10 piezas fabricadas: 23,1; 24,5; 22,8; 25,2; 24,1; 23,7; 24,9; 23,5; 25,0; 24,3. Se desea verificar, con \(\alpha = 0{,}05\), si estos datos son compatibles con una distribución normal.

Paso 1 – Ordenar los datos, calcular la media y la suma de cuadrados

Datos ordenados: 22,8; 23,1; 23,5; 23,7; 24,1; 24,3; 24,5; 24,9; 25,0; 25,2.

\(\bar{x} = \dfrac{22{,}8 + 23{,}1 + \cdots + 25{,}2}{10} = \dfrac{241{,}1}{10} = 24{,}11 \text{ mm}\)

\(\displaystyle SS = \sum_{i=1}^{10}(x_i - \bar{x})^2 = (22{,}8-24{,}11)^2 + \cdots + (25{,}2-24{,}11)^2 = 6{,}069\)

Paso 2 – Calcular los coeficientes \(a_i\)

La calculadora usa el algoritmo de Royston (1992; AS R94, 1995), que parte de los cuantiles normales esperados y aplica correcciones polinómicas a los coeficientes extremos para aproximar los coeficientes de Shapiro–Wilk. No son los coeficientes exactos, cuyo cálculo requiere la matriz de covarianzas de los estadísticos de orden de la normal. Para \(n = 10\), los coeficientes de Royston y los tabulados por Shapiro y Wilk (1965) son:

La aproximación de Royston difiere de los valores tabulados en la cuarta cifra decimal y es la misma que usa R. En los pasos siguientes se usan los valores tabulados, habituales en el cálculo manual; con los de Royston, \(b\) y \(W\) coinciden en sus cuatro primeros decimales.

Paso 3 – Calcular el numerador

Cada \(a_i\) pondera la diferencia entre el \(i\)-ésimo dato por arriba y por abajo de la muestra ordenada:

\(b = \sum_{i=1}^{5} a_i\bigl(x_{(11-i)} - x_{(i)}\bigr)\)

\(b = 0{,}5739\cdot(25{,}2-22{,}8) + 0{,}3291\cdot(25{,}0-23{,}1) + 0{,}2141\cdot(24{,}9-23{,}5) + 0{,}1224\cdot(24{,}5-23{,}7) + 0{,}0399\cdot(24{,}3-24{,}1)\)

\(b = 1{,}3774 + 0{,}6253 + 0{,}2997 + 0{,}0979 + 0{,}0080 = 2{,}4083\)

Paso 4 – Estadístico \(W\)

\(W = \dfrac{b^2}{SS} = \dfrac{2{,}4083^2}{6{,}069} = \dfrac{5{,}800}{6{,}069} \approx 0{,}9557\)

Este valor coincide con el que devuelve R: shapiro.test(...)$statistic = 0.9557.

Paso 5 – p-valor mediante simulación Monte Carlo

La calculadora genera \(B = 5000\) muestras normales \(N(0,1)\) de tamaño \(n = 10\), calcula \(W\) en cada una y cuenta las simulaciones con \(W_{\text{sim}} \leq W_{\text{obs}}\). Con semilla 12345 son \(b = 3663\), así que:

\(\hat p = \dfrac{b+1}{B+1} = \dfrac{3664}{5001} \approx 0{,}7327\)

(La fracción sin corregir, \(3663/5000 = 0{,}7326\), es prácticamente igual: con un p-valor lejos de \(\alpha\) la corrección no cambia nada.) El error estándar de Monte Carlo es \(\sqrt{0{,}7327 \cdot 0{,}2673/5000} \approx 0{,}006\) y el intervalo de Wilson del 95 % para el p-valor es aproximadamente \([0{,}720;\ 0{,}745]\). R obtiene \(p = 0{,}7355\) con la aproximación analítica de Royston (una transformación normalizadora de \(W\)). Ambos son aproximaciones del mismo p-valor y resultan compatibles dentro del error de simulación (0,7355 está dentro del intervalo), pero no cabe esperar que coincidan en todos los decimales. Puedes reproducirlo introduciendo los datos en la calculadora con \(B = 5000\) y semilla 12345.

Conclusión

Dado que \(p \approx 0{,}73 \gg 0{,}05 = \alpha\), no se rechaza \(H_0\). Con un nivel de significación del 5 %, no hay evidencia estadística de que las dimensiones de las piezas no sigan una distribución normal. Esto no demuestra que lo sean: con \(n = 10\) la potencia del test es limitada y solo detectaría desviaciones marcadas. Antes de apoyarse en herramientas paramétricas como los gráficos de control \(\bar{x}\text{-}R\), conviene revisar también un gráfico Q–Q y el conocimiento previo del proceso.

Cómo interpretar el resultado

Rechazar \(H_0\) (p-valor < \(\alpha\)) indica que hay evidencia estadística de que los datos no provienen de una distribución normal. Sin embargo, la causa puede ser diversa: colas más pesadas de lo normal (leptocurtosis), asimetría, o simplemente la presencia de valores atípicos. Para identificar el tipo de desviación, complementa con el gráfico Q-Q (cuantil-cuantil): si los puntos se curvan hacia arriba en la cola derecha hay asimetría positiva; si las colas se separan de la línea de referencia, hay exceso de curtosis.

No rechazar \(H_0\) (p-valor ≥ \(\alpha\)) no demuestra normalidad; solo indica que los datos son compatibles con una distribución normal al nivel elegido. Con muestras pequeñas (\(n < 20\)), el test tiene poca potencia frente a desviaciones moderadas (colas algo más pesadas, asimetría leve) y con frecuencia no rechazará \(H_0\) aunque los datos no sean normales; solo las desviaciones marcadas se detectan con regularidad. Con muestras grandes (\(n > 200\)), en cambio, desviaciones triviales pueden resultar altamente significativas sin impactar en la validez de los métodos paramétricos que asumen normalidad. En la práctica, lo relevante es si la desviación es suficientemente grave para invalidar el test que se pretende aplicar.

El estadístico W de Shapiro-Wilk varía entre 0 y 1: valores próximos a 1 indican que la muestra es consistente con la normalidad; valores cercanos a 0 señalan una desviación severa. A diferencia de los tests basados en chi-cuadrado o KS, Shapiro-Wilk tiene mayor potencia para detectar desviaciones en muestras pequeñas y es el test de normalidad recomendado para \(n \leq 50\). La herramienta muestra el valor W, el p-valor Monte Carlo con su error estándar y su intervalo de confianza Monte Carlo y, en el gráfico, la distribución simulada de W bajo normalidad junto con el W observado y el valor crítico. Para ver cómo se desvían los datos de la normal, complementa el resultado con un gráfico Q–Q.

Preguntas frecuentes

• ¿Qué tamaños muestrales son habituales para Shapiro–Wilk? Se usa mucho en muestras pequeñas y medianas, aunque también puede aplicarse en tamaños mayores según software.
• ¿Un p-valor alto demuestra normalidad? No. Solo indica que no hay evidencia suficiente para rechazar normalidad con esa muestra y ese nivel de significación.
• ¿Debo usar solo el test? No. La práctica recomendada es combinar test formal con gráfico Q–Q e histograma.

Referencias:

• Shapiro, S. S. & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3–4), 591–611. — Paper original del test y de los coeficientes \(a_i\) tabulados.
• Royston, J. P. (1992). Approximating the Shapiro–Wilk W-test for non-normality. Statistics and Computing, 2, 117–119. — Algoritmo de aproximación polinómica de los coeficientes \(a_i\) implementado en esta calculadora.
• Royston, J. P. (1995). Remark AS R94: A remark on Algorithm AS 181: The W-test for normality. Applied Statistics, 44(4), 547–551. — Versión revisada del algoritmo (AS R94), base del cálculo de \(W\) en R y en esta herramienta.
• Acklam, P. J. An algorithm for computing the inverse normal cumulative distribution function (algoritmo publicado en la web del autor). — Fuente de los coeficientes de la aproximación racional de \(\Phi^{-1}(p)\) utilizada en el código, con un error relativo inferior a \(1{,}2 \times 10^{-9}\).