T test: confronto fra medie di due campioni
Introduzione
Confronto fra due medie
Nel capitolo ch:chiquadro, abbiamo introdotto il confronto fra variabili di tipo categoriale. In questo capitolo, affronteremo la statistica che permette di valutare la relazione fra una variabile di tipo categoriale ed una numerica (ad intervalli o rapporti). Nel caso specifico, ci limitiamo alla circostanza in cui la variabile indipendente, categoriale, ha due sole categorie.
In questa circostanza, le osservazioni sulla variabile dipendente, numerica, vengono divise in due insiemi. Il questito inferenziale che ci si pone è di valutare se i valori della variabile dipendente, di tipo numerico, differiscono significativamente da un gruppo all'altro.
Nel confronto fra due campioni si può adottare l'approccio parametrico, utilizzando il t-test, oppure un approccio non parametrico. In questo capitolo, verrà introdotto prima l'approccio non parametrico, in quanto più intuitivo, e dopo l'approccio parametrico. In entrambi i casi, l'argomento verrà affrontato prima con un approccio simulativo, e poi utilizzando la distribuzione teorica di riferimento. Infine, verrà utilizzata la corrispondente funzione di R e ne verranno letti i risultati.
Calcolo non parametrico
Approccio intuitivo
Intuitivamente, possiamo dire che i due gruppi differiscono se le misurazioni di un gruppo sono, in genere, sistematicamente più alte (o più basse) delle misurazioni sull'altro gruppo.
Se tutte le osservazioni di un gruppo (chiamiamolo gruppo B) sono più elevate di tutte le osservazioni dell'altro (gruppo A), possiamo inferire che, relativamente alla variabile misurata, vi è una differenza significativa fra il gruppo B ed il gruppo A.
Seguendo questa intuizione, un metodo per valutare se la variabile indipendente, categoriale, ha una relazione sulla variabile dipendente (numerica), è quello di confrontare ogni elemento del gruppo A con ogni elemento del gruppo B.
Confronto fra elementi
Sempre facendo ricorso all'intuizione, appare chiaro che se il numero di confronti vinti da un gruppo è molto superiore al numero di confronti vinti dall'altro, la differenza è significativa.
In questo conteggio, possiamo arrivare a due condizioni estreme: nella prima, gli elementi di A vincono tutti i confronti nei confronti di tutti gli elementi di B; nella seconda, sono gli elementi di B a vincere tutti i confronti.
L'ipotesi nulla, $H_0$, assume la parità fra il numero di confronti vinti da A e vinti da B.
La simulazione
Errore di campionamento
Sappiamo però che è molto improbabile che il confronto fra i due gruppi sia perfettamente pari. Come sappiamo, infatti, anche nella circostanza in cui i due campioni sono estratti dalla stessa popolazione ed assegnati ad una o all'altra categoria a caso, emergeranno delle differenze dovute all'errore di campionamento.
Generare la distribuzione dell'errore
Per stimare l'entità di questo errore, e misurare la probabilità che una differenza sia attribuibile o meno ad esso, possiamo usare la stessa metodologia vista nei capitoli precedenti:
- generare k coppie di campioni, estratte ed assegnate casualmente;
- calcolare il numero di confronti vinti dall'uno e dall'altro gruppo;
- salvare questi valori in un vettore, che costituisce la distribuzione dell'errore di campionamento.
Confrontare una coppia di campioni
Potendo disporre della distribuzione dell'errore di campionamento, data una coppia di campioni possiamo calcolare il numero di vittorie dell'uno e dell'altro gruppo, e valutare dove si collocano nella distribuzione.
Se il risultato di questi confronti si colloca sulle code della distribuzione, possiamo rifiutare l'ipotesi nulla ed accettare l'ipotesi alternativa, ovvero che vi è una differenza significativa fra i due gruppi.
In caso contrario, non si rifiuta l'ipotesi nulla.
R: genero la popolazione, due campioni
Generiamo una popolazione di 10.000 unità, con media 20, sd 2 e distribuzione normale, usando rnorm.
n <- 10000 # numerosità della popolazione
m <- 100 # numerosità di ognuno dei 2 campioni
k <- 10000 # numero di coppie di campioni generati
media_teorica <- 20 # arbitrario
sd_teorica <- 2 # arbitrario
popolazione <- rnorm (n, media_teorica, sd_teorica)
R: calcoliamo i confronti fra i due campioni
Estraiamo ora 2 campioni dalla popolazione. Instanziamo due contatori, vinceA e vinceB. Con un ciclo for annidato, confrontiamo ogni valore del campione1 con ogni valore del campione2. Incrementiamo il contatore vinceA quando a vincere è l'unità del campione1, incrementiamo vinceB quando vince campione2 (essendo valori a virgola mobile non è importante gestire il pareggio).
campione1 <- sample(popolazione,m,replace=FALSE)
campione2 <- sample(popolazione,m,replace=FALSE)
vinceA <- 0; vinceB <- 0
for (x in 1:m) {
oss1 <- campione1[x]
for (y in 1:m) {
oss2 <- campione2[y]
if (oss1>oss2) {
vinceA <- vinceA + 1 # incremento vinceA
} else {
vinceB <- vinceB + 1 # incremento vinceB
}
}
}
R: risultati
Non dovrebbe sorprenderci il fatto che la somma dei due valori è pari a m*m, ovvero il numero dei confronti.
c(vinceA, vinceB)
## [1] 4484 5516
Questa statistica viene chiamata Mann-Whitney-Wilcoxon U
Calcolo della posizione ordinale
È però possibile ottenere lo stesso risultato con un calcolo diverso: l'insieme di tutte le osservazioni viene ordinato, e ad ogni osservazione viene assegnato un punteggio pari alla sua posizione;
- per ognuno dei due gruppi, si somma il punteggio di ogni osservazione;
- a questi valori, si sottrae quello che è il minimo valore possibile, ovvero m*(m+1)/2.
Il vantaggio di questo algoritmo è che può essere esteso a confronti fra più di due gruppi.
R: il calcolo del ranking
due_rank <- matrix(rank(c(campione1,campione2)),
nrow=2,ncol=m, byrow=TRUE)
due_somma_rank <- apply (due_rank,1,sum)
rank_atteso <- m*(m+1)/2
wilcoxon <- due_somma_rank-rank_atteso
wilcoxon
## [1] 4485 5515
La simulazione
Introdotta la statistica, usiamo la simulazione. Generiamo 2 vettori. Nel primo, distribuzione, inseriremo le coppie di valori che calcoleremo usando la statistica di Wilcoxon. Nel secondo, differenze, salviamo una seconda statistica: la differenza delle medie fra i due campioni. La prima distribuzione ci serve in questa sezione, la seconda nella sezione dedicata al calcolo parametrico.
distribuzione <- vector(mode = "numeric", length = 20000)
differenze <- vector(mode = "numeric", length = length(distribuzione)/2)
A questo punto, usando un ciclo for, possiamo generare k=10000 coppie di campioni, calcolare per ognuno le due statistiche, e salvarla nei due vettori.
R: la generazione delle coppie di campioni
for (i in 1:length(distribuzione)/2) {
due_campioni <- matrix(sample(popolazione,2*m,replace=FALSE),nrow=2,ncol=m, byrow=TRUE)
campione1 <- due_campioni[1,]
campione2 <- due_campioni[2,]
due_rank <- matrix(rank(c(campione1,campione2)),nrow=2,ncol=m, byrow=TRUE)
due_somma_rank <- apply (due_rank,1,sum)
rank_atteso <- m*(m+1)/2
wilcoxon <- due_somma_rank-rank_atteso
distribuzione[i*2-1] <- wilcoxon[1]
distribuzione[i*2] <- wilcoxon[2]
differenze[i]<-mean(campione1)-mean(campione2)
}
R: la distribuzione U
par(mfrow=c(1,2))
hist(distribuzione)
qqnorm(distribuzione)
qqline(distribuzione)
R: valori critici
Possiamo calcolare i valori critici, ad esempio per $\alpha$ = 0.01 e 0.05 bidirezionale
quantile(distribuzione, probs = c(0.005,0.025,0.975,0.995))
## 0.5% 2.5% 97.5% 99.5%
## 3960.00 4206.95 5822.00 6067.00
Calcolare il p-value dalla distribuzione
Ora, generiamo due nuovi campioni, calcoliamo la statistica U, e vediamo dove si colloca rispetto alla distribuzione.
due_campioni <- matrix(sample(popolazione,2*m,replace=FALSE),nrow=2,ncol=m, byrow=TRUE)
campione1 <- due_campioni[1,]
campione2 <- due_campioni[2,]
due_rank <- matrix(rank(c(campione1,campione2)),nrow=2,ncol=m, byrow=TRUE)
due_somma_rank <- apply (due_rank,1,sum)
rank_atteso <- m*(m+1)/2
Calcolare il p-value dalla distribuzione
(wilcoxon <- due_somma_rank-rank_atteso)
## [1] 4522 5478
(p_value_simulazione <- rank(c(wilcoxon,distribuzione))[1:2]/length(distribuzione))
## [1] 0.116125 0.878750
La distribuzione U Mann-Whitney-Wilcoxon
La distribuzione di errore che abbiamo ottenuto grazie al confronto di 10000 coppie di campioni, è nota come distribuzione U Mann-Whitney-Wilcoxon [@WhitleyBall:2002].
R: le funzioni per la distribuzione U
R mette a disposizione il gruppo di funzioni per calcolare la densità, la probabilità, per generare dei numeri secondo la distribuzione. Inoltre, mette a disposizione la funzione wilcox.test per calcolare, automaticamente, la statistica ed il p-value.
- dwilcox(x, m, n) calcola la densità di x
- pwilcox(q, m, n) calcola la probabilità
- rwilcox(nn, m, n) genera nn numeri casuali.
- wilcox.test(gruppoA,gruppoB) calcola il test corrispondente
m e n sono la numerosità del primo e del secondo campione (che, nel nostro esempio, sono uguali -- m)
R: calcolo del p-value con pwilcox
Utilizzando pwilcox calcoliamo i due p-value. Il valore interessante è quello più basso. Se, come nell'esempio, l'ipotesi è a due vie, dobbiamo raddoppiare il p-value
(p_value_Wilcoxon <-pwilcox(wilcoxon, 100,100))
## [1] 0.1219539 0.8785401
# se il test è a due code, il p-value va moltiplicato per 2
p_value_simulazione*2
## [1] 0.23225 1.75750
p_value_Wilcoxon*2
## [1] 0.2439078 1.7570803
wilcox.test(campione1,campione2,exact = TRUE)
##
## Wilcoxon rank sum exact test
##
## data: campione1 and campione2
## W = 4522, p-value = 0.2439
## alternative hypothesis: true location shift is not equal to 0
Come sempre, i risultati della distribuzione generata non sono uguali, ma simili, a quelli della distribuzione teorica.
R: uso di wilcox.test
Vediamo ora il calcolo effettuando la funzione di R wilcox.test
wilcox.test(campione1,campione2,exact = TRUE)
##
## Wilcoxon rank sum exact test
##
## data: campione1 and campione2
## W = 4522, p-value = 0.2439
## alternative hypothesis: true location shift is not equal to 0
La funzione calcola la statistica W e il p-value.
Approccio parametrico
La differenza delle medie
L'algoritmo utilizzato nelle sezioni precedenti costituisce l'approccio non parametrico al confronto fra due campioni.
Come abbiamo visto, il calcolo non parametrico non prende in considerazione né la media né la deviazione standard e nemmeno la distribuzione dei campioni.
Il vantaggio del calcolo non parametrico è che fa pochissime assunzioni:
- le m+n osservazioni devono essere indipendenti;
- m e n devono avere una numerosità di almeno 5 elementi.
Assunzioni
Per applicare il test parametrico, al contrario, è necessario non solo che le osservazioni siano indipendenti (e la numerosità adeguata). È altresì necesario che:
- la distribuzione dei due campioni sia normale;
- la varianza dei due gruppi non sia diversa.
R permette il calcolo del t test anche in caso di varianze differenti (attraverso l'approssimazione di Welch).
La differenza fra le medie
Il test parametrico si basa sulla differenza fra le medie dei due campioni. Questo spiega l'assunto di normalità dei campioni.
In pratica, il test calcola il valore assoluto della differenza fra le due medie, e la confronta con la distribuzione t di Student.
Naturalmente, anche in questo caso possiamo calcolare il p-value ignorando la distribuzione teorica, ma basandoci sulla distribuzione dell'errore di campionamento generata dal confronto delle nostre 10.000 coppie di campioni.
Nel ciclo for che usammo per generare la statistica U, popolammo anche un vettore, differenze, con il codice differenze[i]<-mean(campione1)-mean(campione2). Possiamo ora usare quel vettore di distribuzioni dell'errore.
R: p-value usando la distribuzione
Calcoliamo il p-value confrontando la distanza delle due medie con la distribuzione dell'errore.
distanza<-abs(mean(campione1)-mean(campione2))
(p_value_differenze_simulazione <- 1-rank(c(distanza,differenze))[1]/length(differenze))
## [1] 0.1237
# se test a due code, va raddoppiato
p_value_differenze_simulazione*2
## [1] 0.2474
Distribuzione dell'errore
Usiamo la funzione density per visualizzare la distribuzione dell'errore, e qqnorm-qqline per testarne la normalità.
par(mfrow=c(1,2))
plot(density(differenze))
qqnorm(differenze)
qqline(differenze,col="red")
Varianza dell'errore
Possiamo osservare dal grafico qqnorm che la distribuzione approssima la distribuzione normale.
In realtà, la varianza della distribuzione dell'errore è legata alla numerosità dei due campioni. Più precisamente, la varianza stimata è pari a $$S_{\bar-\bar}^2 = \frac{s_1^2 (m-1) + s_2^2(n-1)}{m+n-2} (\frac{1}+\frac{1})$$
Dunque, l'errore standard della differenza fra le medie può essere calcolato, con R, usando
errore_standard <- sqrt((var(campione1)*(m-1)+var(campione2)*(m-1))/
(m-1+m-1)*(1/m+1/m))
Calcolo di t
Calcoliamo t, e calcoliamo il p-value (che raddoppiamo, se l'ipotesi è bidirezionale).
t <- distanza/errore_standard
(p_value_differenze_t <- (1 - pt(t, df = (m-1+m-1))))
## [1] 0.1419497
# se test a due code, va raddoppiato
p_value_differenze_t*2
## [1] 0.2838993
Uso della funzione t.test
Infine, utilizziamo la funzione t.test
(t_test <- t.test(campione1,campione2))
##
## Welch Two Sample t-test
##
## data: campione1 and campione2
## t = -1.0745, df = 198, p-value = 0.2839
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.9321028 0.2745954
## sample estimates:
## mean of x mean of y
## 19.86102 20.18978
Leggere l'output
La funzione restituisce t = -1.0745156, i gradi di libertà stimati 197.9997053, il p-value = 0.2838993, e l'intervallo di confidenza della differenza fra le medie: -0.9321028, 0.2745954. Quando i due termini dell'intervallo di confidenza hanno segni opposti, la differenza non è significativa.