5.8 Sull’ipotesi di gaussianità

Abbiamo visto nelle due precedenti sezioni come l’introduzione di opportune variabili gaussiane e la successiva applicazione della formula di Bayes o la stima di massima verosimiglianza permetta di intepretare metodi come la PCA o i minimi quadrati in termini probabilistici. Oltre all’intepretazione, questo permette di chiarire le ipotesi sottostanti per garantire, almeno in teoria, una corretta applicazione del metodo. Ad esempio, nel caso dei minimi quadrati, (almeno) i residui devono avere densità gaussiane.

Rimane tuttavia il problema di argomentare che \(n\) osservazioni \((x_i)_{i=1}^n\) di dati dalla realtà possano essere ragionevolmente modellizate tramite eventi del tipo \(X_i =x_i\), dove le variabili aleatorie \(X_i\) siano gaussiane indipendenti, tutte con gli stessi parametri¹⁴.

Vi sono diversi approcci, ma si possono riassumere in essenzialmente due categorie: 1. Approcci qualitativi: si sfruttano teoremi limite, come la legge dei grandi numeri (ossia validi per la numerosità \(n\) molto grande) i quali garantiscono che determinate variabili empiriche, ossia dipendenti dalle osservazioni \((x_i)_{i=1}^n\) siano vicine in un senso opportuno, a variabili gaussiane di opportuni parametri. In questa categoria rientrano il confronto tra l’istogramma dei valori osservati a cui si sovrapponga la densità gaussiana stimata.

# Vediamo un esempio in cui è
# ragionevole l'ipotesi di gaussianità
# (i residui del modello lineare
# applicato al dataset Iris) e
# confrontiamolo con un esempio in cui
# invece l'ipotesi non lo sia (le
# osservazioni della lunghezza dei
# petali).

par(mfrow = c(1, 2))

iris_res <- iris_reg_lin$residuals
hist(iris_res, col = miei_colori[1], freq = FALSE,
  xlab = "valori dei residui", main = "")

# aggiungiamo il grafico della densità
# gaussiana con i parametri stimati
# dalle osservazioni (media e
# deviazione standard campionarie)

z <- seq(min(iris_res), max(iris_res), by = 0.01)
norm_dens <- dnorm(z, mean = mean(iris_res),
  sd = sd(iris_res))
lines(z, norm_dens, lwd = 3, col = miei_colori[2])


## vediamo invece lo stesso con la
## lughezza dei petali

iris_petali <- iris$Petal.Length
hist(iris_petali, col = miei_colori[3], freq = FALSE,
  xlab = "lunghezza dei petali", main = "")
z <- seq(min(iris_petali), max(iris_petali),
  by = 0.01)
norm_dens <- dnorm(z, mean = mean(iris_petali),
  sd = sd(iris_petali))
lines(z, norm_dens, lwd = 3, col = miei_colori[2])

# notiamo una buona aderenza dei dati
# alla densità teorica nel primo caso,
# e invece una notevole differenza nel
# secondo.

Un secondo metodo grafico è il Q-Q plot, in cui si confrontano la funzione quantile della variabile empirica, ossia la variabile uniforme discreta sugli \(n\) valori osservati, con la funzione quantile di una opportuna gaussiana – in tal caso, l’ipotesi di gaussianità è tanto più ragionevole quanto più i punti sul grafico siano allineati.

# per il Q-Q plot in R usiamo il
# comando qqnorm() che automaticamente
# confronta il quantile della variabile
# 'empirica' con quello di una
# gaussiana

par(mfrow = c(1, 2))

qqnorm(iris_res, col = miei_colori[1], pch = 16)

# per aggiungere la linea che si
# dovrebbe ottenere se l'ipotesi fosse
# vera (per opportuni parametri) usiamo
# il comando qqline()

qqline(iris_res, col = miei_colori[2], lwd = 2)


# ripetiamo per la lunghezza dei petali


qqnorm(iris_petali, col = miei_colori[3],
  pch = 16)

qqline(iris_petali, col = miei_colori[2],
  lwd = 2)

# nel primo caso la maggior parte dei
# punti è ben allineata alla retta,
# mentre nel secondo si discostano.

2. Approcci quantitativi: si introducono tuttavia dei test statistici in cui l’ipotesi nulla è l’evento (o meglio l’unione degli eventi al variare dei parametri di media e varianza) \(\mathcal{H}_0 =\) “le osservazioni provengono da \(n\) variabili gaussiane indipendenti, tutte con gli stessi parametri” e l’alternativa semplice \(H_1\) è la sua negazione. La descrizione specifica di questi test, in particolare dell’evento su cui si basa poi la decisione (rifiutare o meno) l’ipotesi nulla è troppo lunga, e la omettiamo. Ricordiamo però dal breve cenno ai test nella Sezione 2.7 che le quantità principali, in particolare il valore \(p\), sono calcolate rispetto alla probabilità condizionata alla validità dell’ipotesi nulla, pertanto si possono calcolare, almeno in linea di principio, perché l’ipotesi nulla riguarda densità gaussiane – il problema in questo caso sarebbe l’alternativa in cui non è chiaro quale densità considerare. In ogni caso, bisogna comunque prestare attenzione al fatto che un test statistico, pur essendo quantitativo non è una “dimostrazione” che l’ipotesi nulla sia vera (si usa appunto la locuzione non viene rifiutata per evitare di cadere in questa trappola concettuale). Ricordiamo infine che, più piccolo il valore \(p\), maggiore sarà il “grado di fiducia” che il test attribuisce nel rifiutare l’ipotesi nulla, quindi se il motivo per cui utilizziamo un test è di confermare, o meglio non smentire, l’ipotesi di gaussianità, il test sarà tanto più utile quanto più grande (ossia vicino ad \(1\)) è il valore \(p\). Vediamo degli esempi.

# Un test di gaussianità è dovuto a
# Shapiro (e Wilk).

shapiro.test(iris_res)

## 
##  Shapiro-Wilk normality test
## 
## data:  iris_res
## W = 0.99298, p-value = 0.6767

shapiro.test(iris_petali)

## 
##  Shapiro-Wilk normality test
## 
## data:  iris_petali
## W = 0.87627, p-value = 7.412e-10

# vediamo come nel caso dei residui il
# p-value sia molto grande (0.67),
# mentre nel secondo il test rifiuta la
# gaussianità con un p-value
# estremamente piccolo.

# Un test basato sulla funzione di
# ripartizione è dovuto a Kolmogorov e
# Smirnov. In realtà questo è un test
# che permette di confrontare anche con
# altre densità, non solo gaussiane,
# quindi dobbiamo specificare il
# parametro 'pnorm' per il test di
# gaussianità (con altri parametri si
# può testare altre densità).

ks.test(iris_res, "pnorm", mean = mean(iris_res),
  sd = sd(iris_res))

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  iris_res
## D = 0.040916, p-value = 0.9632
## alternative hypothesis: two-sided

ks.test(iris_petali, "pnorm", mean = mean(iris_petali),
  sd = sd(iris_petali))

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  iris_petali
## D = 0.19815, p-value = 1.532e-05
## alternative hypothesis: two-sided

# anche in questo caso il valore p è
# estremamente indicativo e conferma
# quanto osservato qualitativamente.

5.8.1 Esercizi

Esercizio 5.13 Si considerino le varie colonne del dataset ‘mtcars’ e si discuta se sia opportuno supporre che siano osservazioni di variabili gaussiane indipendenti. Lo stesso per i residui della regressione lineare della colonna ‘qsec’ rispetto alla colonna ‘hp’.