L’approccio simulativo

Gli errori di campionamento

L’analisi dei dati deve confrontarsi con la gestione degli errori. Se una buona metodologia ed un corretto campionamento possono minimizzare l’impatto degli errori sistematici (i bias), gli errori casuali non possono essere eliminati. L’analisi inferenziale permette al ricercatore di stimare l’entità di questi errori, e di capire quanto le misure e le relazioni emerse siano da imputare a tali errori.

L’analisi si basa sul calcolo di alcune statistiche. Nell’analisi univariata si calcolano gli indici di centralità (moda, mediana, media) e di dispersione (la differenza interquantile, la varianza); nelle statistiche bivariate si calcolano delle statistiche capaci di misurare le relazioni fra variabili. Le statistiche uni che bivariate devono tener conto dell’errore di campionamento.

Facciamo alcuni esempi.

Se il campionamento è stato fatto in maniera corretta la media del campione costituisce la migliore stima della media della popolazione (la media è una stima unbiased); se dalla stessa popolazione, però, estraggo dieci campioni diversi, otterrò dieci stime differenti.

Un tipico disegno sperimentale consiste nel dividere il campione in 2 gruppi, somministrare un trattamento ad un gruppo (sperimentale), somministrare un diverso trattamento (o un placebo) all’altro gruppo, e misurare l’effetto attraverso una variabile numerica; per valutare l’effetto del trattamento, si misura la differenza fra le medie dei due gruppi. Anche in questo caso ci si deve chiedere: questa differenza va attribuita al trattamento o al caso (alla variabilità campionaria)? Infatti, in maniera del tutto paragonabile all’esempio precedente, cosa succederebbe se applicassimo lo stesso trattamento (o nessun trattamento) ai due gruppi? Ci aspettiamo che le medie dei due gruppi siano perfettamente uguali? La risposta è naturalmente no: le medie saranno probabilmente simili, ma non uguali.

Facciamo un terzo esempio: immaginiamo di voler capire se vi è una relazione fra due variabili numeriche. Decidiamo di adottare la statistica della correlazione di Pearson, una misura che si muove nel range $-1 < r < +1$ e dove 0 significa assenza di correlazione. Anche in questo caso, però, nella circostanza di due variabili fra loro indipendenti, non possiamo aspettarci una correlazione esattamente pari a 0.

Distribuzione degli errori

Approccio parametrico

Fortunatamente, gli errori dovuti al caso (e alla varianza campionaria) sono soggetti a delle distribuzioni note (quantomeno per quanto riguarda le statistiche più comuni). La cosiddetta statistica parametrica si basa proprio sul fatto che, se alcuni assunti sono rispettati, la distribuzione dell’errore delle statistiche usate approssima, previo opportuna trasformazione, delle distribuzioni teoriche. Il processo inferenziale utilizza questa proprietà; si calcola la statistica, si opera la trasformazione, e si confronta il risultato con la distribuzione teorica.

Statistiche non parametriche

Il limite dell’approccio parametrico è che fa delle assunzioni sulle variabili; vi sono delle circostanze in cui queste assunzioni non vengono rispettate. In questi casi, le statistiche parametriche possono essere inaffidabili; a questo punto, diventa opportuno affidarsi a delle famiglie di statistiche non parametriche, il cui vantaggio è quello di fare un minore numero di assunzioni.

Generalmente, l’approccio delle statistiche non parametriche consiste nel trasformare la variabile dipendente, numerica, in una variabile ordinale. La trasformazione consiste nel calcolare il rank, ovvero il valore ordinale della misura.

Approccio simulativo (resampling)

Esiste poi un’altra possibilità: utilizzare il calcolatore per generare la distribuzione dell’errore, e basare il processo inferenziale non sulla distribuzione teorica, ma sulla distribuzione generata.

Questo approccio è relativamente recente, in quanto è computazionalmente oneroso, e dunque può essere applicato soltanto con degli strumenti di calcolo potenti. Oggi, però, possono essere applicati agevolmente anche con i comuni computer, e dunque stanno guadagnando crescente popolarità.

L’approccio simulativo ha alcuni vantaggi, il principale dei quali è che fa pochissime assunzioni, e dunque può essere applicato anche nel caso, ad esempio, di distribuzioni che non possono essere ricondotte alle distribuzioni teoriche. Il secondo vantaggio è che può essere applicato anche a statistiche non comuni, per le quali non esiste – o non è nota – una distribuzione teorica.

L’approccio simulativo ha infine il vantaggio di essere particolarmente intuitivo, in quanto permette di mostrare l’errore di campionamento, la sua distribuzione, e i rispettivi parametri. Questa caratteristica lo rende particolarmente indicato ai fini didattici, in quanto è possibile simulare la varianza di campionamento, generare la distribuzione campionaria, e confrontarla con la distribuzione teorica. L’approccio computazionale è inoltre un ottimo modo per giocare con strumenti come R, prendere confidenza con il linguaggio, e capire cosa succede dietro alle quinte quando si utilizzano le funzioni di testing – parametriche e non parametriche.

Introduzione all’approccio simulativo

Per introdurre l’approccio simulativo, utilizziamo R per fare delle simulazioni che ci permettano di riprodurre, in laboratorio, l’errore di campionamento.

Attraverso la simulazione possiamo creare delle circostanze difficilmente riproducibili nella realtà: possiamo generare una popolazione, generare un numero molto alto di campioni, e valutare qualitativamente (graficamente) e quantitativamente l’errore stocastico di campionamento.

Generare popolazione e campioni

Generare la popolazione

Nel contesto della simulazione, generare una popolazione significa generare un vettore di valori casuali. Se si assume che la distribuzione della popolazione sia normale, è possibile utilizzare la funzione rnorm per generare un vettore di numeri distribuiti normalmente intorno ad una media e con una deviazione standard predefinita.

La lunghezza del vettore corrisponde alla numerosità della nostra popolazione virtuale.

Nel nostro esempio, genereremo una popolazione con media teorica 100 e deviazione standard teorica 15 (la scelta di media e deviazione standard è arbitraria).

Generare dei campioni

A partire dal vettore popolazione, è possibile estrarre un vettore campione (di numerosità $m < n$). Per fare questo, R mette a disposizione la funzione sample(x,m,replace=FALSE), dove x è la popolazione e m è la numerosità del campione.

Per visualizzare la distribuzione dell’errore di campionamento, utilizzeremo una popolazione di 10000 valori, e genereremo 200 campioni di numerosità 50.

Dunque n = 10000 (numerosità della popolazione simulata), k = 200 (numero di campioni), m = 50 (osservazioni per campione). Poi, genereremo anche una serie di campioni da 20 osservazioni.

Analisi descrittiva

Una volta generati questi dati, possiamo utilizzare alcune tecniche di analisi univariata per fare delle misurazioni.

In primo luogo possiamo calcolare la media e la deviazione standard della popolazione. Ci aspetteremo che la prima sia prossima a 100 e la seconda a 15. Poi, possiamo visualizzare un istogramma con la distribuzione della popolazione, che ci aspettiamo sia di tipo normale. Per verificarlo, possiamo usare le funzioni qqnorm e qqline.

n_pop <- 10000 # numerosità della popolazione
m50 <- 50 # osservazioni per campione
K <- 200  # numero di campioni da generare
media_teorica <- 100 # media teorica della popolazione
sd_teorica <- 15 # deviazione standard teorica della popolazione
popolazione <- rnorm(n_pop, media_teorica, sd_teorica)
mean(popolazione) # la media della popolazione

## [1] 100.0148

sd(popolazione)   # la deviazione standard della popolazione

## [1] 14.92306

hist(popolazione) # l'istogramma

Utilizzando qqnorm, valutiamo la normalità della distribuzione

qqnorm(popolazione)
qqline(popolazione, col = 2)

Ora, creiamo una matrice 200*50. Ogni riga rappresenta un campione di 50 osservazioni. Popoliamo le righe con la funzione sample, che campiona 50 osservazioni dalla popolazione.

Con medie_campioni50 <- apply(campioni50, 1, mean), calcoliamo la media di ogni campione e la salviamo nel vettore (di lunghezza 200) medie_campioni50. Su questo vettore calcoliamo la media e la deviazione standard (che rappresentano la media delle medie e la deviazione standard delle medie, ovvero l'errore standard.

campioni50 <- matrix(sample(popolazione, m50*K,replace = TRUE),nrow = K, ncol = m50)
## for (k in 1:K) {campioni50[k, ] <- sample(popolazione, m50)}
medie_campioni50 <- apply(campioni50, 1, mean)
mean(medie_campioni50)

## [1] 99.70899

sd(medie_campioni50)

## [1] 2.111001

L'istogramma della distribuzione campionaria

hist(medie_campioni50)

qqnorm(medie_campioni50)
qqline(medie_campioni50, col = 2)

Testiamo la normalità della distribuzione delle medie campionarie, usando lo Shapiro-Wilk normality test.

shapiro.test(medie_campioni50)

## 
## 	Shapiro-Wilk normality test
## 
## data:  medie_campioni50
## W = 0.99176, p-value = 0.3173

Leggere i risultati

La funzione shapiro.test restituisce un :riferimento:. Questo significa che non è rifiutata l'ipotesi di normalità. Dunque, non vi è violazione della normalità della distribuzione.

Campioni di numerosità 20

Ripetiamo la procedura, ma questa volta generiamo campioni di 20 osservazioni. Questo passaggio ci serve per capire se e come cambia la distribuzione campionaria al variare della numerosità del campione.

m20 <- 20
campioni20 <- matrix(nrow = K, ncol = m20)
campioni20 <- matrix(sample(popolazione, m20*K,replace = TRUE),nrow = K, ncol = m20)
medie_campioni20 <- apply(campioni20, 1, mean)
mean(medie_campioni20)

## [1] 99.70112

sd(medie_campioni20)

## [1] 3.60575

Intervallo di confidenza

A partire da queste simulazioni, possiamo introdurre il concetto di intervallo di confidenza.

Conoscendo la popolazione, possiamo prevedere il valore esatto della media di un campione di numerosità m estratto casualmente? La risposta, abbiamo visto, è negativa. Possiamo però stimare un intervallo entro il quale possiamo prevedere dove questa media verrà a cadere.

Nemmeno l'intervallo, però, può garantirci la sicurezza al 100%, in quanto non possiamo escludere di incorrere in campionamenti particolarmente sbilanciati da una parte o dall'altra.

La cosa più ragionevole da fare è quella di stabilire un livello di rischio percentuale accettabile, e di calcolare l'intervallo in base a questo rischio.

Calcolare il range

Detto in altri termini, possiamo calcolare i valori minimo e massimo, e dunque il range, entro il quale, probabilmente, il (100-rischio)% delle medie dei campioni andrà a cadere.

Se, ad esempio, consideriamo accettabile un rischio del 5%, calcoleremo il range entro il quale si collocano le medie del 95% dei campioni estratti. Questo ci permette di tagliare le code estreme, a destra e a sinistra, della distribuzione.

Per fare questo, tagliamo il 2.5% di campioni con media più bassa e il 2.5% di campioni con media più alta.

La media del campione con media più bassa rimanente, e la media del campione con media più alta rimanente, costituiscono il range che cercavamo, ovvero l'intervallo di confidenza. Per calcolare questi valori, sarà sufficiente estrarre i percentili 2.5 e 97.5 dalla distribuzione delle medie dei campioni.

confidenza_campioni50 <- quantile(medie_campioni50, probs = c(0.025, 0.975))
confidenza_campioni20 <- quantile(medie_campioni20, probs = c(0.025, 0.975))
confidenza_campioni50

##      2.5%     97.5% 
##  95.36919 103.36483

confidenza_campioni20

##      2.5%     97.5% 
##  92.75594 106.02945

differenza_campioni50 <- confidenza_campioni50[2] - confidenza_campioni50[1]
differenza_campioni20 <- confidenza_campioni20[2] - confidenza_campioni20[1]

Come possiamo notare, l'intervallo di confidenza della distribuzione campionaria è diverso, cambiando la numerosità dei campioni. Nel caso di campioni di numerosità 20, il range è di 13.2735122, mentre per i campioni di numerosità 50, è di 7.9956332.

Confrontare le due distribuzioni

Ora, usiamo la funzione density per confrontare le due distribuzioni, quella dei campioni di 50 osservazioni, e quella dei campioni di 20. Abbiamo disegnato anche le due medie (le righe verticali) e gli intervalli di confidenza (le righe orizzontali).

density20 <- density(medie_campioni20)
density50 <- density(medie_campioni50)
plot(density20, ylim = c(0, max(density50$y)), col = 2, lty = 2)
# linea verticale con la media dei campioni da 20
abline(v = mean(medie_campioni20), col = 2, lty = 2) 
lines(x = confidenza_campioni20, y = c(0.02, 0.02), col = 2, lty = 2)
lines(density50, col = 3, lty = 4)
abline(v = mean(medie_campioni50), col = 3, lty = 4)
lines(x = confidenza_campioni50, y = c(0.03, 0.03), col = 3, lty = 4)

Possiamo notare come la distribuzione dei campioni di 20 osservazioni sia più larga di quella da 50. Corrispondentemente, anche i due intervalli di confidenza sono diversi.

Numerosità dei campioni e varianza

La varianza della distribuzione delle medie dei campioni costituisce una stima dell'errore di campionamento: più bassa la varianza (e la sd), più basso l'errore, e viceversa.

Dalle nostre simulazioni, si può intuire che l'entità dell'errore è legato alla numerosità del campione. <!-- %Se m == n, ovvero se il campione coincidesse con la popolazione, l'errore sarebbe pari a zero. %L'errore è massimo nella situazione opposta, ovvero se decidiamo di utilizzare campioni di numerosità 1.

%In questo caso, la media delle medie dei campioni (che, nella fattispecie, corrisponde al valore dell'unica osservazione) continuerebbe ad essere pari a 20 e la deviazione standard della distribuzione dei campioni si approssimerebbe a 2, la sd della popolazione. -->

Generare molti campioni da un campione

La simulazione presentata nei paragrafi precedenti, seppur utile da un punto di vista didattico, è irrealistica: il ricercatore non può lavorare sull'intera popolazione, ma soltanto un campione.

Inoltre, il compito del ricercatore è quello di stimare la media della popolazione partendo dal campione, e non il contrario.

Partiamo da due osservazioni

• Se abbiamo a disposizione soltanto il nostro campione, la miglior stima della media della popolazione è la media del campione stesso.
• La distribuzione del campione dovrebbe "assomigliare" (al netto dell'errore statistico) alla distribuzione della popolazione. E, dunque, la distribuzione del campione è la miglior stima che abbiamo della distribuzione della popolazione.

Se assumiamo che la media della popolazione sia pari a quella del campione, e che anche la distribuzione sia paragonabile, possiamo immaginare di generare numerosi campioni fittizzi a partire dal campione noto.

Bootstrapping

Questa tecnica è nota come bootstrapping, e permette di calcolare l'intervallo di confidenza di un parametro (quale, ad esempio, la media).

Per generare un nuovo campione dal campione esistente, basta estrarre a caso m osservazioni dal campione originale.

Naturalmente, l'estrazione dev'essere con ripetizione. In caso contrario, il nuovo campione sarebbe identico. Dunque, alcuni elementi verranno estratti più di una volta, altri nessuna.

Percentili e intervallo di confidenza

In questo modo, possiamo generare dei nuovi campioni dal campione esistente. Anche in questo caso possiamo calcolare la media per ogni nuovo campione, e calcolare la distribuzione delle medie.

A partire da questa distribuzione, possiamo calcolare l'intervallo di confidenza, partendo dai percentili. Useremo i percentili 2.5 e 97.5 per un intervallo di confidenza del 95% (e un errore del 5%).

Per iniziare, prendiamo ora il primo dei 200 campioni generati, ed usiamolo per il bootstrapping. Calcoliamo il vettore delle medie dei bootstrap. Calcoliamo la media delle medie.

campioneA <- campioni50[1, ]
mean(campioneA)

## [1] 94.67507

bootstraps <- matrix(sample(campioneA, size = 10000, replace = TRUE), nrow = k, ncol = m50, byrow = TRUE)
medie_bootstraps <- apply(bootstraps, 1, mean)
confidenza_bootstraps <- quantile(medie_bootstraps, probs = c(0.025, 0.975))
mean(medie_bootstraps)

## [1] 94.5916

confidenza_bootstraps

##     2.5%    97.5% 
## 90.65777 98.90072

Confronto fra le distribuzioni

Nella sezione precedente, abbiamo visto la situazione ideale ma improbabile: conoscere l'intera popolazione, estrarre k campioni, calcolare la media per ognuno dei campioni; in questo modo abbiamo la vera distribuzione campionaria, di cui possiamo calcolare media e varianza.

In questa sezione, vediamo una situazione più realistica: abbiamo un campione, e lavoriamo su quello. Attraverso il bootstrapping, generiamo k campioni virtuali, e calcoliamo la distribuzione campionaria virtuale. Per capire se il secondo algoritmo, realistico ma virtuale, produce risultati robusti, confrontiamo le due medie, i due intervalli di confidenza e le due distribuzioni nel grafico \ref.

densityboot <- density(medie_bootstraps)
plot(density50, ylim = c(0, max(c(density50$y, densityboot$y))), col = 2, lty = 2)
abline(v = mean(medie_campioni50), col = 2, lty = 2)
lines(x = confidenza_campioni50, y = c(0.2, 0.2), col = 2, lty = 2)
lines(densityboot, col = 3, lty = 4)
abline(v = mean(medie_bootstraps), col = 3, lty = 4)
lines(x = confidenza_bootstraps, y = c(0.3, 0.3), col = 3, lty = 4)
abline(v = mean(popolazione))

Le due distribuzioni, seppure non identiche, sono molto simili.
Anche gli intervalli di confidenza sono paragonabili: gli intervalli calcolati sui 200 campioni sono pari a 19.44196 e 20.46475; gli intervalli calcolati attraverso il bootstrapping sono 19.56849 e 20.50061.

Possiamo dunque intuitivamente affermare che il metodo del bootstrapping riesce a simulare, in maniera piuttosto precisa, la distribuzione campionaria.

Usare l'approccio parametrico

Il calcolo parametrico dell'intervallo di confidenza è l'argomento del prossimo capitolo. Qui, ci limitiamo ad anticipare i risultati del test parametrico t.test.

# l'intervallo di confidenza osservato estraendo 200 campioni
confidenza_campioni50

##      2.5%     97.5% 
##  95.36919 103.36483

# il calcolo parametrico dell'intervallo
t_campioneA <- t.test(campioneA)
t_campioneA$conf.int

## [1] 90.48922 98.86093
## attr(,"conf.level")
## [1] 0.95

# l'intervallo di confidenza calcolato con i bootstrap
confidenza_bootstraps

##     2.5%    97.5% 
## 90.65777 98.90072

Come vediamo, il metodo parametrico e il metodo che usa il bootstrap non restituiscono risultati uguali, ma molto simili.