Capitolo 5 Stime ed incertezza

Nel capitolo precedente abbiamo visto che:

  1. I fenomeni biologici seguono una legge di natura (verità ‘vera’), che ne costituisce il meccanismo deterministico fondamentale. Questa legge di natura produce un risultato atteso \(Y_E\).
  2. Quando si organizza un esperimento, i soggetti sperimentali obbediscono a questo meccanismo di fondo, al quale tuttavia si sovrappongono molto altri elementi di ‘confusione,’ altamente incontrollabili, che vanno sotto il nome di errore sperimentale.
  3. L’osservazione sperimentale è quindi un’immagine confusa della verità vera (\(Y_O \neq Y_E\)) e, soprattutto, essa tende ad essere diversa per ogni sforzo di campionamento.
  4. Compito del ricercatore è comprendere come sia la verità ‘vera,’ separata dal ‘rumore di fondo’ generato dall’errore sperimentale.

Questo dualismo tra verità ‘vera’ (inconoscibile) e verità sperimentale (esplorabile tramite un esperimento opportunamente pianificato) è l’aspetto centrale di tutta la biometria ed è schematizzato nella figura 5.1.

Osservazioni sperimentali e meccanismi perturbativi

Figure 5.1: Osservazioni sperimentali e meccanismi perturbativi

5.1 Esempio: una soluzione erbicida

Nel capitolo precedente, abbiamo introdotto un esempio relativo ad un pozzo inquinato da un erbicida a concentrazione pari a 120 \(mg/L\), che viene misurata tramite un gascromatografo. Questo strumento di misura, unitamente a tutte le altre fonti ignote di errore, produce un coefficiente di variabilità del 10% (corrispondente ad una deviazione standard di 12 \(mg/L\)). Abbiamo anche visto che, se immaginiamo di fare le analisi in triplicato, come usuale per questo tipo di lavori, i risultati di questo esperimento possono essere simulati ricorrendo ad un generatore di numeri casuali:

set.seed(1234)
Y <- rnorm(3, 120, 12)
Y
## [1] 105.5152 123.3292 133.0133

Ricordiamo che i numeri casuali, in quanto tali, dovrebbero essere diversi ogni volta che li estraiamo, anche se, utilizzando la funzione set.seed(), è possibile indurre l’algoritmo a produrre sempre gli stessi valori, in modo che i calcoli di questo capitolo siano riproducibili.

5.1.1 Analisi dei dati: stima dei parametri

A questo punto mettiamoci in una situazione reale e, di conseguenza, dimentichiamo di conoscere che la concentrazione ignota è pari a 120 \(mg/L\) e che \(\sigma = 12\). Ipotizziamo quindi che le nostre osservazioni sperimentali siano generate da un modello del tipo:

\[Y_i = \mu + \varepsilon_i\]

con:

\[\varepsilon_i \sim N(0, \sigma)\]

Nelle due equazioni sovrastanti, gli elementi incogniti sono \(\mu\) e \(\sigma\). Guardando il campione, le nostre migliori stime per queste due quantità, che chiameremo rispettivamente \(m\) ed \(s\), sono pari rispettivamente alla media e alla deviazione standard del campione.

m <- mean(Y)
s <- sd(Y)
m; s
## [1] 120.6192
## [1] 13.9479

Questo processo con il quale assegniamo alla popolazione le caratteristiche del campione prende il nome di stima puntuale dei parametri. Vediamo ancora una volta che l’osservazione sperimentale non coincide con la verità ‘vera’ (\(m \neq \mu\) e \(s \neq \sigma)\), ma non siamo molto distanti, considerando il 10% di variabilità dello strumento di analisi. Tuttavia, visto che dobbiamo trarre conclusioni che riguardano la popolazione e non il campione, è giustificato da parte nostra un atteggiamento prudenziale: prima di dire che la concentrazione erbicida nella soluzione è pari 120.6192187, dobbiamo chiederci: che cosa succederebbe se ripetessimo l’esperimento molte altre volte?

5.1.2 La ‘sampling distribution’

In questo caso l’esperimento è solo ‘elettronico’ e possiamo quindi ripeterlo un numero anche molto elevato di volte, seguendo questa procedura:

  1. Ripetiamo l’estrazione precedente per 100’000 volte (ripetiamo l’analisi chimica per 100’000 volte, sempre con tre repliche)
  2. Otteniamo 100’000 medie
  3. Calcoliamo la media delle medie e la deviazione standard delle medie
# Simulazione MONTE CARLO - Esempio 1
set.seed(1234)
result <- rep(0, 100000)
for (i in 1:100000){
  sample <- rnorm(3, 120, 12)
  result[i] <- mean(sample)
}
mean(result)
## [1] 120.0341
sd(result)
## [1] 6.939063

In sostanza, la simulazione Monte Carlo ci consente di fare quello che dovremmo sempre fare, cioè ripetere l’esperimento un numero di volte molto elevato, anche se finito (un numero infinito è chiaramente impossibile!). A questo punto abbiamo in mano una popolazione di medie, che viene detta sampling distribution, un ‘oggetto’ abbastanza ‘teorico,’ ma fondamentale per la statistica frequentista, perché caratterizza la variabilità dei risultati di un esperimento, e quindi la sua riproducibilità.

Notiamo che:

  1. La media delle medie è praticamente coincidente con \(\mu\), la verità ‘vera.’ Ciò conferma che l’unico modo di ottenere risultati totalmente precisi è ripetere infinite volte l’esperimento;
  2. La deviazione standard delle medie è pari a 6.939063. Questo valore prende il nome di errore standard della media (SEM).

Esploriamo meglio la sampling distribution. Con R possiamo provare a discretizzarla e a riportarla su di un grafico a barre (figura 5.2 ).

Sampling distribution empirica e teorica

Figure 5.2: Sampling distribution empirica e teorica

5.1.3 L’errore standard

La sampling distribution che abbiamo ottenuto con la simulazione Monte Carlo è puramente empirica. Sarebbe interessante capire con più esattezza se esista una funzione di densità che permetta di descriverla con esattezza. In effetti, il grafico precedente mostra che la sampling distribution assomiglia molto ad una distribuzione normale, con media pari a 120 e deviazione standard pari all’errore standard.

Formalmente, il problema si può risolvere grazie alla legge di propagazione degli errori, che stabilisce tre importanti elementi:

  1. Se ho due variabili normalmente distribuite e le sommo tra di loro, la variabile risultante è ancora normale. Se ho una variabile normalmente distribuita e la moltiplico per una costante, la variabile risultante è ancora normale.
  2. Per variabili indipendenti, la varianza della somma è uguale alla somma delle varianze.
  3. La varianza del prodotto di una variabile per una costante \(k\) è pari alla varianza della variabile originale moltiplicata per \(k^2\).

Consideriamo che, quando preleviamo alcuni individui da una popolazione, ognuno di essi porta con sé una sua componente di incertezza, che egli ‘eredita’ dalla popolazione di cui fa parte. In questo caso, la popolazione ha una varianza pari a \(12^2 = 144\) e quindi ognuno dei tre soggetti campionati eredita tale varianza. Quando calcolo la media di tre osservazioni, in prima battuta io le sommo. A questo punto, dato che si tratta di osservazioni indipendenti, la propagazione degli errori (punto 2) ci dice che la varianza della somma è uguale a \(144 \times 3 = 432\).

Dopo aver sommato, il calcolo della media richiede che il risultato venga diviso per 3. La legge di propagazione degli errori (punto 3) ci dice quindi che la varianza viene divisa per \(3^2 = 9\). Insomma la popolazione delle medie è normale (punto 1), ha media pari a 120 e varianza pari a \(432/9 = 48\) e, di conseguenza, deviazione standard pari a \(\sqrt{48} = 6.928\), cioè \(12/\sqrt{3}\). In generale, l’errore standard di una media è:

\[\sigma_m = \frac{\sigma }{\sqrt n }\]

dove n è la dimensione del campione.

5.2 Stima per intervallo

Che cosa ci insegna questo esperimento? Ci insegna che, se prendiamo una distribuzione normale con media \(\mu\) e deviazione standard \(\sigma\) e cominciamo ad estrarre campioni, le medie dei campioni sono variabili, secondo una distribuzione normale con media \(\mu\) e deviazione standard \(\sigma / \sqrt{n}\). Questo concetto è interessante e può essere utilizzato per caratterizzare l’incertezza dei risultati di un esperimento. Riassumiamo:

  1. Abbiamo fatto un esperimento con tre repliche campionando da una distribuzione normale incognita.
  2. Abbiamo ottenuto i tre valori 105.5152, 123.3292 e 133.0133.
  3. In base alle osservazioni in nostro possesso, m = 120.6192 \(mg/L\), e, considerando che la cosa più probabile è che la media del campione sia uguale a quella della popolazione, concludiamo che \(\mu = m\).
  4. Dobbiamo adottare un atteggiamento prudenziale in relazione alla media, dato che non sappiamo il valore vero di \(\mu\). Sappiamo che le medie campionarie producono una sampling distribution caratterizzata da una deviazione standard pari a \(\sigma / \sqrt{n}\). Non conoscendo \(\sigma\), utilizziamo la sua miglior stima \(s\) e concludiamo che \(ES = 13.95/\sqrt{3} = 8.05\).
  5. Concludiamo quindi che \(\mu\) è pari a 120.6192 \(\pm\) 8.053.

Abbiamo caratterizzato l’incertezza del risultato attraverso un intervallo di valori (stima per intervallo).

5.3 L’intervallo di confidenza

La stima per intervallo fu uno degli interessi di ricerca del matematico polacco Jerzy Neyman (1894-1981), che definì la teoria degli intervalli di confidenza, ancora molto seguita anche ai giorni nostri. Partendo dal presupposto che le medie campionarie sono distribuite normalmente con media \(\mu\) e deviazione standard \(\sigma/\sqrt{n}\), è possibile calcolare la probabilità \(\alpha\) di trovare un campione la cui media era contenuta in un certo intervallo. Prendiamo, ad esempio, la nostra popolazione iniziale, con \(\mu\) = 120 e \(\sigma\) = 12. È facile vedere che c’è il 68% circa di probabilità di ottenere un campione con media inclusa nell’intervallo \(120 \pm 12a/\sqrt{3}\):

pnorm(120 + 12/sqrt(3), 120, 12/sqrt(3)) - 
  pnorm(120 - 12/sqrt(3), 120, 12/sqrt(3))
## [1] 0.6826895

Aumentando l’ampiezza dell’intervallo, è possibile aumentare la probabilità di campionare al suo interno. Ad esempio, se moltiplichiamo l’errore standard per due, la probabilità di ottenere un campione con una media inclusa nell’intervallo \(120 \pm 2 \times 12/\sqrt{3}\) supera (di poco) il 95%:

mult <- 2
pnorm(120 + mult * 12/sqrt(3), 120, 12/sqrt(3)) - 
  pnorm(120 - mult * 12/sqrt(3), 120, 12/sqrt(3))
## [1] 0.9544997

È possibile ottenere esattamente una probabilità del 95% utilizzando come moltiplicatore dell’errore standard il 97.5-esimo quantile della distribuzione normale standardizzata:1

mult <- qnorm(0.975, mean = 0, sd = 1)
pnorm(120 + mult * 12/sqrt(3), 120, 12/sqrt(3)) - 
  pnorm(120 - mult * 12/sqrt(3), 120, 12/sqrt(3))
## [1] 0.95

Se approssimiamo il 97.5-esimo quantile della distribuzione normale standardizzata alla seconda cifra decimale (per semplicità), possiamo scrivere:

\[ P \left[ \mu - 1.96 \times \frac{\sigma}{\sqrt{n} } \leq m \leq \mu + 1.96 \times \frac{\sigma}{\sqrt{n} } \right] = 0.95 \]

Per il nostro esempio possiamo concludere che la probabilità di estrarre un campione di acqua dal nostro pozzo inquinato, analizzarlo in triplicato ed ottenere una concentrazione media all’interno dell’intervallo \(\mu \pm 1.96 \times\sigma/\sqrt{n}\) è del 95%.

Con semplici passaggi algebrici, possiamo ottenere l’intervallo di confidenza:

\[ P \left[ m - 1.96 \times \frac{\sigma}{\sqrt{n} } \leq \mu \leq m + 1.96 \times \frac{\sigma}{\sqrt{n} } \right] = 0.95 \]

Proviamo a leggere l’espressione sovrastante: se prendiamo la stima \(m\) e il suo errore standard \(\sigma_m\) e calcoliamo un intervallo di confidenza come 1.96 volte \(\sigma_m\), esiste una probabilità pari a 0.95 che questo intervallo contenga \(\mu\), cioè la media vera e ignota della popolazione. In pratica, l’intervallo di confidenza potrebbe essere approssimato con il doppio dell’errore standard.

Il problema è che, nella pratica sperimentale, \(\sigma\) non è noto. Neyman propose di utilizzare \(s\) al posto di \(\sigma\), cioè la deviazione standard del campione, invece che quella della popolazione. Quindi, nel nostro caso, i limiti dell’intervallo di confidenza sono pari a:

m - qnorm(0.975) * s/sqrt(3)
## [1] 104.836
m + qnorm(0.975) * s/sqrt(3)
## [1] 136.4025

più semplicemente:

m - 2 * s/sqrt(3)
## [1] 104.5136
m + 2 * s/sqrt(3)
## [1] 136.7249

Questo che abbiamo calcolato è l’intervallo di confidenza per una probabilità del 95%. Aumentando opportunamente il moltiplicatore possiamo calcolare gli intervalli di confidenza per una probabilità del 99%, del 99.9% e così via. Di fatto, l’intervallo di confidenza del 95% è il più utilizzato in pratica.

5.4 Qual è il senso dell’intervallo di confidenza?

E’ utile ricordare il nostro punto di partenza e il nostro punto di arrivo:

  1. PUNTO DI PARTENZA: una distribuzione normale con \(\mu\) = 120 e \(\sigma\) = 12. Nella realtà assumiamo che la distribuzione di partenza sia normale, mentre i suoi parametri sono totalmente ignoti.
  2. PUNTO DI ARRIVO: una stima di \(\mu\) ed un intervallo di confidenza.

Che cosa significa questo intervallo? Esso fornisce:

  1. una misura di precisione: più piccolo è l’intervallo, maggiore è la precisione della stima;
  2. un’espressione di confidenza nel fatto che, se ripetessimo molte volte l’esperimento, nel 95% dei casi l’intervallo calcolato conterrebbe \(\mu\).

Insomma, l’intervallo di confidenza serve ad esplicitare la nostra incertezza sulla media vera della popolazione in studio.

5.5 Come presentare i risultati degli esperimenti

Dopo aver letto questo capitolo e quelli precedenti, dovrebbe essere chiaro che la presenza dell’errore sperimentale crea incertezza in relazione alle caratteristiche della popolazione da cui abbiamo estratto il campione. Pertanto, è sempre obbligatorio associare alle nostre stime un indicatore di incertezza, la cui assenza non è, in linea di principio, accettabile. Possiamo considerare le seguenti possibilità:

  1. riportare la media associata alla deviazione standard, per descrivere la variabilità originale del fenomeno in studio;
  2. riportare la media associata all’errore standard, per descrivere l’incertezza associata alla stima della media;
  3. riportare l’intervallo di confidenza ottenuto sottraendo/aggiungendo alla media il doppio dell’errore standard, per descrivere l’incertezza associata alla stima della media;

5.6 Alcune precisazioni

5.6.1 Campioni numerosi e non

Calcolare l’intervallo di confidenza utilizzando il doppio dell’errore standard costituisce un’approssimazione che è valida solo quando abbiamo esperimenti con un numero elevato di soggetti (maggiore di 15-20 circa). Per gli esperimenti piccoli, è necessario incrementare opportunamente il moltiplicatore. A questo fine viene utilizzato il valore corrispondente al 97.5-esimo percentile della distribuzione t di Student (non alla normale standardizzata, come abbiamo fatto finora), con un numero di gradi di libertà pari a quelli del campione studiato. Con tre soggetti il moltiplicatore sarebbe:

qt(0.975, 2)
## [1] 4.302653

Il moltiplicatore diminuisce all’aumentare della numerosità e, con 20 soggetti, diviene molto vicino a 2.

qt(0.975, 20)
## [1] 2.085963

Nel nostro esempio, con tre soggetti, l’intervallo di confidenza sarebbe:

m - qt(0.975, 2) * s/sqrt(3)
## [1] 85.97071
m + qt(0.975, 2) * s/sqrt(3)
## [1] 155.2677

Quindi, ben più alto di quello approssimato calcolato in precedenza. Chi ne volesse sapere di più, trova ulteriori informazioni in appendice.

5.6.2 Popolazioni gaussiane e non

In questo esempio siamo partiti da una popolazione con distribuzione gaussiana. In altri casi potrebbe non essere così. Ad esempio, immaginiamo di avere una popolazione di insetti, nella quale il rapporto tra maschi e femmine è ignoto. Campioniamo 40 insetti e contiamo 15 femmine. Qual è la proporzione di femmine nella popolazione originaria?

In questo caso stiamo studiando una grandezza che, almeno nel principio, non può essere gaussiana, ma è binomiale (vedi il capitolo precedente). Nonostante questo, possiamo utilizzare la stessa tecnica per la stima dell’intervallo di confidenza: sappiamo che la media di una distribuzione binomiale è \(p = 14/40 = 0.375\), mentre la deviazione standard è \(\sigma = \sqrt{0.375 \times (1 - 0.375)} = 0.484\). Di conseguenza, l’errore standard è \(0.484 / \sqrt{40} = 0.077\). L’intervallo di confidenza sarà dato quindi da:

0.375 - 2 * 0.077
## [1] 0.221
0.375 + 2 * 0.077
## [1] 0.529

Chi fosse interessato ad approfondire questi aspetti può proseguire nella lettura, dopo gli esercizi. Gli altri potranno fermarsi agli esercizi sottostanti.

5.7 Analisi statistica dei dati: riassunto del percorso logico

Considerando quanto finora detto, possiamo riassumere la logica dell’inferenza tradizionale nel modo seguente:

  1. Un esperimento è solo un campione di un numero infinito di esperimenti simili che avremmo potuto/dovuto eseguire, ma che non abbiamo eseguito, per mancanza di risorse;
  2. Assumiamo che i dati del nostro esperimento sono generati da un modello matematico probabilistico, che prende una certa forma algebrica e ne stimiamo i parametri utilizzando i dati osservati;
  3. Costruiamo la sampling distribution per i parametri stimati o per altre statistiche rilevanti, in modo da caratterizzare i risultati delle infinite repliche del nostro esperimento, che avremmo dovuto fare, ma che non abbiamo fatto.
  4. Utilizziamo la sampling distribution per l’inferenza statistica.

5.8 Da ricordare

  1. La natura genera i dati
  2. Noi scegliamo un modello deterministico che simula il meccanismo di generazione dei dati attuato dalla natura.
  3. Stimiamo i parametri.
  4. Confrontiamo le previsioni con i dati osservati. Determiniamo \(\epsilon\) e la sua deviazione standard (\(\sigma\))
  5. Assumiamo un modello stocastico ragionevole per spiegare \(\epsilon\), quasi sempre di tipo gaussiano, con media 0 e deviazione standard pari a \(\sigma\), indipendente dalla X (omoscedasticità)
  6. Qualunque stima sperimentale deve essere associata ad un indicatore di variabilità (errore standard o intervallo di confidenza).

5.9 Per approfondire un po’…

5.10 Coverage degli intervalli di confidenza

Abbiamo visto che un metodo semplice per costruire un intervallo di confidenza è utilizzare il doppio dell’errore standard. Questo intervallo, se viene utilizzato come misura di precisione/incertezza, è sempre accettabile. Tuttavia, da un punto di vista strettamente probabilistico, è lecito chiedersi: ma è proprio vero che se io ripeto l’esperimento molte volte e calcolo sempre l’intervallo di confidenza, riesco a centrare la media \(\mu\) nel 95% dei casi? È bene sapere che, con termine inglese, l’effettiva percentuale di campioni per i quali l’intervallo di confidenza, calcolato per un certo P nominale (es. P = 0.95), contiene effettivamente la media \(\mu\) della popolazione, viene detta coverage. Quindi la nostra domanda è: qual è il coverage dell’intervallo di confidenza calcolato con il doppio dell’errore standard?

Proviamo a rispondere a questa domanda con una simulazione Monte Carlo. Prendiamo la solita popolazione normalmente distribuita con \(\mu = 120\) e \(\sigma = 12\) ed estraiamo centomila campioni. Per ogni campione, calcoliamo l’intervallo di confidenza della media (P = 0.95) considerando il doppio dell’errore standard. Verifichiamo poi se questo intervallo contiene il valore 120: se si, assegniamo al campionamento il valore 1 (successo), altrimenti assegniamo il valore 0.

result <- rep(0, 100000)
set.seed(1234)
for (i in 1:100000){
  sample <- rnorm(3, 120, 12)
  limInf<- mean(sample) - sd(sample)/sqrt(3) * 2 
  limSup<- mean(sample) + sd(sample)/sqrt(3) * 2
  if (limInf<= 120 & limSup>= 120) result[i] = 1
}
sum(result)/100000
## [1] 0.81708

La simulazione mostra che la risposta alla domanda precedente è no: il nostro intervallo di confidenza non è riuscito a centrare la media nel 95% dei casi; ciò è avvenuto in poco più dell’80% dei casi (coverage dell 81.7%, circa). In realtà, possiamo facilmente verificare, con altre simulazioni di Monte Carlo, che il coverage si avvicina al 95% solo se abbiamo campioni di numerosità superiore a 15-20 circa.

result <- rep(0, 100000)
set.seed(1234)
for (i in 1:100000){
  n <- 15
  sample <- rnorm(n, 120, 12)
  limInf<- mean(sample) - sd(sample)/sqrt(n) * 2 
  limSup<- mean(sample) + sd(sample)/sqrt(n) * 2
  if (limInf<= 120 & limSup>= 120) result[i] = 1
}
sum(result)/100000
## [1] 0.93594

Insomma, quando gli esperimenti sono piccoli, con poche repliche, dovremmo trovare un metodo di calcolo degli intervalli di confidenza un po’ più affidabile, se veramente volessimo ottenere un coverage pari a quello nominale (P = 0.95).

Il problema, già accennato, nasce dal fatto che \(\sigma_m\) viene sostituito con \(s_m\), cioè il valore di errore standard stimato nel campione. Come tutte le stime, anche \(s_m\) è soggetto ad incertezza, il che aggiunge un elemento ulteriore di imprecisione nella sampling distribution di \(T\). Insomma ci chiediamo, la sampling distribution di T, calcolata con \(s\) invece che \(\sigma\) è ancora normale? Verifichiamo questo aspetto empiricamente, con una nuova simulazione Monte Carlo. Questa volta facciamo la seguente operazione:

  1. campioniamo tre individui
  2. Calcoliamo il valore di \(T\) con la statistica precedente, utilizzando la deviazione standard del campione e lo salviamo
  3. Con un po’ di pazienza, ripetiamo il tutto 100’000 volte.
#SIMULAZIONE MONTE CARLO - t di Student
set.seed(435)
result <- c()
for (i in 1:100000){
  sample3 <- rnorm(3, 120, 12)
  Ti <- (mean(sample3) - 120) / (sd(sample3)/sqrt(3))
  result[i] <- Ti
  }

Se riportiamo i valori ottenuti su una distribuzione di frequenze otteniamo il grafico in Figura 5.3.

Sampling distribution empirica per le medie campionarie, insieme ad una distribuzione gaussiana (blue) e t di Student con 2 gradi di libertà (rossa)

Figure 5.3: Sampling distribution empirica per le medie campionarie, insieme ad una distribuzione gaussiana (blue) e t di Student con 2 gradi di libertà (rossa)

Vediamo che la sampling distribution di \(T\), determinata utilizzando \(s\) invece che \(\sigma\) è solo approssimativamente normale. E’ facile vedere che questa approssimazione è sufficientemente buona solo se la numerosità del campione diviene abbastanza grande (es. \(n > 30)\), ma non certamente quando \(n\) = 3. In questo caso, la sampling distribution che osserviamo è più ‘dispersa’ di quella normale, con un maggior numero di valori sulle code.

Neyman scoprì che la sampling distribution di T poteva essere perfettamente descritta utilizzando la distribuzione t di Student, con un numero di gradi di libertà pari a quelli del campione (in questo caso 2), come vediamo nella figura 5.3. In realtà questa conclusione era stata già raggiunta da William Sealy Gosset (1876 - 1937), uno statistico impiegato presso la fabbrica londinese della famosa birra Guinness, dove elaborava i dati relativi all’andamento del processo di maltazione. Egli, avendo definito questa nuova funzione di densità, per aggirare il divieto di pubblicazione imposto dal suo datore di lavoro, pubblicò i risultati sotto lo pseudonimo Student, da cui deriva il nome della distribuzione di densità.

Quindi, quando i campioni sono piccoli, il modo giusto di calcolare l’intervallo di confidenza è quello di utilizzare l’espressione seguente:

\[P \left( m - \textrm{qt}(0.975,n - 1) \cdot s_m \le \mu \le m + \textrm{qt}(0.975,n - 1) \cdot s_m \right) = 0.95\]

dove \(\textrm{qt}(0.025,n - 1)\) e \(\textrm{qt}(0.975,n - 1)\) sono rispettivamente il 2.5-esimo e il 97.5-esimo percentile della distribuzione t di Student, con n-1 gradi di libertà.

E’ facile osservare che, se l’intervallo di confidenza è calcolato in questo modo, il suo coverage effettivo e pari al 95%.

result <- rep(0, 100000)
set.seed(1234)
for (i in 1:100000){
  sample <- rnorm(3, 120, 12)
  limInf<- mean(sample) + sd(sample)/sqrt(3) * qt(0.025, 2) 
  limSup<- mean(sample) + sd(sample)/sqrt(3) * qt(0.975, 2) 
  if (limInf<= 120 & limSup>= 120) result[i] = 1
}
sum(result)/100000
## [1] 0.94936

Tuttavia, la formula di Neyman, anche se assicura un coverage pari a quello nominale, si presta a cattive letture, che sono insensate da un punto di vista probabilistico, ma tuttavia molto frequenti nella pratica operativa. Ad esempio, immaginiamo di aver effettuato un campionamento dalla nostra popolazione con \(\mu = 120\) e \(\sigma = 12\). Il risultato è:

set.seed(1234)
x <- rnorm(3, 120, 12)
m <- mean(x)
s <- sd(x)
m; s
## [1] 120.6192
## [1] 13.9479

Di conseguenza, l’intervallo di confidenza va da 85.9707117 a 155.2677256

  1. NON E’ VERA l’affermazione che: c’è il 95% di probabilità che la media ‘vera’ sia tra 85.0 e 155.3. La media vera della popolazione è sempre fissa e pari a 120 e non cambia seguendo una distribuzione di probabilità. L’affermazione probabilistica deve essere riferita alla possibilità che l’intervallo di confidenza la centri, non al valore della media.
  2. E’ DUBBIA l’affermazione che: c’è il 95% di probabilità che l’intervallo di confidenza contenga \(\mu\). In questa affermazione, la probabilità è legata all’intervallo di confidenza, ma l’affermazione è ugualmente sbagliata se si riferisce al singolo campionamento. Infatti, il singolo e specifico intervallo di confidenza (86.0 - 155.3) può contenere o no \(\mu\), ma non abbiamo alcun elemento per sapere se effittavamento lo contiene, neanche in termini di probabilità. Invece l’affermazione è corretta se si riferisce al futuro, cioè alle ipotetiche repliche dell’esperimento.
  3. NON E’ VERA: l’affermazione che ripetendo infinite volte l’esperimento, il 95% delle stime che otteniamo cadono nell’intervallo 86.0 e 155.3. Una semplice simulazione mostra che tutte le medie campionate cadono in quell’intervallo
result <- rep(0, 100000)
set.seed(1234)
for (i in 1:100000){
  sample <- rnorm(3, 120, 12)
  if (mean(sample) <= 155.0 & mean(sample) >= 85.3) result[i] = 1
}
sum(result)/100000
## [1] 1

Insomma, l’intervallo di confidenza vale per la sampling distribution e non vale per ogni singolo campionamento (esperimento). L’unica cosa che possiamo lecitamente affermare è che, se ripetessimo l’esperimento un numero molto elevato di volte e calcolassimo l’intervallo di confidenza sempre con la formula di Neyman, nel 95% dei casi saremmo in grado di ‘catturare’ la media vera della popolazione all’interno del nostro intervallo di confidenza, che diviene in questo modo una sorta di polizza assicurativa, per contenere la nostra probabilità d’errore, nel lungo periodo, al disotto del 5%.

5.10.1 Intervalli di confidenza per fenomeni non-normali

Nel sottocapitolo precedente abbiamo presentato un esempio in cui avevamo campionato da una distribuzione normale, riscontrando una sampling distribution per la media campionaria anch’essa normale (almeno approssimativamente). Ma che succede se la distribuzione di partenza è non-normale? La sampling distribution di uno stimatore è ancora normale? Vediamo un nuovo esempio.

Immaginiamo di avere 4’000’000 di semi ben mischiati (in modo che non ci siano raggruppamenti non casuali di qualche tipo), che costituiscono la nostra popolazione di partenza. Vogliamo appurare la frequenza relativa (p) dei semi dormienti. Questa informazione, nella realtà, esiste (\(\pi\) = 0.25), ma non è nota.

Dato l’elevato numero di ‘soggetti,’ non possiamo testare la germinabilità di tutti i semi, ma dobbiamo necessariamente prelevare un campione casuale di 40 soggetti; ogni seme viene saggiato e, dato che la popolazione è molto numerosa, l’estrazione di un seme non modifica sensibilmente la proporzione di quelli dormienti nella popolazione (esperimenti indipendenti).

Dopo aver descritto la popolazione e l’esperimento, ci chiediamo quale sia il modello matematico che genera i nostri dati (numero di successi su 40 semi estratti). Il disegno sperimentale ci assicura che ogni estrazione è totalmente indipendente dalla precedente e dalla successiva ed ha due soli risultati possibili, cioè successo (seme dormiente), o insuccesso (seme germinabile). Di conseguenza, ogni singola estrazione si configura come un esperimento Bernoulliano, con probabilità di successo pari a \(\pi\), il cui valore ‘vero’ esiste, è fisso, pre-determinato (esiste ancor prima di organizzare l’esperimento), anche se incognito e inconoscibile, a meno di non voler/poter esaminare tutti i semi disponibili. L’insieme delle 40 estrazioni (40 esperimenti Bernoulliani) può produrre un ventaglio di risultati possibili, da 40 successi a 40 insuccessi, per un totale di 41 possibili ‘outcomes.’

E’ evidente che i 41 possibili risultati non sono ugualmente probabili e si può dimostrare che la probabilità di ottenere k successi (con k che va da 0 ad n; n è al numero delle estrazioni) dipende da \(\pi\) ed è descrivibile matematicamente con la distribuzione binomiale \(\phi\):

\[\phi(k, n, p) = \frac{n!}{(n-k)!k!} p^k (1 - p)^{(n-k)}\]

Abbiamo quindi definito il modello matematico che descrive la probabilità di tutti i possibili risultati del nostro esperimento e quindi può in qualche modo essere considerato il ‘meccanismo’ che ‘genera’ i dati sperimentali osservati. Si tratta di un meccanismo puramente ‘stocastico’ nel quale è solo il caso che, attraverso il campionamento, determina il risultato dell’esperimento.

Con queste informazioni, possiamo simulare un esperimento con R, ottenendo i seguenti risultati:

set.seed(236)
rbinom(1, 40, 0.25)
## [1] 9

Abbiamo ottenuto 9 successi su 40, cioè 9 semi dormienti su 40 saggiati.La proporzione osservata è p = 9/40 = 0.225. Concludiamo (stima puntuale) che \(\pi\) = 0.225. Anche in questo caso vi è chiara discrasia tra la verità ‘vera’ e l’osservazione sperimentale (tra \(\pi\) e \(p\)).

Cosa succede se ripetiamo l’esperimento? Come abbiamo imparato a fare, possiamo cercare una risposta attraverso la simulazione Monte Carlo, ricorrendo ad un generatore di numeri casuali da una distribuzione binomiale con n = 40 e \(\pi\) = 0.25 (in R si usa la funzione ‘rbinom(numeroDatiCasuali, n, p)’). Il codice è più semplice, in quanto non è necessario impostare un ciclo iterativo:

set.seed(1234)
result <- rbinom(10000000, 40, 0.25)

Esploriamo i risultati ottenuti:

result_p <- result/40
mean(result_p)
## [1] 0.2499795
sd(result_p)
## [1] 0.06849423

Osserviamo subito che, anche se i singoli esperimenti portano a stime diverse da \(\pi\) vero, la media di \(p\) tende ad essere uguale a \(\pi\). L’errore standard (deviazione standard della sampling distribution) è 0.0685. Fino a qui, non vi è nulla di diverso dall’esempio precedente, se teniamo presente che la deviazione standard della popolazione originale (che è binomiale) è pari a \(\sqrt{p \times (1 - p)}\), quindi l’errore standard è \(\sqrt{0.25 \times 0.75 / 40} = 0.0685\).

Rimane da stabilire se la sampling distribution di \(p\) è normale. Possiamo utilizzare i 10’000’000 di valori ottenuti per costruire una distribuzione empirica di frequenze, come nel codice sottostante.

breaks <- seq(0, 0.7, by=0.025)
freqAss <- as.numeric( table(cut(result_p, breaks) ) ) 
freqRel <- freqAss/length(result_p)
density <- freqRel/0.025
p_oss <- breaks[2:length(breaks)]

La distribuzione empirica della proporzione campionaria è visibile in Figura @rif(fig:figName2541).

plot(density ~ p_oss, type = "h",
     xlab = expression(paste(bar(p))),
     ylab="Density", 
    xlim=c(0,0.6) )
b <- seq(0, 1, by=0.1)
curve(dnorm(x, 0.25, 0.0685), add=TRUE, col="red")
Sampling distribution per le proporzioni campionarie

Figure 5.4: Sampling distribution per le proporzioni campionarie

Vediamo che sampling distribution è approssimativamente normale con media pari a 0.25 e deviazione standard pari a 0.0685. Lo percepiamo chiaramente dal grafico soprastante, ma c’è una spiegazione scientifica per questo, basata sul TEOREMA DEL LIMITE CENTRALE:

  1. La sampling distribution di una statistica ottenuta da campioni casuali e indipendenti è approssimativamente normale, indipendentemente dalla distribuzione della popolazione da cui i campioni sono stati estratti.
  2. La media della sampling distribution è uguale al valore della statistica calcolata sulla popolazione originale, la deviazione standard della sampling distribution (errore standard) è pari alla deviazione standard della popolazione originale divisa per la radice quadrata della numerosità di un campione.

5.11 Altre letture

  1. Hastie, T., Tibshirani, R., Friedman, J., 2009. The elements of statistical learning, Springer Series in Statistics. Springer Science + Business Media, California, USA.

  1. Tra tutte le distribuzioni normali, ce n’è una particolare, che ha media 0 e deviazione standard 1. Questa si chiama distribuzione normale standardizzata↩︎