I test di verifica di normalità

Numerose volte nei posts precedenti abbiamo proceduto con l'analisi assumendo che i dati seguissero una distribuzione normale (applicando quindi la statistica parametrica). Vediamo adesso come verificare che i propri dati seguano una distribuzione gaussiana. I test a disposizione sono molto numerosi, e la maggior parte di questi sono disponibili nei vari pacchetti di R. Vediamo qui di seguito l'applicazione di alcuni di essi, ed in particolare:

• D'agostino's K-squared test


• Jarque-Bera test


• Anderson-Darling test for normality


• Cramer-von Mises normality test


• Lilliefors (Kolmogorov-Smirnov) test for normality


• Shapiro-Francia normality test


• Pearson chi-square normality test


• Shapiro-Wilk normality test


• Kolmogorov-Smirnov test

Per tutti gli esempi che analizzeremo, userò un set di dati ottenuto applicando la funzione rnorm(), la quale ci fornisce 50 numeri che seguono una distribuzione normale, nota la media (10) e la deviazione standard (2). Ci aspettiamo quindi che tutti i test che vedremo affermino che i dati seguono una distribuzione normale.

x <- rnorm(50, 10, 2)

---

D'agostino's K-squared test (D'Agostino-Pearson normality test)

Che io sappia questo test non è stato implementato in nessun pacchetto di R. Tuttavia, conoscendo le formule (disponibili su WikiPedia), è facilmente ricostruibile una funzione. Qui di seguito riporto il codice:

dagostino.pearson.test <- function(x) {
    # from Zar (1999), implemented by Doug Scofield, scofield at bio.indiana.edu
    DNAME <- deparse(substitute(x))
    n <- length(x)
    x2 <- x * x
    x3 <- x * x2
    x4 <- x * x3
    # compute Z_g1
    k3 <- ((n*sum(x3)) - (3*sum(x)*sum(x2)) + (2*(sum(x)^3)/n)) /
          ((n-1)*(n-2))
    g1 <- k3 / sqrt(var(x)^3)
    sqrtb1 <- ((n - 2)*g1) / sqrt(n*(n - 1))
    A <- sqrtb1 * sqrt(((n + 1)*(n + 3)) / (6*(n - 2)))
    B <- (3*(n*n + 27*n - 70)*(n+1)*(n+3)) / ((n-2)*(n+5)*(n+7)*(n+9))
    C <- sqrt(2*(B - 1)) - 1
    D <- sqrt(C)
    E <- 1 / sqrt(log(D))
    F <- A / sqrt(2/(C - 1))
    Zg1 <- E * log(F + sqrt(F*F + 1))
    # compute Z_g2
    G <- (24*n*(n-2)*(n-3)) / (((n+1)^2)*(n+3)*(n+5))
    k4 <- (((n*n*n + n*n)*sum(x4)) - (4*(n*n + n)*sum(x3)*sum(x)) -
          (3*(n*n - n)*sum(x2)^2) + (12*n*sum(x2)*sum(x)^2) - 
          (6*sum(x)^4)) /(n*(n-1)*(n-2)*(n-3))
    g2 <- k4 / var(x)^2
    H <- ((n-2)*(n-3)*abs(g2)) / ((n+1)*(n-1)*sqrt(G))
    J <- ((6*(n*n - 5*n + 2)) / ((n+7)*(n+9))) * sqrt((6*(n+3)*(n+5)) /(n*(n-2)*(n-3)))
    K <- 6 + (8/J)*(2/J + sqrt(1 + 4/(J*J)))
    L <- (1 - 2/K) / (1 + H*sqrt(2/(K-4)))
    Zg2 <- (1 - 2/(9*K) - (L^(1/3))) / (sqrt(2/(9*K)))
    K2 <- Zg1*Zg1 + Zg2*Zg2
    pk2 <- pchisq(K2, 2, lower.tail=FALSE)
    RVAL <- list(statistic = c(K2 = K2), p.value = pk2, method =
"D'Agostino-Pearson normality test\n\nK2 is distributed as Chi-squared
with df=2", alternative = "distribution is not normal", data.name =
DNAME)
    class(RVAL) <- "htest"
    return(RVAL)
}

L'ipotesi nulla è che i nostri dati seguono una distribuzione normale. Applichiamo la funzione ai nostri dati:

dagostino.pearson.test(x)

        D'Agostino-Pearson normality test
        
        K2 is distributed as Chi-squared with df=2

data:  x 
K2 = 1.2788, p-value = 0.5276
alternative hypothesis: distribution is not normal 

Essendo p-value > 0.05, accettiamo l'ipotesi nulla: siamo di fronte a una distribuzione normale.

---

Jarque-Bera normality test

E' disponibile nel pacchetto tseries, che deve essere scaricato, ed installato.

library(tseries)
jarque.bera.test(x)

        Jarque Bera Test

data:  x 
X-squared = 1.1181, df = 2, p-value = 0.5717

Anche per questo test, l'ipotesi nulla è la normalità del vettore x. Essendo p-value > 0.05, accettiamo l'ipotesi nulla.

---

Anderson-Darling test for normality

Test disponibile nella libreria nortest (libreria che raccoglie numerosi test di normalità).

library(nortest)
ad.test(x)

        Anderson-Darling normality test

data:  x 
A = 0.1931, p-value = 0.8898

Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0

---

Cramer-von Mises normality test

library(nortest)
cvm.test(x)

        Cramer-von Mises normality test

data:  x 
W = 0.0289, p-value = 0.8576

Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0

---

Lilliefors (Kolmogorov-Smirnov) test for normality

library(nortest)
lillie.test(x)

        Lilliefors (Kolmogorov-Smirnov) normality test

data:  x 
D = 0.0808, p-value = 0.5715

Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0

---

Shapiro-Francia normality test

library(nortest)
sf.test(x)

        Shapiro-Francia normality test

data:  x 
W = 0.9918, p-value = 0.9454

Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0

---

Pearson chi-square test for normality

library(nortest)
pearson.test(x)

        Pearson chi-square normality test

data:  x 
P = 1.2, p-value = 0.991

Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0

---

Shapiro-Wilk normality test

Questo è uno dei test più utilizzati, ed è già disponibile nel pacchetto di base di R (stats), quindi non necessita di alcuna installazione.

shapiro.test(x)

        Shapiro-Wilk normality test

data:  x 
W = 0.9865, p-value = 0.834

Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0

---

Kolmogorov–Smirnov test

Anche questo test è disponibile nel pacchetto di base. Può essere utilizzato per effettuare due diverse analisi. Possiamo ad esempio chiederci "i dati X e Y provengono dalla stessa distribuzione?". In questo caso la verifica si effettua in questo modo:

x <- rnorm(50)
y <- rnorm(30)

ks.test(x, y)

        Two-sample Kolmogorov-Smirnov test

data:  x and y 
D = 0.1533, p-value = 0.7195
alternative hypothesis: two-sided 

Essendo p-value > 0.05, accettiamo l'ipotesi nulla: i due vettori provengono dalla stessa distribuzione.

Un'altra verifica che possiamo effettuare è la seguente: "il vettore X proviene da una distribuzione normale, o da una distribuzione chi-quadro a 1 grado di libertà?". In R si procede così:

ks.test(x, "pnorm")

        One-sample Kolmogorov-Smirnov test

data:  x 
D = 0.0931, p-value = 0.7436
alternative hypothesis: two-sided 

ks.test(x, "pchisq", 1)

        One-sample Kolmogorov-Smirnov test

data:  x 
D = 0.54, p-value = 3.564e-14
alternative hypothesis: two-sided 

In questo caso, oltre a specificare il vettore da analizzare, si deve pure specificare la funzione di confronto, messa tra virgolette. Per testare altre distribuzioni, potete trovare la sintassi precisa digitando ??distribution: così leggete l'elenco delle distribuzioni implementate in R; quindi, richiamando l'help relativo alla distribuzione di vostro interesse (ad esempio help(Normal)), potete vedere la corretta sintassi da utilizzare (generalmente è il nome della funzione, preceduto dalla lettera p: pnorm, pgamma, pchisq).
Come potevamo prevedere, il test di normalità dà p-value > 0.05, mentre il confronto con la distribuzione chi-quadro dà p-value < 0.05. Quindi il vettore x segue una distribuzione normale.