Diagnóstico de normalidad con R1

2
Mario Alfonso Morales Rivera

October 22, 2019

1R Development Core Team (2006). R: A language and environment for statistical computing. R
Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-
project.org.
2 Profesor titular, Departamento de Matemáticas y Estadı́stica Universidad de Córdoba
Introducción

La mayorı́a de procedimientos estadı́sticos que se usan habitualmente suponen que los datos
observados proceden de una población con distribución normal. Una razón de ello es que
muchas variables asociadas a fenómenos naturales y sociales siguen aproximadamente esta
distribución, otra razón (quizás la más importante) del uso extendido del supuesto de nor-
malidad es la facilidad y elegancia con que se obtienen los estimadores y los procedimientos
para la inferencia y prueba de hipótesis. Aunque muchas de las técnicas estadı́sticas son poco
sensibles a la violación del supuesto de normalidad, (en general este supuesto puede obviarse
cuando se cuenta con un tamaño de muestra grande – resultados asintóticos–) es recomend-
able contrastar siempre si se puede asumir o no una distribución normal. El diagnóstico del
supuesto de normalidad incluye desde la simple exploración visual de los datos hasta técni-
cas estadı́sticas sofisticadas que ayudan a decidir si es razonable suponer la normalidad del
conjunto de datos en cuestión.
En este documento se hace una revisión de la literatura y se resumen las técnicas más ex-
tendidas para diagnosticar normalidad. El documento se divide en dos grandes secciones:
Caso univariado y caso multivariado. En cada caso se estudian procedimientos gráficos y las
estadı́sticas para realizar la prueba formal del supuesto de normalidad.
El único aporte significativo de mi parte es la revisión y recopilación de las diferentes opciones
(funciones y librerı́as) que tiene R para la realización de los diferentes gráficos y pruebas, todos
los comandos se transcriben en el documento y pueden ser ejecutados por el lector.

1
Caso univariado

Sabemos que si un vector aleatorio X es normal p variado, entonces cada una de las vari-
ables aleatorias componentes es normal univariada. Por tanto si por lo menos una de las
componentes de un vector aleatorio no es normal, podemos asegurar que éste no es nor-
mal multivariado, de ahı́ la importancia de conocer las técnicas usadas para diagnosticar
normalidad en el caso univariado.

Métodos gráficos
Para el diagnóstico de normalidad en el caso univariado se han desarrollado varias estrategias
gráficas que de manera visual alertan sobre la normalidad o no de los datos.

Histograma
Es una herramienta sencilla de implementar, porque todos los paquetes estadı́sticos tienen
programas para elaborarlos. Si los datos son normales, el histograma deberı́a mostrar la bien
conocida forma acampanada. En la figura 1 se muestran 4 histogramas elaborados a partir de
datos simulados, ¿ podrı́as señalar los que provienen de datos normales? El siguiente código
genera el gráfico que se muestra en la figura 2 a partir de datos simulados. Para unos datos
en particular, solo tiene que introducirlos en el vector x de la forma x<-c(4.5,3.9,...)
x<-rnorm(100,mean=10,sd=1) # introduzca sus datos en x
# el histograma
histo<-hist(x,prob=T,main="",ylab="",ylim=range(hist(x)$density))
# la densidad
z<-pretty(histo$breaks,n=50)
y<-dnorm(z,mean=mean(x),sd=sd(x))
lines(z,y,lty=3,lwd=2,col="blue")

# para más detalles de la función hist

# digite en la consola ?hist

Gráficos Q − Q
Otro procedimiento gráfico para verificar normalidad univariada es el gráfico de proba-
bilidad normal. Este es un gráfico de los cuantiles empı́ricos contra los cuantiles teóricos

2
 Figure 1: Histogramas a partir de datos generados de varias distribuciones

(de ahı́ el nombre Q–Q plot) de la distribución normal estándar. Cuando los puntos en el
gráfico de probabilidad normal quedan cerca de una linea recta, el supuesto de normalidad es
razonable. El patrón de la desviación de los puntos de una linea recta indica la naturaleza de
la separación de la normalidad tal como asimetrı́a, apuntamiento, datos extremos o multiples
modas. Para varios gráficos tı́picos que muestran separación de la normalidad véase Rencher
(1995) página 106.
La mayorı́a de los paquetes estadı́sticos facilitan la elaboración de este gráfico, sin embargo
aquı́ se ilustrará paso a paso su construcción y luego se darán los comandos de R para
generarlo.
1. Ordene las observaciones y denote los valores ordenados por y(1) , y(2) , · · · , y(n) ; de esa
forma y(1) ≤ y(2) ≤ · · · ≤ y(n) entonces el punto y(i) es cuantil muestral i/n. A menudo
la fracción i/n se cambia por (i − 21 )/n como una corrección por continuidad, de esta
forma y(i) se designa como el (i − 12 )/n cuantil muestral.
2. Calcule los cuantiles poblacionales q1 , q2 , · · · , qn , donde qi es el valor para el cual la
probabilidad de obtener una observación menor o igual que él es igual a (i − 21 )/n, es

3
 Figure 2: Histograma con la densidad normal superpuesta

decir, qi es tal que

1
i− 2
P (Z < qi ) =
n
con Z normal estándar.

3. Grafique los pares (qi , y(i) ) y examine la linealidad de los puntos.

Para ilustrar este procedimiento, y todos los de esta sección, usaremos los datos 1.38 1.40
1.42 1.54 1.30 1.55 1.50 1.60 1.41 1.34 . En la tabla 1 se muestran los resultados de los
cálculos. Con el código de R que está a continuación se obtiene el gráfico cuantil–cuantil que
se muestra en la figura 3.

y<-c(1.38,1.40,1.42,1.54,1.30,1.55,1.50,1.60,1.41,1.34)
qqnorm(y)
qqline(y) # pasa la linea

4
 Figure 3: Gráfico de cuantil–cuantil para verificar normalidad.
1
y(i) (i − 2
)/10 qi
1.30 0.05 −1.645
1.34 0.15 −1.036
1.38 0.25 −0.674
1.40 0.35 −0.385
1.41 0.45 −0.126
1.42 0.55 0.126
1.50 0.65 0.385
1.54 0.75 0.674
1.55 0.85 1.036
1.60 0.95 1.645

Tabla 1: Datos ordenados, cuantiles muestrales y cuantiles poblacionales

Contrastes de normalidad
En la literatura estadı́stica se han propuesto varios procedimientos analı́ticos para probar
la normalidad de datos univariados, aquı́ revisaremos la prueba Ji–cuadrado de bondad de
5
ajuste, la prueba de Kolmogorov–Smirnov, la prueba de Shapiro–Wilk y las pruebas basadas
en asimetrı́a y curtosis.

Prueba Ji–cuadrado
Esta prueba es útil para probar el ajuste de un conjunto de datos a cualquier distribución.
Se basa en el estadı́stico
k
X (Oi − Ei )2
χ20 =
i=1
Ei
donde Oi son las frecuencias observadas en las k clases [x0 , x1 ), · · · [xk−1 , xk ] y Ei son las
frecuencias esperadas según el modelo probabilı́stico propuesto, para el caso normal se tiene
Ei = npi con pi = P (xi−1 ≤ X ≤ xi ).
La estadı́stica χ20 se distribuye aproximadamente como una Ji–cuadrado con k − r − 1 grados
de libertad, donde r es el numero de parámetros que se estiman, en el caso de la normal r = 2
porque se estima la media y la varianza.
Con los siguientes comandos de R se realiza esta prueba, tenga en cuenta que hay que instalar
la librerı́a nortest1 .

y<-c(1.38,1.40,1.42,1.54,1.30,1.55,1.50,1.60,1.41,1.34)
library(nortest)
pearson.test(y)

La salida se muestra a continuación, el p−valor indica que no hay suficiente evidencia para
rechazar la hipótesis de normalidad. Tenga en cuenta, sin embargo, que esta prueba es poco
potente con tamaños de muestra pequeños.

Pearson chi-square normality test

data: y
P = 3.2, p-value = 0.3618

Kolmogorov–Smirnov
Se asume que tenemos una muestra aleatoria X1 , X2 , · · · , Xn de alguna distribución continua
con función de distribución acumulada F (·). Denotamos la función de distribución acumulada
empı́rica por
1
FN (x) = (numero de obs ≤ x)
N
1
Si tu computador está conectado a internet, es muy fácil instalar una librerı́a: en el menú Packages click
en install package(s)..., primero aparece una ventana para que selecciones el servidor (se recomienda
que escojas uno que esté cercano a tu pais), click en OK y luego aparece una ventana para que selecciones
la librerı́a (package) que deseas instalar, click on OK. Después de unos segundos, si no hay problemas con la
conexión, el paquete estará disponible y puedes cargarlo con el comando library()

6
la prueba de Kolmogorov–Smirnov se utiliza para probar H0 : F (x) = F0 (x) para todo x
contra H1 : F (x) 6= F0 (x) para algún x, donde F0 a una distribución N (µ, σ 2 ). El estadı́stico
de Kolmogorov–Smirnov es
DN = sup |FN (x) − F0 (x)|
x

y es grande si los datos no son consistentes con H0 . La distribución asintótica de DN , bajo

H0 cierta es √
lim P { N DN ≤ x} = Q(z)
x→∞
con ∞
X
Q(z) = 1 − 2 (−1)k−1 exp{−2k 2 z 2 }
k=1

para cada z > 0. Q(z) es la función de distribución acumulada de una distribución continua
conocida como la distribución de Kolmogorov. En general, los parámetros µ y σ 2 son de-
sconocidos y se pueden reemplazar por su contraparte muestral. Con el siguiente código de
R se prueba, usando Kolmogorov–Smirnov , si los datos del ejemplo son normales con media
µ = 1.4 y σ = 0.1

y<-c(1.38,1.40,1.42,1.54,1.30,1.55,1.50,1.60,1.41,1.34)
ks.test(y,"pnorm", 1.4 , 0.1 )

La salida de la función ks.test() se muestra a continuación , el p–valor mayor que 0.05 indica
que no hay evidencia para rechazar la hipótesis que los datos pertenecen a la distribución
indicada.
data: y 1

D = 0.2413, p-value = 0.5285 2

alternative hypothesis: two.sided 3

Shapiro–Wilks
Esta prueba se basa en la comparación de los valores muestrales ordenados con su localización
esperada bajo la hipótesis nula de normalidad. Sea Z(1) , Z(2) , · · · , Z(n) una muestra ordenada
de una distribución normal estándar, y sea mi = E(Z(i) ), i = 1, · · · , N . Bajo la hipótesis de
normalidad,
E(X(i) ) = µ + σmi
es decir, se espera que las observaciones ordenadas X(i) ’s estén linealmente relacionadas a los
mi ’s.
El estadı́stico de Shapiro–Wilks es
" n #2
1 X
W = 2 aj,n (x(n−j+1) − x(j) )
ns j=1

7
donde s2 es la varianza muestral y
n

2
, si n es par
h= n−1
2
si n es impar

los coeficientes aj,n se consiguen en tablas y x(j) es el j−ésimo valor ordenado de la muestra.
Los comandos R para realizar esta prueba son los siguientes:

y<-c(1.38,1.40,1.42,1.54,1.30,1.55,1.50,1.60,1.41,1.34)
shapiro.test(y)

la salida de la función es la siguiente, como p−valor es mayor que 0.05 no rechazamos la

hipótesis nula H0 : los datos son normales.

Shapiro-Wilk normality test

data: y
W = 0.9519, p-value = 0.6911

Prueba basada en asimetrı́a y curtosis

Esta es una prueba clásica basada en las siguientes medidas de asimetrı́a y curtosis
√ P
n

p n (yi − y)3
i=1
b1 =  3/2
n
P
(yi − y)2
i=1
n
n (yi − y)4
P

b2 =  i=1 2
Pn
(yi − y)2
i=1
√
los cuales son estimadores de los coeficientes de asimetria y curtosis
√ poblacionales β1 y β2
respectivamente. Cuando la población es normal se tiene que β1 = 0 y β2 = √3
La prueba de normalidad basada en la asimetrı́a se lleva a cabo comparando b1 con valores
tabulados o alternativamente, cuando n ≥ 8, la función g definida por
√
p −1 b1
g( b1 ) = δ sinh
λ
√
tiene aproximadamente una distribución normal estándar, donde sinh−1 (x) = ln(x+ x2 + 1)
y los valores de λ y δ se obtienen de tablas.
Si los datos son normales, b2 tiene, de manera asintótica, distribución N (3, 24/n) y por
tanto tenemos una prueba de normalidad basada en la kurtosis: se rechaza la hipótesis de
normalidad si
|b2 − 3|
p > zα/2
24/n

8
Una prueba que usa simultáneamente la asimetrı́a y la kurtosis se basa en la estadı́stica

nb1 n(b2 − 3)2

X2 = +
6 24
la cual, bajo normalidad y asintóticamente se distribuye como una Ji–cuadrado con dos
grados de libertad. Rechazamos la hipótesis de normalidad si X 2 > χ22 (α). Con el siguiente
código de R se realizan las dos primeras pruebas, tenga en cuenta que la librerı́a moments no
pertenece al paquete básico de R ası́ que hay que instalarla manualmente.

y<-c(1.38,1.40,1.42,1.54,1.30,1.55,1.50,1.60,1.41,1.34)
library(moments)
agostino.test(y)
anscombe.test(y)

La salida de los comandos se muestra a continuación (respectivamente)

D'Agostino skewness test

data: y
skew = 0.1886, z = 0.2221, p-value = 0.8242
alternative hypothesis: data have a skewness

Anscombe-Glynn kurtosis test

data: y
kurt = 1.8164, z = -0.8954, p-value = 0.3706
alternative hypothesis: kurtosis is not equal to 3

En ambos casos no se rechaza la hipótesis de normalidad. Con siguiente código se implementa

la prueba omnibus.

library(moments)
b1<-skewness(y)^2
b2<-kurtosis(y)
n<-length(y)
X2<-(n*b1)/6+n*(b2-3)^2/24
pvalor<-pchisq(X2,2,lower.tail=FALSE)
cat("X2=",X2,"p_valor=",pvalor, "\n")

Para los datos de nuestro ejemplo la prueba arroja

X2= 0.6429844 p_valor= 0.7250663

con lo que se concluye que no hay evidencia para rechazar la hipótesis que los datos vienen
de una distribución normal.

9
Caso multivariado

La normalidad multivariante implica la normalidad de las distribuciones marginales, pero la

normalidad de las marginales no garantiza que la distribución conjunta sea normal, a menos
que las variables sean no correlacionadas, situación que es poco común. Ası́ que las pruebas
univariadas individuales sirven para demostrar que no hay normalidad, cuando al menos una
de ellas reporta no normalidad. Como en el caso univariado, estudiaremos procedimientos
gráficos y analı́ticos para el diagnostico de la multinormalidad.

Procedimiento gráfico
Como en el caso univariado, para el caso multivariado existen alternativas gráficas para
decidir si los datos provienen de una distribución normal, aquı́ estudiaremos dos de ellos.

Gráfico del tipo Q × Q

Suponga que x ∼ Np (µ; Σ), entonces (x − µ)0 Σ−1 (x − µ) tiene distribución Ji–cuadrado con
p grados de libertad.
Si se tiene una muestra aleatoria x1 , x2 , · · · , xn de una distribución Np (µ; Σ) entonces di =
(xi − µ)0 Σ−1 (xi − µ) se puede considerar una muestra aleatoria de una distribución Ji–
cuadrado con p grados de libertad, como µ y Σ no se conocen, se estiman. Se puede demostrar
que
nDi2
ui =
(n − 1)2
con Di2 = (xi − x)0 S −1 (xi − x) tiene distribución beta. Los valores ui se ordenan y situan
en una gráfica contra sus contrapartes teóricas dadas por
i−α p−2 n−p−2
vi = , con α = yβ=
n−α−β+1 2p 2(n − p − 1)

es decir, se grafican los pares (vi , u(i) ), si los puntos tienden a caer a lo largo de una linea
recta se concluye que los datos son normales multivariados, si los puntos no presentan esa
tendencia, entonces se concluye que los datos no son normales.
A continuación se transcribe el código para realizar este gráfico usando los datos de la tabla
3.4 de la página 102 de Dias Monroy (2002), los cuales se usarán para todos los ejemplos de
esta sección.

10
 Lectura de datos
Y1<-scan()
72 66 76 77 60 53 66 63 56 57 64 58 41 29 36 38 32 32 35 36 30 35 34
26 39 39 31 27 42 43 31 25 37 40 31 25 33 29 27 36 32 30 34 28 63 45
74 63 54 46 60 52 47 51 52 43 91 79 100 75 56 68 47 50 79 65 70 61
81 80 68 58 78 55 67 60 46 38 37 38 39 35 34 37 32 30 30 32 60 50 67
54 35 37 48 39 39 36 39 31 50 34 37 40 43 37 39 50 48 54 57 43

Y<-matrix(Y1,ncol=4,byrow=TRUE)
Y<-as.data.frame(Y)
Calculos
p<-ncol(Y)
n=nrow(Y)
S<-cov(Y)
Xbar<-apply(Y,2,mean)

Figure 4: Gráfico de tipo Q × Q para verificar multinormalidad.

11
 Di<-mahalanobis(Y,Xbar,S) # los D_i al cuadrado
Ui<-(n*Di)/(n-1)^2 # los U_i
Uio<-sort(Ui) # los U_i ordenados
alpha<-(p-2)/(2*p)
beta<-(n-p-2)/(2*(n-p-1))
vi<-(1:n-alpha)/(n-alpha-beta+1)
# genera el gráfico
plot(vi,Uio,xlab="Teóricos",ylab="Muestrales")

En la figura 4 se muestra el gráfico generado por el código anterior, como los puntos ajustan
pobremente a una linea recta, el gráfico está mostrando un posible desvı́o de la normalidad.

Gráficos por pares

Sabemos que si      
X1 µ1 Σ11 Σ12
X= ∼ Np ;
X2 µ2 Σ21 Σ22
entonces X 1 |X 2 = x2 se distribuye normal p1 variante con media µX 1 |X 2 = µ1 +Σ12 Σ−1 22 (x2 −
µ2 ) lo que significa que la media de X 1 , dado X 2 es una función lineal de X 2 por lo tanto
se garantiza que si la distribución conjunta es normal,

cada par de variable se ajusta a una
linea recta, si la nube de puntos para alguno de los p2 gráficos no muestra ajuste a una linea
recta, es una señal de no normalidad.
En la figura 5 se muestra, en el panel inferior, la nube de puntos para cada par de variables
junto con la linea de regresión ajustada y en el panel superior, el valor de R2 . Según el
gráfico, parece que no hay evidencia fuerte en contra de la normalidad de los datos, sin
embargo reiteramos que la normalidad de las marginales, en este caso bivariadas no implica
la normalidad conjunta.
A continuación se transcribe el código de programación para realizar el gráfico de la figura 5.
panel inferior
panel.lajus <- function(x,y, ...)
{
points(x,y)
abline(lm(y~x))
}
panel superior
panel.R2 <- function(x, y, digits=2, prefix="", cex.cor)
{
usr <- par("usr");
on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r<-summary(lm(y~x))$r.squared
txt <- format(c(r, 0.123456789),digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) cex <- 0.8/strwidth(txt)

12
 Figure 5: Gráfico de dispersión por pares, con linea de regresión ajustada y R2 .

text(0.5, 0.5, txt, cex = cex)

}
el gráfico
pairs(Y, upper.panel=panel.R2, lower.panel=panel.lajus,
cex.labels = 1, font.labels=1)

Contrastes de multinormalidad
En esta sección se discutirán algunas pruebas formales para probar si un conjunto de datos
multivariados viene de una distribución normal.

13
Prueba basada en la distancia de Mahalanobis
Esta prueba se basa en el máximo de la distancia de Mahalanobis de cada observación al
centroide de los datos. La estadı́stica de prueba es
2
D(n) = max{Di2 }
i

con Di2 = (xi − x)0 S −1 (xi − x) para tomar la decision, se compara D(n)
2
con valores que
se encuentran tabulados (tabla C.3 de Dias Monroy (2002) ). El código para calcular esta
estadı́stica es el siguiente:

S<-cov(Y)
Xbar<-apply(Y,2,mean)
D2_n<-max(mahalanobis(Y,Xbar,S))

Prueba de Mardia
Sea X 1 y X 2 independientes e identicamente distribuidos Np (µ; Σ), Mardia define los coefi-
cientes de asimetrı́a y curtosis multivariados mediante:
h 3 i
0 −1
β1,p = E (X 1 − µ) Σ (X 2 − µ)
h 2 i
β2,p = E (X 1 − µ)0 Σ−1 (X 1 − µ)

ya que los terceros momentos centrales para la distribución normal multivariada son cero se
tiene que β1,p = 0 cuando X es Np (µ; Σ) también se puede demostrar que para X normal
multivariado
β2,p = p(p + 2)
si definimos
b −1 (xj − x)
gij = (xi − x)0 Σ

donde Σ b −1 es el estimador de máxima verosimilitud para Σ−1 . Entonces los estimadores de

β1,p y β2,p son
n n
1 XX 3
b1,p = 2 g
n i=1 j=1 ij
n
1X 2
b2,p = g
n i=1 ii

Para decidir si rechazamos o no la hipótesis de multinormalidad se comparan b1,p y b2,p con

los valores tabulados (ver tabla A.5 de Rencher). Si n > 50 y con p ≥ 5, se tiene que

(p + 1)(n + 1)(n + 3)
z1 = b1,p
6[(n + 1)(p + 1) − 6]

14
tiene aproximadamente una distribución χ2 con q = 61 p(p + 1)(p + 2) grados de libertad.
Rechazamos la hipótesis de normalidad debido a la falta de simetrı́a si z1 ≥ χ2q (α)
Para b2,p se quiere rechazar para valores grandes (distribución demasiado aguda) o para
valores pequeños (distribución demasiado plana). Para el primer caso rechazamos si
b2,p − p(p + 2)
z2 = p , (1)
8p(p + 2)/n
que tiene aproximadamente distribución normal estándar, es mayor que zα/2 . En el segundo
caso (valores pequeños de b2,p ) se presentan dos situaciones:
ˆ Cuando 50 ≤ n ≤ 400 usamos
b2,p − p(p + 2)(n + p + 1)/n
z3 = p
8p(p + 2)/(n − 1)
es aproximadamente normal estándar, por lo tanto rechazamos si z3 < z1−α/2
ˆ Cuando n ≥ 400 usamos z2 dada por (1) y por tanto rechazamos si z2 < z1−α/2 .
A continuación se presenta la versión para R de la prueba de Mardia implementada en
SAS/IML por Dı́az Monroy (2002) y presentada en la sección 2.10, página 76.

Y<-as.matrix(Y) # convierte Y data frame a matriz

n<-nrow(Y) # numero de filas de Y
p<-ncol(Y) # numero de columnas de Y
gl_chi<-p*(p+1)*(p+2)/6 # grados de libertad
Q<-diag(n)-(1/n)*matrix(1,n,n) # I_p-(1/n)1_n1'_n
S<-(1/n)*t(Y)%*%Q%*%Y # matriz de covarianzas muestral
G_MATRIZ<- Q%*%Y%*%solve(S)%*%t(Y)%*%Q # Matriz g_hi de la ecuación 2.12
b_1<-sum(G_MATRIZ^3)/(n^2) # cálculo de la simetı́a
b_2<-sum(diag(G_MATRIZ^2))/n # calculo de la curtosis b_(2,p)
EST_b_1<-n*b_1/6 # calculo de la estadı́stica B1 ec. (2.13a)
EST_b_2<-(b_2-p*(p+2))/sqrt(8*p*(p+2)/n) # calculo de la estadı́stica B2
PVAL_ses<-1-pchisq(EST_b_1,gl_chi)
PVAL_cur<-2*(1-pnorm(abs(EST_b_2)))
cat("b_1=",b_1,"b_2=",b_2,"EST_b_1=",EST_b_1,"EST_b_2=",EST_b_2,"\n")
cat("PVAL_ses=",PVAL_ses,"PVAL_cur=",PVAL_cur,"\n")

La salida de la rutina se muestra a continuación.

b_1= 4.476382 b_2= 22.95687 EST_b_1= 20.88978 EST_b_2= -0.3983518
PVAL_ses= 0.4036454 PVAL_cur= 0.6903709

Los valores p indican que no se rechaza la hipótesis de normalidad, sin embargo estos valores
hay que tomarlos con cautela porque, según se dijo en la parte teórica, estos valores son
válidos para n ≥ 50.
Otra opción para rea,lizar la prueba de Mardia en R es la función mvm, peteneciente a la
librerı́a MVM.

15
library(MVN)
mvn(Y)

Prueba de Shapiro–Wilk multivariada

Una generalización multivariada de la prueba de Shapiro–Wilk, define zi = c0 y i , i = 1, 2, · · · , n
donde c es un vector de constantes y
n
ai (z(i) − z)2
P

W (c) = i=1Pn
(zi − z)2
i=1

donde z(1) ≤ z(2) ≤ · · · ≤ z(n) son los valores ordenados de z1 , z2 , · · · zn y ai son los coeficientes
que se encuentran en tablas. La hipótesis de normalidad multivariada no es rechazada si
max{W (c)} ≥ k
c

donde k = α es el nivel de significancia deseado.

El siguiente código de R realiza la prueba de Shapiro–Wilk multivariada, la librerı́a mvnormtest
no se encuentra en el paquete básico por tanto debe instalarse previamente.
library(mvnormtest)
U<-as.matrix(Y)
mshapiro.test(t(U))

La salida se muestra a continuación

Shapiro-Wilk normality test

data: Z
W = 0.9123, p-value = 0.02251

El p−valor bajo (< 0.05) indica que se rechaza la hipótesis nula H0 : los datos son normales

Estadı́stica E
La prueba E para normalidad multivariada fue propuesta e implementada por Szekely y Rizzo
(2005). La estadistica de prueba para normalidad p−variada está dada por
n n X n
!
2X 1 X
E =n Ekyi − Zk − EkZ − Z 0 k − 2 kyi − yj k ,
n i=1 n i=1 j=1

donde y1 , . . . , yn es la muestra estandarizada, Z, Z 0 son independientes e idénticamente dis-

tribuidos normales p−variantes estandarizados, y k · k denota la norma euclidiana.
La prueba E está implementada por bootstrap paramétrico con 999 réplicas (opción por
defecto)

16
 library(energy)
mvnorm.etest(Y)

Energy test of multivariate normality: estimated parameters

data: x, sample size 28, dimension 4, replicates 999

E-statistic = 1.2766, p-value = 0.007007

Nuevamente el p−valor indica que se rechaza la hipótesis nula, los datos no parecen provenir
de una distribución normal. Otra forma de realizar esta prueba es mediante el código

library(MVN)
mvn(Y, mvnTest = c("energy") )

17
Bibliography

Dı́az, L. G. (2002), Estadı́stica Multivariada: Inferencia y Métodos, Universidad Nacional de

Colombia.

Gross, J. (2003), The nortest Package.

Jarek, S. (2005), The mvnormtest Package.

Johnson, D. (2000), Métodos multivariados aplicados al análisis de datos, John Wiley & Sons.

Komsta, L. (2005), The moments Package.

Murrell, P. (2006), R Graphics, Chapman & Hall / CRC Sons.

N., R. & Dipak, D. (2002), A first course in linear models theory, Chapman & Hall / CRC.

Peña, D. (2000), Análisis de datos multivariantes, Thomson editores.

Rencher, A. (1995), Methods of Multivariate Analysis, John Wiley & Sons.

Rizzo, M. (2006), The energy Package.

Team, R. D. C. (2006), R: A Language and Environment for Statistical Computing, R Foun-

dation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0.

18