Capitolo 6 Decisioni ed incertezza

Nel capitolo precedente abbiamo visto come è possibile esprimere l’incertezza che il campionamento e, in genere, l’errore sperimentale producono sulle nostre stime. In particolare, abbiamo visto che, per una certa statistica rilevata su un campione, è possibile costruire una sampling distribution (o sample space o distribuzione campionaria), che descrive la variabilità della statistica stessa tra un campionamento e l’altro. La sampling distribution è, almeno approssimativamente, gaussiana e la sua deviazione standard (detta errore standard) può essere utilizzata per definire una ‘banda’ di incertezza per la nostra stima, con un procedimento che prende il nome di inferenza statistica (stima per intervallo). Analogamente, la sampling distribution può essere utilizzata per prendere decisioni in presenza di incertezza, con un procedimento che si chiama test d’ipotesi. Anche in questo capitolo, vediamo alcuni semplici, ma realistici esempi.

6.1 Confronto tra due medie: il test t di Student

Un ricercatore ha testato due genotipi (A e P) in un disegno sperimentale a randomizzazione completa, con cinque repliche (dieci parcelle in totale). Alla fine dell’esperimento ha determinato la produzione parcellare, ottenendo quindi dieci valori, che debbono essere considerati come un campione di tutti quelli possibili; per ogni genotipo, i valori sono diversi, a seguito degli effetti stocastici di varia natura, agronomica ed ambientale, che costituiscono il meccanismo perturbativo di fondo.

I risultati sono i seguenti:

A <- c(65, 68, 69, 71, 78)
P <- c(80, 81, 84, 88, 94)

Nel campione A la media è pari a 70.2, mentre la deviazione standard è pari a 4.87. L’errore standard è pari a 2.18 e quindi l’intervallo di confidenza della media è 70.2 \(\pm\) 6.04. Invece, nel campione P, la media è 85.4, mentre la deviazione standard è pari a 5.72. L’errore standard è pari a 2.56, mentre l’intervallo di confidenza per la media è 85.4 \(\pm\) 7.11.

Dopo aver completato questo esperimento, ci chiediamo se sia possibile concludere che il genotipo P è più produttivo del genotipo A, coerentemente con i risultati osservati. Nel rispondere a questa domanda bisogna tener presente che i due campioni ottenuti sono totalmente irrilevanti, dato che il nostro interesse è rivolto alle popolazioni che hanno generato i campioni; vogliamo infatti che le nostre conclusioni abbiano carattere di universalità e non siano specifiche per il nostro esperimento. Dobbiamo quindi trovare un metodo per decidere se il genotipo P, in generale, è più produttivo del genotipo A, pur in presenza dell’incertezza legata all’errore sperimentale.

Un primo approccio intuitivo potrebbe essere basato sugli intervalli di confidenza delle due medie. Possiamo notare che il limite di confidenza superiore per A (70.2 + 6.04 = 76.24) è inferiore al limite di confidenza inferiore per P (85.4 - 7.11 = 78.29), in modo che gli intervalli di confidenza non si sovrappongono (Figura 6.1). In base a questo criterio, quindi, potremmo concludere che le popolazioni da cui provengono i due campioni sono diverse e, di conseguenza, P è un genotipo più produttivo.

Medie ed intervalli di confidenza calcolati sue due campioni. Si può notare che gli intervalli di confidenza non si sovrappongono.

Figura 6.1: Medie ed intervalli di confidenza calcolati sue due campioni. Si può notare che gli intervalli di confidenza non si sovrappongono.

Anche se questo criterio è, in qualche modo, accettabile, in pratica si preferisce utilizzare un altro criterio più rigoroso, che illustreremo di seguito.

6.1.1 L’ipotesi nulla e alternativa

Innanzitutto, ricordiamo la logica Popperiana illustrata nel primo capitolo, secondo la quale nessun esperimento può dimostrare che un’ipotesi scientifica è vera, mentre è possibile dimostrare che essa è falsa. E’quindi conveniente porre la nostra ipotesi scientifica in modo che essa possa essere falsificata; dovendo dimostrare che la produzione di A è diversa da quella di P, possiamo formulare l’ipotesi scientifica (\(H_0\)) in questo modo:

\[H_0: \mu_A = \mu_P = \mu\]

In altre parole, la nostra ipotesi di lavoro è che i due campioni siano in realtà estratti da due distribuzioni normali con la stessa media e la stessa deviazione standard, il che equivale a dire che essi provengono da un’unica distribuzione normale con media \(\mu\) e deviazione standard \(\sigma\). Questa ipotesi si chiama ipotesi nulla e, se riuscissimo a falsificarla, avremmo conseguito il nostro obiettivo, in totale coerenza con la logica Popperiana. Vi chiediamo di fare nuovamente attenzione al fatto che l’ipotesi nulla riguarda le medie delle popolazioni che hanno generato i campioni, non i campioni stessi, per i quali già sappiamo che le medie sono diverse.

Oltre all’ipotesi nulla possiamo anche definire l’ipotesi alternativa, che abbracceremmo se dovessimo riuscire a falsificare quella nulla. L’ipotesi alternativa semplice è:

\[H_1 :\mu_A \ne \mu_P\]

Se avessimo elementi sufficienti già prima di effettuare l’esperimento (e non dopo aver visto i risultati), potremmo anche adottare ipotesi alternative complesse, del tipo

\[H_1 :\mu_A > \mu_P\]

oppure:

\[H_1 :\mu_A < \mu_P\]

Abbiamo già anticipato che le ipotesi (nulla ed alternativa) debbono essere stabilite prima di effettuare l’esperimento. Che cosa è possibile attendersi prima di un esperimento, in relazione alla produttività di due genotipi? In mancanza di informazioni preliminari attendibili, dovremmo essere indotti ad abbracciare l’ipotesi alternativa semplice \(\mu_A \neq \mu_P\), mentre una delle due ipotesi alternative complesse può essere adottata soli in quei casi in cui l’altra può essere esclusa ‘a priori’, già prima di aver visto i risultati dell’esperimento.

6.1.2 La statistica T

Ma torniamo a lavorare sull’ipotesi nulla. Possiamo intuire che l’idea che \(\mu_A = \mu_P\) non è impossibile; infatti, se avessimo una sola popolazione gaussiana con media \(\mu\), potremmo comunque campionare cinque parcelle sulla coda sinistra (e quindi con bassa produzione) e cinque parcelle sulla coda destra (con elevata produzione), in modo che \(m_A \neq m_P\). Il problema è capire quanto questo sia probabile.

Cerchiamo di definire una statistica che, posto \(\mu_A = \mu_P\), ci permetta di capire che possibilità ho di trovare \(m_A \neq m_P\). Due sono i valori di cui tenere conto:

  1. l’entità della differenza tra le medie dei campioni: più la differenza tra le due medie è alta e più è improbabile che essa si produca quando \(\mu_A = \mu_P\);
  2. l’entità dell’errore standard. Maggiore l’incertezza di stima delle due medie, maggiore la probabilità di trovare ampie differenze tra \(m_A\) ed \(m_P\) anche quando \(\mu_A = \mu_P\).

Una buona statistica è quella indicata di seguito:

\[T = \frac{m_A - m_P}{SED}\]

Si può osservare che T, in realtà, non è altro che il rapporto tra le quantità indicate in precedenza ai punti 1 e 2: infatti la quantità al numeratore è la differenza tra le medie dei due campioni, mentre la quantità al denominatore è il cosiddetto errore standard della differenza tra due medie (SED). Quest’ultima quantità si può ottenere pensando che i due campioni sono estratti in modo indipendente e, pertanto, la varianza della somma (algebrica) è uguale alla somma delle varianze. La varianza di ogni media è data dal quadrato dell’errore standard (SEM) e, pertanto, possiamo scrivere:

\[SED^2 = SEM_1^2 + SEM_2^2\]

Sappiamo anche che il SEM si ottiene dividendo la deviazione standard di ogni campione per la radice quadrata del numero dei dati, quindi:

\[SED^2 = \frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}\]

cioè:

\[SED = \sqrt{ \frac{s_1^2}{n_1} + \frac{s_2^2}{n_2} }\]

Possiamo anche scrivere:

\[SED = \sqrt{ \frac{s_1^2 \, n_2 + s_2^2 \, n_1}{n_1 \, n_2} }\]

e, se le varianze sono uguali (\(s_1^2 = s_2^2 = s^2\)), segue che:

\[SED = \sqrt {s^2 \frac{n_1 + n_2}{n_1 \, n_2 } }\]

Se fosse anche \(n_1 = n_2 =n\), potremmo scrivere:

\[SED = \sqrt{2 \, \frac{s^2}{n} } = \sqrt{2} \times SEM\]

Il valore da noi osservato per T è:

\[T = \frac{70.2 - 85.4}{3.361547} = -4.5217\]

dove il denominatore è ottenuto come:

\[SED = \sqrt{ 2.18^2 + 2.56^2 } = 3.361547\]

6.1.3 Simulazione Monte Carlo

Ipotizzando che l’ipotesi nulla sia vera (\(\mu_A = \mu_P\)), che valori può assumere la statistica T? E’facile intuire che il valore più probabile sia 0, ma, dato che la media di un campione non necessariamente coincide con la media della popolazione da cui proviene, possiamo anche aspettarci che T sia diverso da zero. Tuttavia, trovare valori di T molto alti o molto bassi dovrebbe essere via via meno probabile.

Nel nostro esperimento, il valore di T che abbiamo ottenuto è abbastanza diverso da zero, il che indica un certo grado di discrepanza tra l’osservazione è l’ipotesi nulla. Possiamo affermare che ciò sia imputabile solo alla variabilità di campionamento e che quindi il nostro esperimento confermi l’ipotesi nulla?

Per rispondere a questa domanda, supponiamo che l’ipotesi nulla sia vera. In questo caso, immaginiamo che le nostre dieci parcelle siano tutte estratte da una sola popolazione di parcelle con media e deviazione standard stimate (stima puntuale) come segue:

media <- mean(c(A, P))
devSt <- sd(c(A, P))
media
## [1] 77.8
devSt
## [1] 9.44928

Prendiamo quindi questa popolazione normale, con \(\mu = 77.8\) e \(\sigma = 9.45\), ed utilizziamo un generatore di numeri casuali gaussiani per estrarre numerose (100’000) coppie di campioni, calcolando, per ogni coppia, il valore T, come abbiamo fatto con la nostra coppia iniziale.

Il codice da utilizzare in R per le simulazioni è il seguente:

T_obs <- -4.521727
set.seed(34)
result <- rep(0, 100000)
for (i in 1:100000){
  sample1 <- rnorm(5, media, devSt)
  sample2 <- rnorm(5, media, devSt)
  SED <- sqrt( (sd(sample1)/sqrt(5))^2 +
                 (sd(sample2)/sqrt(5))^2 )
  result[i] <- (mean(sample1) - mean(sample2)) / SED
}

# Valutazione dei valori di T
mean(result)
## [1] -0.001230418
max(result)
## [1] 9.988315
min(result)
## [1] -9.993187
# Quanti valori sono più discrepanti del mio?
length(result[result < T_obs]) / 100000
## [1] 0.00095
length(result[result > - T_obs]) /100000
## [1] 0.00082

Come risultato (vettore ‘result’) otteniamo una lista di 100’000 valori di T, tutti compatibili con l’ipotesi nulla vera (\(\mu_A = \mu_P\)); in altre parole, otteniamo una sampling distribution per T, in quanto le variazioni tra un valore e l’altro sono solo dovute al campionamento, visto che l’ipotesi nulla è vera. In mezzo a questi 100’000 valori ne troviamo alcuni piuttosto alti (\(> 9\)) e piuttosto bassi (\(< - 9\)). Negli anni 20 del 1900, Fischer propose di utilizzare come ‘forza dell’evidenza scientifica’ proprio la probabilità di ottenere un risultato uguale o più estremo di quello osservato, supponendo vera l’ipotesi nulla. Per applicare questo criterio, dobbiamo partire dalla nostra osservazione (T = -4.521727) e considerare che il valore è negativo, ma solo perché abbiamo scritto la differenza come \(m_A - m_P\) invece che come \(m_P - m_A\). Quindi dobbiamo andarci a cercare nel vettore ‘result’ i valori che risultano minori di -4.5217 e maggiori di 4.5217, che sono più discrepanti di quello da noi osservato.

Vediamo che, dei 100’000 valori di T simulati assumendo vera l’ipotesi nulla, poco meno dell’uno per mille sono inferiori a -4.5217 e altrettanti sono superiori al suo reciproco (4.5217). In totale, la probabilità di osservare un valore di T uguale o più estremo di quello da noi osservato è molto bassa a pari allo 0.18% circa. Questo valore di probabilità è detto P-level (o P-value) e viene utilizzato come criterio decisionale: se esso è inferiore a 0.05 (5% di probabilità), come in questo caso, rifiutiamo l’ipotesi nulla ed abbracciamo l’alternativa, concludendo che i due trattamenti sono significativamente diversi tra loro (in termini di risposta prodotta nei soggetti trattati, ovviamente).

6.1.4 Soluzione formale

Eseguire una simulazione di Monte Carlo per costruire una sampling distribution per T non è sempre agevole e, pertanto, ci dobbiamo chiedere se esista una funzione di densità per rappresentare questa sampling distribution. Nel grafico sottostante abbiamo trasformato il vettore ‘result’ in una distribuzione di frequenze relative, considerando intervalli di ampiezza pari a 2.5, in modo da poter rappresentare la sampling distribution di T con un istogramma.

#Sampling distribution per T 
max(result);min(result)
## [1] 9.988315
## [1] -9.993187
b <- seq(-10, 10, by = 0.25)
hist(result, breaks = b, freq=F, 
  xlab = "T", ylab="Density", 
  xlim=c(-10,10), ylim=c(0,0.45), main="")
curve(dnorm(x), add=TRUE, col="blue")
curve(dt(x, 8), add=TRUE, col="red")
Sampling distribution empirica a confronto con una distribuzione normale (in rosso) e una distribuzione t di Student con otto gradi di libertà

Figura 6.2: Sampling distribution empirica a confronto con una distribuzione normale (in rosso) e una distribuzione t di Student con otto gradi di libertà

Vediamo che l’istogramma non è rappresentabile fedelmente con una curva gaussiana (curva blu in Figura 6.2), in quanto le code sono leggermente più alte. Oltre alla funzione di densità gaussian esiste un’altra funzione simile, detta ‘t di Student’, che, rispetto alla gaussiana, è caratterizzata da una più elevata densità nelle code e quindi è in grado di fornire una descrizione più precisa del nostro istogramma (curva rossa in Figura 6.2). La forma della funzione t di Student dipende dal numero di gradi di libertà; in questo caso ne abbiamo otto, in quanto abbiamo due campioni con cinque individui e, quindi, quattro gradi di libertà ciascuno.

La curva di densità di ‘t di Student’ ci permette di calcolare il P-level molto velocemente, senza dover ricorrere ad una simulazione Monte Carlo. Chiaramente, dato che l’ipotesi alternativa è quella semplice dobbiamo prendere entrambe le code (test a due code), utilizzando il codice sottostante.

2 * pt(-4.5217, 8, lower.tail=T) 
## [1] 0.00194554

Abbiamo moltiplicato per 2 il risultato, in quanto la funzione dt() fornisce la probabilità di trovare individui inferiori a -4.5217 (‘lower.tail = T’), ma, essendo la funzione simmetrica, la probabilità di trovare soggetti superiori a 4.5217 è esattamente la stessa.

Vediamo che il P-level ottenuto formalmente è simile a quello ottenuto empiricamente, ma, rispetto a quest’ultimo, è più preciso, in quanto con la simulazione di Monte Carlo non abbiamo potuto considerare, come avremmo dovuto, un numero infinito di repliche dell’esperimento.

Allo stesso valore, più semplicemente, si perviene utilizzando la funzione t.test():

t.test(A, P, var.equal = T)
## 
##  Two Sample t-test
## 
## data:  A and P
## t = -4.5217, df = 8, p-value = 0.001945
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -22.951742  -7.448258
## sample estimates:
## mean of x mean of y 
##      70.2      85.4

Gli argomenti di questa funzione sono i due vettori, oltre all’argomento ‘var.equal’, che in questo caso è stato impostato su ‘TRUE’, considerando i due campioni estratti da distribuzioni gaussiane con la stessa varianza. Torneremo su questo aspetto tra poco.

6.1.5 Interpretazione del P-level

Abbiamo detto che quando il P-level è inferiore a 0.05, concludiamo che vi sono prove scientifiche sufficientemente forti per rifiutare la nostra ipotesi di partenza.

Bisogna sottolineare come il P-level nella statistica tradizionale sia stato inizialmente proposto da Fisher come criterio di comportamento e non come un vero e proprio criterio inferenziale-probabilistico. Successivamente, Jarzy Neyman ed Egon Pearson, intorno al 1930, proposero di utilizzare il P-level come probabilità di errore di I specie, cioè come probabilità di rifiutare erroneamente l’ipotesi nulla. Tuttavia, trattandosi di una probabilità calcolata a partire da una sampling distribution, cioè da un’ipotetica infinita ripetizione dell’esperimento, essa non ha alcun valore in relazione al singolo esperimento effettivamente eseguito, come i due autori menzionati in precedenza hanno esplicitamente chiarito.

Di conseguenza, nel caso in esempio, affermare che abbiamo una probabilità di errore pari a 0.0019 nel rifiutare l’ipotesi nulla, rappresenterebbe un abuso: le nostre conclusioni potrebbero essere false o vere, ma non abbiamo alcun elemento per scegliere tra le due opzioni. Possiamo solo affermare che, se ripetessimo infinite volte l’esperimento e se l’ipotesi nulla fosse vera, otterremmo un risultato estremo come il nostro o più estremo solo in 2 casi (circa) su 1000. In altre parole, basando sempre le nostre conclusioni sul criterio anzidetto (rifiuto l’ipotesi nulla se il P-value è inferiore a 0.05) commettiamo un errore in non più del 5% dei casi, cosicché il P-level può essere guardato come la probabilità di ‘falso-positivo’ nel lungo periodo, non in relazione ad ogni singolo sforzo sperimentale di campionamento.

6.1.6 Tipologie alternative di test t di Student

Chiunque abbia utilizzato un computer per eseguire un test di t, si sarà accorto che è necessario scegliere tra tre procedure alternative. Infatti, abbiamo:

  1. t-test appaiato. In questo caso le misure sono prese a coppia sullo stesso soggetto e non sono quindi indipendenti.
  2. t-test omoscedastico. Le misure sono prese su soggetti diversi (indipendenti) e possiamo supporre che i due campioni provengano da due popolazioni con la stessa varianza.
  3. t-test eteroscedastico. Le misure sono prese su soggetti diversi, ma le varianze non sono omogenee.

Quello che abbiamo appena utilizzato con la funzione t.test() è un test omoscedastico, che suppone che le varianze dei campioni siano simili tra loro, almeno dello stesso ordine di grandezza (omogeneità delle varianze od omoscedasticità).

Se consideriamo invece la coppia di campioni riportata di seguito, possiamo notare che le varianza sono molto diverse e l’ipotesi di omoscedasticità non è difendibile.

D1
## [1] 12.06608 11.79470 11.85008 12.14712 12.14591
D2
## [1] 35.14014 35.20918 31.96391 34.51307 33.91213
mean(D1)
## [1] 12.00078
mean(D2)
## [1] 34.14769
var(D1)
## [1] 0.02798071
var(D2)
## [1] 1.767402

Pertanto, dobbiamo utilizzare un t-test eteroscedastico (test di Welch), nel quale il numero di gradi di libertà non è ottenuto per semplice somma dei gradi di libertà di ogni campione, ma è approssimato con la formula di Satterthwaite:

\[DF_s \simeq \frac{ \left( s^2_1 + s^2_2 \right)^2 }{ \frac{(s^2_1)^2}{DF_1} + \frac{(s^2_2)^2}{DF_2} }\]

Vediamo che se le varianze e i gradi di libertà sono uguali, la formula precedente riduce a:

\[DF_s = 2 \times DF\]

Nel nostro caso:

dfS <- (var(D1) + var(D2))^2 / 
  ((var(D1)^2)/4 + (var(D2)^2)/4)
dfS
## [1] 4.126621

Con R, il test di Welch si esegue utilizzando l’argomento ‘var.equal = F’:

t.test(D1, D2, var.equal = F)
## 
##  Welch Two Sample t-test
## 
## data:  D1 and D2
## t = -36.959, df = 4.1266, p-value = 2.326e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -23.79070 -20.50312
## sample estimates:
## mean of x mean of y 
##  12.00078  34.14769

Il t-test appaiato è invece quello in cui le misure sono state prese negli stessi soggetti, prima e dopo un certo trattamento sperimentale. Ad esempio, se le misure D1 e D2 fossero state rilevate sugli stessi soggetti (due misure per soggetto, in posizioni corrispondenti dei vettori D1 e D2), allora avremmo cinque soggetti invece che dieci e, di conseguenza, avremmo solo solo 4 gradi di libertà invece che 8:

t.test(D1, D2, var.equal=F, paired=T)
## 
##  Paired t-test
## 
## data:  D1 and D2
## t = -38.002, df = 4, p-value = 2.864e-06
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -23.76497 -20.52885
## sample estimates:
## mean difference 
##       -22.14691

6.2 Confronto tra due proporzioni: il test \(\chi^2\)

Il test di t è molto utile, ma soltanto nel caso in cui si abbia a che fare con caratteri quantitativi, cioè con variabili misurate su una scala continua, per le quali sia possibile calcolare statistiche descrittive, come appunto la media. Talvolta, i ricercatori sono interessati a rilevare caratteristiche qualitative, come ad esempio lo stato di una pianta in seguito ad un trattamento (morta o viva), il colore dei semi (si ricordino i piselli verdi e gialli di Mendel) ed altre caratteristiche che non sono misurabili su una scala continua.

Avendo a che fare con variabili qualitative, l’unica statistica rilevabile è il numero di soggetti che presentano le diverse modalità. Ad esempio, immaginiamo un esperimento per verificare se un coadiuvante aumenta l’efficacia di un insetticida. In questo esperimento, utilizziamo l’insetticida da solo e miscelato con il coadiuvante su due gruppi di insetti diversi. Nel primo gruppo (trattato con insetticida) contiamo 56 morti su 75 insetti trattate, mentre nel secondo gruppo (trattato con insetticida e coadiuvante) otteniamo 48 morti su 50 insetti trattati.

Nel capitolo 3 abbiamo visto che i risultati di questo esperimento si riducono ad una tabella di contingenze, per la quale sappiamo già calcolare una statistica che misuri la connessione tra variabili (trattamento e mortalità), detta \(\chi^2\):

counts <- c(56, 19, 48, 2)
tab <- matrix(counts, 2, 2, byrow = T)
row.names(tab) <- c("I", "IC")
colnames(tab) <- c("M", "V")
tab
##     M  V
## I  56 19
## IC 48  2
summary( as.table(tab) )
## Number of cases in table: 125 
## Number of factors: 2 
## Test for independence of all factors:
##  Chisq = 9.768, df = 1, p-value = 0.001776

La connessione è la statistica giusta per rispondere alla nostra domanda di ricerca; infatti ci stiamo chiedendo se la proporzione dei morti è indipendente dal tipo di trattamento oppure no. Il valore di \(\chi^2\) osservato è pari a 9.768, il che indica un certo grado di connessione. Infatti, ricordiamo che, in caso di indipendenza tra le variabili, \(\chi^2\) dovrebbe essere zero (Capitolo 3).

Tuttavia, noi non siamo interessati ai due campioni, in quanto i 125 soggetti osservati sono tratti da due popolazioni più ampie. Considerando queste due popolazioni, poniamo l’ipotesi nulla in questi termini:

\[H_o :\pi_1 = \pi_2 = \pi\]

Vediamo che, come negli altri esempio, l’ipotesi nulla riguarda i parametri delle popolazioni (\(\pi_1\) e \(\pi_2\)), non quelli dei campioni (\(p_1\) e \(p_2\)). Ci chiediamo: se l’ipotesi nulla fosse vera (\(\pi_1 = \pi_2\)), quale sarebbe la sampling distribution per \(\chi^2\)? E soprattutto, quanto sarebbe probabile ottenere un valore alto come il nostro o più alto?

Vediamo che l’output della funzione summary() ci da già il P-level, che, essendo inferiore a 0.05, ci consente di rigettare l’ipotesi nulla, affermando che esiste una differenza significativa tra l’effetto dell’insetticida quando è utilizzato da solo e quando è utilizzando in abbinamento con un coadiuvante.

6.2.1 Simulazione Monte Carlo

Per chi fosse interessato, proviamo a vedere come si costruisce la sampling distribution per \(\chi^2\) con una simulazione Monte Carlo. Possiamo utilizzare la funzione r2dtable(), che, partendo da una situazione in cui l’ipotesi nulla è vera (cioè \(\pi_1 = \pi_2\)) e quindi i due caratteri sono indipendenti, produce tante tabelle di contingenza (nel nostro caso 10’000), rispettando i totali marginali della nostra tabella di partenza. Le tabelle prodotte sono restituite come lista, quindi possiamo utilizzare la funzione lapply() per applicare ad ogni elemento della lista la funzione che restituisce il \(\chi^2\) (‘chiSim’).

chiSim <- function(x) summary(as.table(x))$stat
set.seed(1234)
tabs <- r2dtable(10000, apply(tab, 1, sum), apply(tab, 2, sum))
chiVals <- as.numeric( lapply( tabs, chiSim) )
length(chiVals[chiVals > 9.768])
## [1] 19

Vediamo che vi sono 19 valori più alti di quello da noi osservato, quindi il P-value è 0.0019.

6.2.2 Soluzione formale

In modo formale, si può dimostrare che, se \(n\) è sufficientemente grande (n > 30), il valore osservato di \(\chi^2\) segue appunto la distribuzione di probabilità \(\chi^2\), con un numero di gradi di libertà \(\nu\) pari al numero dei dati indipendenti, che, in questo caso, è pari ad 1. Infatti, una volta fissata una frequenza, le altre sono automaticamente definite, dovendo restituire i totali marginali. In R, possiamo utilizzare la funzione pchisq() per calcolare la probabilità di ottenere valori pari o superiori a 9.768:

pchisq(9.76801, 1, lower.tail=F)
## [1] 0.001775746

Allo stesso risultato, ma in modo più semplice, è possibile pervenire utilizzando la funzione chisq.test(), applicata alla tabella di contingenza:

chisq.test(tab, correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  tab
## X-squared = 9.768, df = 1, p-value = 0.001776

Come nel caso del ‘t’ di Student, abbiamo diversi tipi di test di chi quadro. In particolare, possiamo applicare o no la correzione per la continuità, che è necessaria quando il numero dei soggetti è piccolo (minore di 30, grosso modo). Nel nostro caso, non lo abbiamo ritenuto necessario ed abbiamo quindi aggiunto l’argomento ‘correct = F’.

6.4 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.