Análisis Multivariante - Máster en técnicas estadísticas Curso 2010-2011

Tema 3. Contraste de la normalidad multivariante César A. Sánchez Sellero

Índice 1. Introducción

2

2. Contraste de la normalidad univariante

2

2.1. Test basado en la asimetría . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

2.2. Test basado en la kurtosis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.3. Test de Shapiro-Wilk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

3. Contraste de la normalidad multivariante

10

3.1. Conceptos previos: invariancia, estandarización, y distancias y ángulos de Mahalanobis 10 3.2. Tests basados en medidas de asimetría y kurtosis . . . . . . . . . . . . . . . . . . .

11

3.3. Extensiones multivariantes del test de Shapiro-Wilk . . . . . . . . . . . . . . . . . . .

15

3.4. Aplicación práctica de los test de normalidad multivariante . . . . . . . . . . . . . . .

15

1

2

Análisis Multivariante

1. Introducción En el tema anterior se han estudiado los procedimientos naturales de inferencia sobre el vector de medias y la matriz de covarianzas, bajo suposición de normalidad en la distribución del vector de variables observadas. Esta hipótesis de normalidad es esencial en muchos de aquellos procedimientos. Por ejemplo, la forma elíptica de las regiones de confianza es consustancial a la forma de la distribución normal multivariante. Si la forma en que se presentaran las observaciones en el espacio no fuera elíptica, tampoco lo debería ser una región de confianza para el vector de medias. En este tema veremos una breve introducción a los métodos de contraste de normalidad en observaciones multivariantes. Empezaremos con un repaso de los métodos de contraste de la normalidad en el caso univariante. Este será el objeto de la Sección 2, en la cual se podrán fijar algunos conceptos y prestar atención singular a los métodos susceptibles de ser extendidos al caso multivariante. En la Sección 3 se abordan algunos métodos (los más representativos) para el contraste de normalidad multivariante. Están basados en la extensión de las medidas de asimetría, kurtosis, así como del test de Shapiro-Wilk al caso multivariante. En este tema, además del objetivo de contraste de normalidad, conseguiremos una mejor comprensión de cómo se pueden disponer los datos en el espacio multidimensional.

2. Contraste de la normalidad univariante Se han desarrollado muchos procedimientos para el contraste de la normalidad en el caso univariante. No es el objetivo de esta sección el recorrerlos todos. Por el contrario, nos centraremos únicamente en aquellos que se extienden de manera más natural al caso multivariante. Serán los métodos basados en la asimetría, la kurtosis, y el test de Shapiro-Wilk.

2.1. Test basado en la asimetría Es bien conocido que la distribución normal es simétrica en torno a su media. Como además, la forma más común de desviarse respecto de la normalidad es por falta de simetría, por ejemplo en variables positivas, parece lógico construir un método de contraste en base a cierta medida de asimetría. Si X1 , . . . , Xn es una muestra aleatoria simple (observaciones independientes e idénticamente

Tema 3. Contraste de la normalidad multivariante

3

3 0

1

2

Frequency

4

5

6

Histogram of z

−1

0

1

2

3

z

Figura 1: Representación de 30 datos simulados con distribución normal estándar.

distribuidas), el coeficiente de asimetría muestral se define como

A= ¯ = (1/n) siendo X

∑n

i=1 Xi

(1/n)

∑n ( i=1

¯ Xi − X

S3

)3

1∑ = n n

i=1

la media muestral y S 2 = (1/n)

(

∑n

¯ Xi − X S

i=1 (Xi

)3

¯ 2 la varianza muestral. − X)

Es la media de las potencias cúbicas (también llamado momento de orden tres) de los valores estandarizados. Al estar estandarizados, habrá observaciones a ambos lados de la media, manteniendo un equilibrio (en orden uno). Sin embargo, la potencia tres que se emplea en el coeficiente de asimetría, rompe el equilibrio en caso de distribución asimétrica. Lo vemos en las Figuras 1 y 2. A pesar de haber sido simulados con distribución normal, en la Figura 1 hay una leve asimetría hacia la derecha, lo cual da lugar a un índice de asimetría muestral de 0’3384. En la Figura 2 se representa a la izquierda el histograma asociado a 30 datos simulados con distribución exponencial de media 1. En la derecha se muestra el histograma de los mismos datos, tras cambiarles el signo. El coeficiente de asimetría muestral vale 1’718 a la izquierda y -1’718 a la derecha. Si la distribución de los datos muestrales es normal, entonces el coeficiente de asimetría tiene distribución asintótica normal de media cero y varianza 6/n, por lo que se puede emplear como estadístico de contraste el siguiente:

√

n A ≈ N (0, 1) 6

4

Análisis Multivariante

8 10 6 4 0

0

2

4

6

Frequency

8 10

Histogram of x

2

Frequency

Histogram of x

0

1

2

3

4

5

6

7

x

−7 −6 −5 −4 −3 −2 −1

0

x

Figura 2: Representación de 30 datos simulados con distribución exponencial de media 1 (gráfico de la izquierda) y sus opuestos (gráfico de la derecha).

Se rechazará la normalidad, en base al coeficiente de asimetría, cuando el estadístico anterior sea muy grande (en positivo o negativo), en comparación con los cuantiles de la N (0, 1).

Ejemplo 1. Vamos a presentar los resultados del contraste sobre los mismos datos simulados que se han utilizado para ilustrar las explicaciones anteriores. El código R se encuentra en el anexo al final del tema, bajo el título de Ejemplo 1. A continuación se muestran los valores de los coeficientes de asimetría, los estadísticos de contraste y los niveles críticos asociados, para los datos simulados de la distribución normal, de la exponencial y los opuestos de los datos exponenciales. Como cabía esperar, la pequeña asimetría muestral de los datos normales no resulta significativa, mientras que sí lo es con los datos exponenciales.

> c(a_normal,a_exponencial,-a_exponencial) [1] 0.338375 1.717739 -1.717739 > c(estadistico_normal,estadistico_exponencial,-estadistico_exponencial) [1] 0.7566295 3.8409814 -3.8409814

Tema 3. Contraste de la normalidad multivariante

5

> c(pvalor_normal,pvalor_exponencial,pvalor_exponencial) [1] 0.4492718469 0.0001225434 0.0001225434

2.2. Test basado en la kurtosis Partiendo de la muestra X1 , . . . , Xn , el coeficiente de apuntamiento o kurtosis se define como

1∑ K= n n

i=1

∑n

(

¯ Xi − X S

)4 −3 ∑n

¯ 2 la varianza muestral. − X) En esta definición se resta tres para que el coeficiente valga cero en la distribución normal. ¯ = (1/n) siendo X

i=1 Xi

la media muestral y S 2 = (1/n)

i=1 (Xi

Recordemos que las características de posición se calculan mediante un momento de orden uno (la media), las de dispersión a través de un momento de orden dos (la varianza) y la asimetría con un momento de orden tres (el coeficiente de asimetría). Pues bien, la kurtosis se calcula a través del momento de orden cuatro, pues de esta manera se pueden medir características de forma, que son un refinamiento todavía más detallado en el análisis de la distribución. Lo vemos en la Figura 3, en la cual se representan tres distribuciones con la misma media, desviación típica y asimetría (asimetría cero, pues son simétricas). Son la T de Student con 5 grados de libertad, la normal y la uniforme, ambas con la misma media y desviación típica que la T5 . Para ello, se toma la normal con media cero y desviación típica

√

√ √ 5/3, y la uniforme entre − 5 y 5.

En la Figura 3 vemos que la distribución T de Student es la que presenta una mayor concentración de probabilidad en las proximidades de la media, así como en las colas (lejanías extremas)1 . Esto hace que tenga un coeficiente de apuntamiento más grande. Por el contrario, la distribución uniforme es la que presenta menor densidad en las proximidades de la media, y carece de colas. En este caso, la densidad es más grande en los valores con desviaciones intermedias (µ − σ) y (µ + σ). Como resultado, el coeficiente de apuntamiento será más pequeño para la distribución uniforme. La distribución normal presenta un comportamiento intermedio entre la T5 y la uniforme. 1

El fenómeno de la cola es difícilmente perceptible en la Figura 3. Aún así, se puede observar que el trazo rojo, correspondiente a la T5 , se encuentra ligeramente más elevado que el azul, correspondiente a la normal, en las colas. Es una pequeña diferencia a la vista, porque ambas densidades convergen a cero en las colas. Sin embargo, la consecuencia es que la T5 presenta una frecuencia no despreciable de observaciones muy grandes, cosa que no ocurre con la normal. Estas observaciones grandes (o extremas), al elevarlas a la potencia cuarta, aportan una contribución considerable a la kurtosis.

Análisis Multivariante

0.2 0.0

0.1

yt

0.3

6

−6

−4

−2

0

2

4

6

Figura 3: Funciones de densidad normal (trazo azul), T de Student (trazo rojo) y uniforme (trazo verde). La T de Student tiene cinco grados de libertad, y las distribuciones normal y uniforme han sido tomadas con la misma media y varianza que la T5 .

Si una distribución tiene un apuntamiento similar al de la normal decimos que es mesocúrtica, si el apuntamiento es mayor se dice leptocúrtica, y si es menor se dice platicúrtica. La T de Student es una distribución leptocúrtica y la uniforme platicúrtica. Volviendo al propósito de contraste de la normalidad, si estamos ante una distribución simétrica, podemos valorar una posible desviación del modelo normal en función de su apuntamiento o kurtosis. Rechazaremos la normalidad si la kurtosis muestral es mucho mayor o mucho menor que cero, que es lo que corresponde a la distribución normal. Ciertos cálculos de probabilidad, que omitimos aquí, permitieron determinar la distribución de la kurtosis muestral, bajo la hipótesis nula de normalidad. Así, si la distribución de los datos muestrales es normal, entonces la kurtosis tiene distribución asintótica normal de media cero y varianza 24/n, por lo que se puede emplear como estadístico de contraste el siguiente:

√

n K ≈ N (0, 1) 24

Se rechazará la normalidad, en base a la kurtosis, cuando el estadístico anterior sea muy grande (en positivo o negativo), en comparación con los cuantiles de la N (0, 1).

Ejemplo 2. Vamos a extraer muestras simuladas de tamaño 100 de cada una de las tres distribuciones que se han comentado en esta sección, y cuyas densidades se representan en la Figura 3:

Tema 3. Contraste de la normalidad multivariante

7

T5 de Student, normal y uniforme. El código R figura en el anexo al final de tema. A continuación mostramos los valores obtenidos para el coeficiente de apuntamiento, el estadístico de contraste, que no es más que la estandarización del coeficiente, y el nivel crítico para el contraste. Los resultados son los que cabía esperar. Así, el coeficiente de apuntamiento muestral es próximo a cero para la normal, negativo para la uniforme y positivo para la T5 . Al estandarizarlos, se obtienen cantidades comparables con la N (0, 1), que arrojan los niveles críticos que aparecen en la última línea. Las discrepancias más significativas se obtienen con la T5 de Student.

> c(k_normal,k_uniforme,k_tstudent) [1] -0.1864935 -1.1646089 1.6039029 > c(estadistico_normal,estadistico_uniforme,estadistico_tstudent) [1] -0.3806784 -2.3772480 3.2739530 > c(pvalor_normal,pvalor_uniforme,pvalor_tstudent) [1] 0.703441930 0.017442358 0.001060542

2.3. Test de Shapiro-Wilk Consideremos los datos estandarizados

Zi =

¯ Xi − X S

i ∈ {1, . . . , n}

El estadístico de Shapiro-Wilk se construye de la siguiente manera:

∑

[n/2]

W =

( ) ai,n Z(n−i+1):n − Zi:n

i=1

siendo Z1:n < . . . < Zn:n la muestra ordenada de los datos estandarizados y ai,n ciertas constantes. Consiste en calcular las distancias entre los datos de la muestra ordenada, simétricos respecto de la mediana, esto es, la distancia entre el primero y el último, el segundo y el penúltimo, y así sucesivamente; en general el Zi:n y el Z(n−i+1):n . El propósito es comparar estas distancias con las que habría en una muestra de observaciones normales. Al final, el coeficiente de determinación tiene una relación muy estrecha con el QQ-Plot, pues se puede interpretar como el coeficiente de determinación de los puntos representados en el QQ-Plot.

8

Análisis Multivariante

En consecuencia, se rechazará la normalidad cuando el estadístico de Shapiro-Wilk sea pequeño, pues eso será debido a que los puntos del QQ-Plot se alejan de la recta esperada.

Ejemplo 3. Vamos a extraer muestras simuladas de tamaño 100 de cada una de las distribuciones vistas en las secciones anteriores: normal, exponencial, uniforme y T de Student. Ciertos parámetros de estas distribuciones, como media, desviación típica u otros, que en las secciones anteriores se determinaron con el propósito de facilitar la comparación, ahora son irrelevantes, porque el test de Shapiro-Wilk comienza con una estandarización de los datos. La idea es que el test de Shapiro-Wilk es más general, pues permite detectar tanto defectos de simetría como de kurtosis, o incluso de otro tipo no considerado por los tests anteriores. En la Figura 4 se representan los cuatro QQ-Plots. El primero de ellos, situado en la esquina superior izquierda, corresponde a los datos simulados de una distribución normal. Hay leves desviaciones respecto de la recta en la cual estarían los datos que coinciden con los cuantiles esperados bajo normalidad. Estas desviaciones son perfectamente justificables por el azar. Los resultados del test de Shapiro-Wilk figuran debajo, donde se obtuvo como valor del estadístico 0’9938 y como nivel crítico 0’931. Se acepta por tanto la normalidad.

> shapiro.test(z) Shapiro-Wilk normality test data: z W = 0.9938, p-value = 0.931 En el gráfico situado en la esquina superior derecha se presenta el QQ-Plot asociado a la muestra de observaciones exponenciales. Los puntos deberían encontrarse sobre la línea, si procedieran de la distribución normal. Sin embargo, la falta de simetría hace que los cuantiles grandes estén muy por encima de la recta. Los cuantiles pequeños también se encuentran por encima de la media, aunque no tan distanciados. Los dos gráficos inferiores representan distribuciones simétricas. En ambos casos, los cuantiles de las colas no coinciden con los esperados bajo distribución normal. La diferencia es que la distribución uniforme tiene un comportamiento más regular, con un conjunto amplio de observaciones que discrepan de la normalidad. Por el contrario, en la distribución T de Student casi todas las observaciones se asemejan bastante bien a la distribución normal, mientras que de vez en cuando unas pocas observaciones pueden adoptar valores extremos incompatibles con la normalidad. En los tres casos los niveles críticos conducen al rechazo de la normalidad. No reproducimos todas las salidas, que puede obtener el lector ejecutando el código R, que se acompaña en el anexo a este tema.

Tema 3. Contraste de la normalidad multivariante

9

Datos exponenciales

4 3 2 0

1

Sample Quantiles

1 0 −1 −2

Sample Quantiles

2

5

Datos normales

−1

0

1

2

−2

−1

0

1

Theoretical Quantiles

Theoretical Quantiles

Datos uniformes

Datos T de Student

2

2 0 −4

−2

Sample Quantiles

0.8 0.6 0.4 0.2 0.0

Sample Quantiles

1.0

−2

−2

−1

0

1

Theoretical Quantiles

2

−2

−1

0

1

2

Theoretical Quantiles

Figura 4: QQ-Plots de datos simulados de distintas distribuciones.

10

Análisis Multivariante

3. Contraste de la normalidad multivariante En esta sección veremos cómo se pueden extender los conceptos de asimetría y kurtosis al caso multivariante. Asimismo, también se verá la manera de extender el test de Shapiro-Wilk con datos multivariantes. Empezamos estudiando ciertos conceptos previos en relación con este tipo de contrastes.

3.1. Conceptos previos: invariancia, estandarización, y distancias y ángulos de Mahalanobis Sabemos que la normalidad se conserva frente a traslaciones y transformaciones lineales, esto es: Para cualesquiera b vector de dimensión d y A matriz d × d no singular,

X tiene distribución normal ⇐⇒ AX + b tiene distribución normal Por tanto, si estamos contrastando que el vector X tenga distribución normal, tanta evidencia habrá a favor o en contra de esta hipótesis, si empleamos los datos de X como si empleamos los datos transformados AX + b. Definición . Diremos que un test de normalidad multivariante es invariante si proporciona los mismos resultados con los datos originales que con datos transformados por traslaciones o transformaciones lineales. Nos parece que la invariancia es una propiedad deseable para un test de normalidad. Ahora, si somos libres de transformar los datos, se nos ocurre aplicar la estandarización multivariante, pues es una transformación que los simplifica al suprimirles la media y la matriz de covarianzas. Pasamos entonces de los datos originales X1 , . . . , Xn a los datos estandarizados:

( ) ¯ Zi = S −1/2 Xi − X ¯ la media muestral y S la matriz de covarianzas muestral. siendo X La idea es que los datos Z1 , . . . , Zn deberían presentar una distribución normal estándar, esto es, deberían mostrar un aspecto esférico, sin especial predilección por unas direcciones frente a otras, y con mayor concentración en la zona central. Las distancias de Mahalanobis de las observaciones al vector de medias se pueden expresar así:

( ) ( ) ¯ S −1 Xi − X ¯ = Z ′ Zi = ∥Zi ∥2 rii = Xi − X i Ahora consideramos los valores

( ) ( ) ¯ S −1 Xj − X ¯ = Zi′ Zj = ∥Zi ∥ · ∥Zj ∥ · cos(θij ) rij = Xi − X

11

2

Tema 3. Contraste de la normalidad multivariante

1

(a)

0

(b)

−2

−1

−(b)

−3

−2

−1

0

1

2

Figura 5: Representación de los valores rij .

siendo θij el ángulo que forman las observaciones estandarizadas Zi y Zj en el espacio usual. En este momento definimos θij como el ángulo de Mahalanobis que forman las observaciones Xi y Xj en el espacio original. En la Figura 5 se representa una muestra de observaciones, que han sido estandarizados. Por tanto, las distancias y los ángulos de Mahalanobis coinciden con las distancias y ángulos ordinarias en esta representación. El punto central, representado en color azul, es la posición del vector de medias. En el gráfico se han destacado también dos puntos, denotados por (a) y (b), respectivamente. Las longitudes de las flechas, que son sus distancias al vector de medias, serían las raíces cuadradas de sus valores rii y rjj . El valor rij para este par de puntos sería el producto de las longitudes, multiplicado por el coseno del ángulo que forman las dos flechas. Así, si el ángulo es de 90 grados,

rij vale cero, si es menor de 90 grados rij es positivo, y si es mayor de 90 grados es negativo. De este modo, rij refleja si los dos puntos (a) y (b) están del mismo lado o si se encuentran en lados opuestos. Esto será utilizado en la sección siguiente para definir una medida de asimetría.

3.2. Tests basados en medidas de asimetría y kurtosis Existen diversas formas de definir medidas de asimetría y kurtosis con datos multivariantes. Consideraremos aquí las introducidas por Mardia.

12

Análisis Multivariante

Coeficiente de asimetría multivariante Así, calcularemos la asimetría de la siguiente manera:

Am =

n n 1 ∑∑ 3 rij n2 i=1 j=1

En la Figura 5 se ofrece una representación que resulta muy útil para interpretar los valores rij . Como vimos en la sección anterior, rij es positivo si los dos puntos en cuestión, (a) y (b) en el gráfico, están del mismo lado (ángulo inferior a 90 grados). Por el contrario, rij será negativo si se encuentran en lados opuestos respecto del vector de medias. Siempre se cumple que

n ∑

rij = 0 =

i=1

n ∑

rij

j=1

esto es, se cancelan los valores positivos con los negativos, lo cual es muy comprensible, pues el vector de media ocupa la posición central. Esto también se cumplía en el caso univariante. Sin embargo, al elevar al cubo los rij , se desequilibra este balance. Para mantener el equilibrio, por cada punto, como el (b), tendría que haber otro punto, que se representa como -(b) en el gráfico, situado en la posición contraria a (b) respecto del vector de medias. Los puntos (b) y -(b) tendrían los mismos valores rij con el resto de la muestra, pero con signos cambiados. Al efectuar cualquier potencia impar, como el cubo, se cancelarían. Así, si cada punto de la muestra tuviera otro punto simétrico en la propia muestra, el coeficiente de asimetría valdría cero. En otro caso, si la simetría se cumple aceptablemente, el coeficiente sería sólo algo mayor que cero, mientras que si hay un comportamiento claramente distinto en lados opuestos, el coeficiente de asimetría sería más grande. En la Figura 6 se presentan seis ejemplos de muestras bivariantes de tamaño 40. Mediante puntos huecos de color negro se representan los datos muestrales, con un punto relleno de color azul se representa el vector de medias y con una elipse de color azul se representa la matriz de covarianzas, siguiendo el procedimiento del documento de Preliminares. El gráfico (b) contiene 40 puntos situados en una circunferencia centrada en cero. El coeficiente de asimetría vale cero. El gráfico (a) representa datos simulados de la distribución normal. A nivel teórico, la distribución normal multivariante tiene un coeficiente de asimetría teórico igual a cero. Sin embargo, a nivel muestral dicho coeficiente puede ser algo mayor que cero. Los datos del gráfico (a) presenta un valor de 0’299 para el coeficiente de asimetría multivariante, Am . En el gráfico (c) se representan datos bivariantes con marginales exponenciales e independientes. La asimetría del modelo se refleja en un coeficiente de asimetría muestral de 15’04. En el gráfico (d) se han simulados tres grupos de datos normales multivariantes. Esto hace que la distribución de la muestra completa no sea normal. De hecho, la presencia de varios grupos es una

Tema 3. Contraste de la normalidad multivariante

13

causa muy común de falta de normalidad. También es común que produzca un efecto de asimetría, que permite detectar esa falta de normalidad. En este caso, el coeficiente Am ha tomado el valor 2’45. En el gráfico (e) se representó una muestra simulada cuyas marginales tienen distribución de Laplace y son independientes. Entonces la distribución teórica es simétrica. Sin embargo el coeficiente de asimetría muestral Am tomó el valor 6’14, lo cual es grande. La causa se encuentra en algunos valores extremos, que son frecuentes en un modelo Laplace, y que desequilibran la muestra. Recordemos que la distribución de Laplace se puede ver como diferencia de dos exponenciales independientes. En el gráfico (f) se han simulado observaciones uniformes en el cuadrado unidad. El modelo teórico es simétrico y de hecho el valor muestral de Am ha sido 0’484. Como test de normalidad, y dado que la distribución normal es simétrica, rechazaremos la normalidad, por falta de simetría, si el coeficiente Am es demasiado grande, comparado con ciertos valores tabulados correspondientes a muestras normales. Aproximaremos estos valores por simulación. Los detalles sobre la aplicación práctica del test se encuentran en la sección 3.4. Coeficiente de kurtosis multivariante El coeficiente de kurtosis propuesto por Mardia adopta la forma siguiente:

1∑ 2 = rii n n

Km

i=1

A la vista de esta definición, el coeficiente de kurtosis está basado únicamente en los valores rii . 2 que Como rii es el cuadrado de la distancia del dato i-ésimo al vector de medias, los valores rii

figuran en la definición de Km son las potencias cuartas de estas distancias. Por tanto, el coeficiente Km generaliza la medida de kurtosis al caso multivariante como estas potencias cuartas de la distancia a la media. De esta manera, para la kurtosis no importará la dirección o sentido de la desviación respecto de la media, sino únicamente la magnitud de dicha desviación. Se trata, pues, de detectar si los datos se agrupan en torno a la Media±Desviación típica, que en el caso multivariante sería la elipse en azul de la Figura 6, o si por el contrario se agrupan o muy cerca de la media o muy lejos de ésta. Si están cerca de la elipse la kurtosis será pequeña (distribución platicúrtica), y si están o muy cerca de la media o muy lejos de ella, la kurtosis será grande y diremos que la distribución es leptocúrtica. Igual que en el caso univariante, la distribución normal multivariante presenta un valor intermedio de kurtosis (diríamos que es mesocúrtica). Así, en el gráfico (a) de la Figura 6 el valor del coeficiente de kurtosis Km es 7’56. Los gráficos (b), (d) y (f) presentan valores menores, serían platicúrticos, mientras que los gráficos (c) y (e) son leptocúrticos. Respecto del contraste de normalidad, el test basado en la kurtosis rechazará la normalidad tanto

Análisis Multivariante

−1.0

−2

−1

0.0

0

0.5

1

1.0

14

−2

−1

0

1

2

−1.0 −0.5

0.5

1.0

(b) Circunferencia

8

0.0

10

1.0

12

14

2.0

16

3.0

(a) Normales

0.0

0

2

4

6

8

14

18

(d) Tres grupos

0.0

−4 −2

0

0.4

2

4

0.8

6

(c) Exponenciales

10

−2

0

2

4

6

(e) Distribuciones de Laplace

0.0

0.2

0.4

0.6

0.8

1.0

(f) Uniforme en el cuadrado

Figura 6: Ejemplos de muestras bivariantes para su contraste de normalidad.

Tema 3. Contraste de la normalidad multivariante

15

para valores demasiado grandes como para valores demasiado pequeños de Km . De nuevo, los valores críticos están tabulados y los aproximaremos por simulación siguiendo las pautas de la sección 3.4.

3.3. Extensiones multivariantes del test de Shapiro-Wilk Existen diversas maneras de extender el test de Shapiro-Wilk. Nos centraremos en una idea muy sencilla, que surge del concepto de invariancia. Así, el vector aleatorio es normal multivariante si y sólo si su estandarización lo es, esto es:

X ∈ Nd (µ, Σ) ⇐⇒ Z = Σ−1/2 (X − µ) ∈ Nd (0, Id ) Lo interesante es que Z tiene sus d componentes (las d variables) incorrelacionadas, y entonces Z es normal si y sólo si sus componentes son normales. Ante tal resultado, basta con comprobar si las d variables resultantes de la estandarización multivariante son normales. Volvemos entonces a los datos estandarizados:

( ) ¯ Zi = S −1/2 Xi − X ¯ la media muestral y S la matriz de covarianzas muestral. siendo X Si W1 , . . . , Wd son los estadísticos de Shapiro-Wilk de cada componente de Z1 , . . . , Zn , entonces podemos considerar el estadístico de Shapiro-Wilk multivariante:

1∑ = Wj d d

WV G

j=1

Este estadístico ha sido propuesto por Villasenor Alva y González Estrada (2009). Se rechazará la normalidad cuando WV G sea pequeño, en comparación con ciertos valores tabulados, que de nuevo podemos aproximar por simulación, según se indica en la sección 3.4, que viene a continuación.

3.4. Aplicación práctica de los test de normalidad multivariante Para la aplicación práctica de los tests de normalidad multivariante hemos programado tres funciones en código R, que figuran en los anexos al final del tema: representa (x,texto) A partir de una matriz x de observaciones bivariantes (dos columnas), obtiene un gráfico de dispersión, acompañado del vector de media y la matriz de covarianzas.

16

Análisis Multivariante

estadisticos (x) A partir de una matriz x de observaciones multivariantes (con cualquier número de columnas), calcula los estadísticos Am , Km y WV G . pvalores (d,n,ns,am0,km0,w0) A partir de los valores d=número de variables, n=tamaño muestral, ns=número de muestras simuladas, am0=valor de Am en la muestra, km0=valor de Km en la muestra y w0=valor de WV G en la muestra, proporciona una tabla aproximada por simulación con los cuantiles de estos estadísticos bajo la hipótesis de normalidad, y aproxima también por simulación los niveles críticos para el contraste de normalidad. Para poder emplear estas funciones, habría que ejecutar el código correspondiente a cada una de ellas. La ejecución de ese código no genera salidas ni cálcuo alguno, sino que “compila" el código, como paso previo para poder llamar a las funciones. Para la ejecución de los tests sólo son necesarias las funciones estadisticos y pvalores, mientras que la función representa sirve para representar los datos en el caso bivariante. En el último anexo se encuentra el código para la obtención de los seis ejemplos representados en la Figura 6, así como para el contraste de normalidad sobre ellos. Este código contiene una parte dedicada a la simulación de las muestras y definición de variables. Para realizar los tests en otra muestra, no necesariamente obtenida por simulación, bastaría con ejecutar las funciones:

> estadisticos(x) > est > pvalores(d,n,ns,est[1],est[2],est[3]) donde las variables x, d, n y ns deben contener los valores correspondientes. Sobre los seis ejemplos simulados, comentaremos la salida de resultados para el último ejemplo. Los demás siguen la misma estructura. El vector est contiene los valores del coeficiente de asimetría Am , el coeficiente de kurtosis Km y el estadístico de Shapiro-Wilk multivariante WV G . Los vemos a continuación.

> est # Estadísticos Asimetría 0.4844482

Kurtosis 6.2222052

Shapiro-Wilk 0.9538382

Para saber si estos valores son grandes o pequeños, de cara al contraste de normalidad, se comparan con los cuantiles en 10.000 muestras simuladas de una distribución normal. Los resultados figuran debajo. En primer lugar se presentan los cuantiles. Los órdenes de los cuantiles son 0.01,

Tema 3. Contraste de la normalidad multivariante

17

0.025, 0.05, 0.1, 0.5, 0.9, 0.95, 0.975, 0.99. La mediana es el cuantil 0.5, que sería un valor de referencia para saber si el estadístico en cuestion (Asimetría, Kurtosis o Shapiro-Wilk) es grande o pequeño. Así, la asimetría de esta muestra es 0’484 lo cual es superior a 0’427, que es la mediana bajo normalidad. Por tanto, esta muestra es algo más asimétrica que lo que corresponde bajo normalidad. Para saber si este resultado es significativo como para rechazar la normalidad, recordamos que se rechaza la normalidad si la asimetría es grande. Por tanto, el nivel crítico sería la probabilidad que deja 0’484 a su derecha. Esa probabilidad, nivel crítico, vale 0’4386. Se aproximó por simulación considerando que 4.386 de las 10.000 muestras simuladas tenían una asimetría superior al valor 0’484 obtenido en la muestra original. Para la kurtosis, su valor en la muestra original ha sido de 6’22, lo cual es inferior a 7’49, que es la mediana bajo normalidad. Por tanto, esta distribución es algo más platicúrtica que lo que correspondería bajo normalidad. Para efectuar el contraste recordamos que se rechaza la normalidad tanto por exceso como por defecto de kurtosis. El valor muestral 6’22 se encuentra entre el cuantil 0’025 y el cuantil 0’05. El nivel crítico unilateral se hallará entre esos dos órdenes, mientras que el bilateral estará entre 0’05 y 0’10. El cálculo más preciso haciendo el recuento exacto sobre las muestras simuladas arroja un nivel crítico bilateral aproximado de 0’0948. Para el estadístico de Shapiro-Wilk multivariante, el valor en la muestra original es de 0’954. Cuanto menor sea este estadístico, peor es el ajuste de normalidad, y se rechaza la normalidad cuando el estadístico es muy pequeño. El valor 0’954 es inferior a la mediana bajo normalidad, que vale 0’973. Además, se encuentra entre el cuantil 0’05 y el cuantil 0’1. Por tanto el nivel crítico se hallará entre estos dos órdenes. Su valor aproximado de manera más precisa con las 10.000 muestras simuladas es de 0’0566, pues sólo 566 de las 10.000 muestras simuladas presentaban un valor del estadístico inferior al de la muestra original.

> pvalores(d,n,ns,est[1],est[2],est[3]) $tabla 0.01 0.025 0.05 0.1 0.5 0.9 Asimetría 0.03549275 0.05503973 0.0847581 0.1253795 0.4266827 1.0577071 Kurtosis 5.87278811 6.04901296 6.2419949 6.4783751 7.4851531 8.8914292 Shapiro-Wilk 0.94019144 0.94735872 0.9529298 0.9582334 0.9733307 0.9829626 0.95 0.975 0.99 Asimetría 1.327830 1.607430 1.9728774 Kurtosis 9.429676 9.986294 10.6976297 Shapiro-Wilk 0.985076 0.986636 0.9882852 $Niveles_críticos Asimetría Kurtosis Shapiro-Wilk 0.4386 0.0948 0.0566

18

Análisis Multivariante

Los resultados de los estadísticos y los niveles críticos para los seis ejemplos se obtienen ejecutando el código R del último anexo. Sólo la muestra (a) procede de una distribución normal. Los tests rechazan la normalidad de las demás muestras, excepto en el caso (f), cuya desviación de la normalidad es menor. Al igual que en el caso univariante, el test de Shapiro-Wilk es el mejor como método general de contraste de normalidad, mientras que los tests basados en la asimetría o la kurtosis son eficaces para detectar algunas alternativas más específicas.

Bibliografía. Henze, N. (2002). Invariant tests for multivariate normality: a critical review. Statistical Papers, 43, 467–506. Mardia, K.V. (1975). Assessment of multinormality and the robustness of Hotelling’s T 2 test. Applied Statistics, 24, 163–171. Villasenor Alva, J.A. y González Estrada, E. (2009). A generalization of Shapiro-Wilk test for multivariate normality. Communications in Statistics – Theory and Methods, 38, 1870–1883.

Tema 3. Contraste de la normalidad multivariante

Código R para los ejemplos Ejemplo 1. Test de asimetría en datos univariantes set.seed(123456) #--- Generación de observaciones normales n=30 z=rnorm(n) hist(z) points(z,rep(0,n)) #--- Cálculo de la media, la varianza y el coeficiente de asimetría med=mean(z) s2=var(z)*(n-1)/n s=sqrt(s2) a_normal=mean(( (z-med)/s )^3) points(med,0,col="blue",pch=19) #--- Test basado en el coeficiente de asimetría estadistico_normal=sqrt(n/6)*a_normal pvalor_normal = 2 *( 1 - pnorm( abs(estadistico_normal) ) ) windows() par(mfrow=c(1,2)) #--- Generación de observaciones exponenciales u=runif(n) x=-log(u) hist(x) points(x,rep(0,n)) #--- Cálculo de la media, la varianza y el coeficiente de asimetría med=mean(x) s2=var(x)*(n-1)/n s=sqrt(s2) a_exponencial=mean((x-med)^3)/s^3 points(med,0,col="blue",pch=19) #--- Test basado en el coeficiente de asimetría estadistico_exponencial=sqrt(n/6)*a_exponencial

19

20

Análisis Multivariante

pvalor_exponencial = 2 *( 1 - pnorm( abs(estadistico_exponencial) ) ) #--- Se cambia el signo a las observaciones x=-x hist(x) points(x,rep(0,n)) #--- Cálculo de la media, la varianza y el coeficiente de asimetría med=mean(x) s2=var(x)*(n-1)/n s=sqrt(s2) a=mean(( (x-med)/s )^3) points(med,0,col="blue",pch=19) par(mfrow=c(1,1)) #--- Presentación de resultados c(a_normal,a_exponencial,-a_exponencial) c(estadistico_normal,estadistico_exponencial,-estadistico_exponencial) c(pvalor_normal,pvalor_exponencial,pvalor_exponencial)

Ejemplo 2. Test de kurtosis en datos univariantes set.seed(123456) n=100 #--- Generación de observaciones normales z=rnorm(n) #--- Test basado en el coeficiente de kurtosis med=mean(z) s2=var(z)*(n-1)/n s=sqrt(s2) k_normal=mean( ( (z-med)/s )^4) - 3 estadistico_normal=sqrt(n/24)*k_normal pvalor_normal = 2 *( 1 - pnorm( abs(estadistico_normal) ) ) #--- Generación de observaciones uniformes z=runif(n)

Tema 3. Contraste de la normalidad multivariante

#--- Test basado en el coeficiente de kurtosis med=mean(z) s2=var(z)*(n-1)/n s=sqrt(s2) k_uniforme=mean( ( (z-med)/s )^4) - 3 estadistico_uniforme=sqrt(n/24)*k_uniforme pvalor_uniforme = 2 *( 1 - pnorm( abs(estadistico_uniforme) ) ) #--- Generación de observaciones T de Student z=rt(n,5) #--- Test basado en el coeficiente de kurtosis med=mean(z) s2=var(z)*(n-1)/n s=sqrt(s2) k_tstudent=mean( ( (z-med)/s )^4) - 3 estadistico_tstudent=sqrt(n/24)*k_tstudent pvalor_tstudent = 2 *( 1 - pnorm( abs(estadistico_tstudent) ) ) #--- Presentación de resultados c(k_normal,k_uniforme,k_tstudent) c(estadistico_normal,estadistico_uniforme,estadistico_tstudent) c(pvalor_normal,pvalor_uniforme,pvalor_tstudent)

Ejemplo 3. Test de Shapiro-Wilk con datos univariantes set.seed(123456) n=100 windows() par(mfrow=c(2,2)) #--- Generación de observaciones normales z=rnorm(n) qqnorm(z,main="Datos normales") qqline(z) shapiro.test(z) #--- Generación de observaciones exponenciales

21

22

Análisis Multivariante

z=rexp(n) qqnorm(z,main="Datos exponenciales") qqline(z) shapiro.test(z) #--- Generación de observaciones uniformes z=runif(n) qqnorm(z,main="Datos uniformes") qqline(z) shapiro.test(z) #--- Generación de observaciones T de Student z=rt(n,5) qqnorm(z,main="Datos T de Student") qqline(z) shapiro.test(z) par(mfrow=c(1,1))

representa: Función R que representa datos bivariantes, junto con el vector de medias y la matriz de covarianazas representa=function(x,texto){ # Los datos se encuentran en la matriz x #d=ncol(x) # Número de variables=2 n=nrow(x) # Número de individuos #--- Cálculo del vector de medias y matriz de covarianzas med=apply(x,2,mean) sc=cov(x) s=sc*(n-1)/n #--- Diagonalización de s auto=eigen(s) v=auto$vectors lambda=auto$values #--- Representación

Tema 3. Contraste de la normalidad multivariante

23

plot(x,xlab=texto,ylab="") points(med[1],med[2],pch=19,col="blue") tita=seq(0,2*pi,length=101) # Vector con los ángulos medr=matrix(rep(med,101),byrow=TRUE,nrow=101) # Truco para repetir el vector de medias diez veces, en 101 filas elipse0=medr+t(sqrt(lambda[1])*v[,1]%*%t(cos(tita)) +sqrt(lambda[2])*v[,2]%*%t(sin(tita))) lines(elipse0,col="blue") }

estadisticos: Función R para el cálculo de la asimetría, la kurtosis y el estadístico de Shapiro-Wilk en sus versiones multivariantes estadisticos=function(x){ # Los datos se encuentran en la matriz x d=ncol(x) # Número de variables n=nrow(x) # Número de individuos #--- Cálculo del vector de medias y matriz de covarianzas med=apply(x,2,mean) sc=cov(x) s=sc*(n-1)/n #--- Diagonalización de s auto=eigen(s) v=auto$vectors lambda=auto$values si12=v%*%diag(1/sqrt(lambda))%*%t(v) #--- Estandarización multivariante xs=(x-rep(1,n)%*%t(med))%*%si12 #--- Similaridades r=xs%*%t(xs) #--- Cálculo de los estadísticos am=sum(r^3)/(n*n) km=sum(diag(r)^2)/n

24

Análisis Multivariante

w=0 for (j in 1:d){s=shapiro.test(xs[,j]) w=w+s$statistic} w=w/d salida=c(am,km,w) names(salida)=c("Asimetría","Kurtosis"," Shapiro-Wilk") salida}

pvalores: Función R que aproxima por simulación los cuantiles de los estadísticos multivariantes, y los niveles críticos pvalores=function(d,n,ns,am0,km0,w0){ set.seed(123456) am=c() # Inicialización de vectores km=c() w=c() for(is in 1:ns){ #--- Generación de un vector normal multivariante x=matrix(rnorm(d*n),ncol=d) #--- Cálculo del vector de medias y matriz de covarianzas med=apply(x,2,mean) sc=cov(x) s=sc*(n-1)/n #--- Diagonalización de s auto=eigen(s) v=auto$vectors lambda=auto$values si12=v%*%diag(1/sqrt(lambda))%*%t(v) #--- Estandarización multivariante xs=(x-rep(1,n)%*%t(med))%*%si12 #--- Similaridades

Tema 3. Contraste de la normalidad multivariante

r=xs%*%t(xs) #--- Cálculo de los estadísticos am[is]=sum(r^3)/(n*n) km[is]=sum(diag(r)^2)/n w[is]=0 for (j in 1:d){s=shapiro.test(xs[,j]) w[is]=w[is]+s$statistic} w[is]=w[is]/d } p=c(0.01,0.025,0.05,0.1,0.5,0.9,0.95,0.975,0.99) # Órdenes de los cuantiles am=sort(am) qam=am[p*ns] km=sort(km) qkm=km[p*ns] w=sort(w) qw=w[p*ns] tabla=rbind(qam,qkm,qw) colnames(tabla)=p rownames(tabla)=c("Asimetría","Kurtosis","Shapiro-Wilk") pvalor_am0=sum(am>am0)/ns pvalor_km0=sum(km>km0)/ns pvalor_km0=2*min(pvalor_km0,1-pvalor_km0) pvalor_w0=sum(w