domenica 16 novembre 2008

Test di Wilcoxon dei ranghi con segno

Confronto delle medie tra 2 gruppi di campioni appaiati, metodo non parametrico.

Il sindaco di una città vuole verificare se i livelli di inquinamento si riducono chiudendo le strade del centro al traffico. A tal fine viene misurato il tasso di inquinamento ogni 60 minuti (dalle 8 alle 22: totale di 15 rilevazioni) in una giornata in cui il traffico è aperto, e in una giornata di chiusura al traffico. Ecco di seguito i valori di inquinamento atmosferico:

Con traffico: 214, 159, 169, 202, 103, 119, 200, 109, 132, 142, 194, 104, 219, 119, 234
Senza traffico: 159, 135, 141, 101, 102, 168, 62, 167, 174, 159, 66, 118, 181, 171, 112


Siamo di fronte a dati appaiati, perché esiste un vincolo tra le rilevazioni, consistente nel fatto che stiamo considerando la stessa città (con le sue peculiarità atmosferiche, ventilazione, etc.) seppure in due differenti giorni. Non potendo supporre una distribuzione gaussiana per i valori rilevati, dobbiamo procedere con un test non parametrico, il test di Wilcoxon dei ranghi con segno.

> a = c(214, 159, 169, 202, 103, 119, 200, 109, 132, 142, 194, 104, 219, 119, 234)
> b = c(159, 135, 141, 101, 102, 168, 62, 167, 174, 159, 66, 118, 181, 171, 112)
>
> wilcox.test(a,b, paired=TRUE)

Wilcoxon signed rank test

data: a and b
V = 80, p-value = 0.2769
alternative hypothesis: true location shift is not equal to 0


Essendo il p-value maggiore di 0.05, significa che le medie delle rilevazioni sono sostanzialmente rimaste inalterate (accettiamo l’ipotesi H0), quindi bloccare il traffico per un solo giorno non ha portato a nessun miglioramento per quanto riguarda l’inquinamento della città.
Il valore V = 80 corrisponde alla somma dei ranghi assegnati alle differenze con segno positivo. Con qualche acrobatismo matematico, possiamo calcolare manualmente la somma dei ranghi assegnati alle differenze con segno positivo, e la somma dei ranghi assegnati alle differenze con segno negativo, per confrontare questo intervallo con l’intervallo tabulato sulle tavole di Wilcoxon per campioni appaiati, e confermare la nostra decisione statistica. Ecco come fare per calcolare le due somme.

> a = c(214, 159, 169, 202, 103, 119, 200, 109, 132, 142, 194, 104, 219, 119, 234)
> b = c(159, 135, 141, 101, 102, 168, 62, 167, 174, 159, 66, 118, 181, 171, 112)
>
> diff = c(a - b)
#calcolo il vettore contenente le differenze tra tutti i valori di a e i valori di b
> diff.rank = rank(abs(diff)) #assegno i ranghi alle differenze calcolate prima, prese in valore assoluto; la funzione abs(x) dà in output il valore di x senza segno
> diff.rank.sign = diff.rank * sign(diff) #do il segno ai ranghi, richiamando i segni dei valori delle differenze
> ranghi.pos = sum(diff.rank.sign[diff.rank.sign > 0]) #calcolo la somma dei ranghi assegnati alle differenze con segno positivo, ossia maggiori di zero
> ranghi.neg = -sum(diff.rank.sign[diff.rank.sign < 0]) #calcolo la somma dei ranghi assegnati alle differenze con segno negativo, ossia minori di zero, cambiando di segno il risultato, che sarebbe logicamente negativo
> ranghi.pos #chiamo il risultato: esso corrisponde al valore V calcolato con il wilcox.test
[1] 80
> ranghi.neg #chiamo il risultato
[1] 40


L’intervallo calcolato è (40, 80). L’intervallo tabulato sulle tavole di Wilcoxon per campioni appaiati, con 15 differenze è (25, 95). Poiché l’intervallo calcolato è contenuto nell’intervallo tabulato, accettiamo l’ipotesi H0 di uguaglianza delle medie. Come preannunciato dal p-value, la chiusura delle strade al traffico non ha apportato alcun miglioramento in termini di tasso d’inquinamento.

NOTA.
Il test di Wilcoxon dei ranghi con segno prevede la NON assegnazione del rango a quelle differenze pari a zero. Pertanto con il metodo illustrato, si avrebbero risultati falsati, in quanto anche a differenze pari a zero, viene assegnato un rango. Si consiglia pertanto di verificare sempre i contenuti dei vari vettori che vengono calcolati. Se qualcuno ha tempo, può provare a costruire una funzione con le sequenza di calcoli da me utilizzata, aggiungendo la condizione di eliminare dall'assegnazione dei ranghi quelle differenze pari a zero. Io per ora non ho tempo, ma magari in futuro scriverò una funzione completamente automatizzata.

EDIT.
La risposta al dubbio espresso nella nota è molto più semplice di quanto avessi previsto! Per eliminare le differenze pari a zero dal calcolo, è sufficiente aggiungere la seguente istruzione:
diff <- diff[ diff != 0 ]
e il gioco è fatto!

Confronto tra due proporzioni: metodo parametrico e non parametrico

Verifica di ipotesi tra due proporzioni: metodo parametrico (test Z su proporzioni) e non parametrico (test del chi quadro).

È possibile effettuare la verifica di ipotesi anche su dati che seguono una distribuzione binomiale. Consideriamo ad esempio il seguente problema.
Il proprietario di una società di scommesse vuole verificare se un suo clienti sta barando oppure no. Per fare questo vuole confrontare il numero di successi del suo cliente col numero di successi di un suo dipendente, del quale ha la certezza che non bara. In un mese di tempo, il suo dipendente effettua 74 scommesse e ne vince 30; il suo cliente nello stesso arco di tempo effettua 103 scommesse, vincendone 65. Il suo cliente è un baro, oppure no?


Un problema del genere può essere risolto in due diversi modi: utilizzando un metodo parametrico, e uno non parametrico.

  • Risoluzione col metodo parametrico: test z.

Si ricorre a un test z nel caso in cui si possano fare le seguenti due assunzioni: la probabilità comune di successo è approssimabile a 0.5, e il numero di giocate è molto elevato (assunzioni indispensabili per poter approssimare una binomiale ad una Gauss). Supponiamo che questo sia il nostro caso. In R non esiste una funzione per calcolare il valore di z; perciò riprendiamo la formula matematica, e scriviamo la nostra funzione:



> z.prop = function(x1,x2,n1,n2){
+ numeratore = (x1/n1) - (x2/n2)
+ p.comune = (x1+x2) / (n1+n2)
+ denominatore = sqrt(p.comune * (1-p.comune) * (1/n1 + 1/n2))
+ z.prop.ris = numeratore / denominatore
+ return(z.prop.ris)
+ }


La funzione z.prop calcola il valore di z, ricevendo in input il numero di successi (x1 e x2), e il numero di giocate totale (n1 e n2). Applichiamo la funzione appena scritta coi dati del nostro problema:

> z.prop(30, 65, 74, 103)
[1] -2.969695


Abbiamo ottenuto un valore di z superiore al valore di z-tabulato (1.96), il che ci porta a concludere che il cliente che il direttore stava controllando è effettivamente un baro, in quanto le sue probabilità di successo sono più elevate rispetto ad un utente che non bari al gioco.

  • Risoluzione col metodo non parametrico: il test del chi-quadro.

Supponiamo adesso di non poter fare alcuna assunzione sui dati del problema, e quindi di non poter approssimare la binomiale ad una Gauss. Risolviamo il problema con il test del chi-quadro applicato ad una tabella di contingenza 2x2. In R esiste la funzione prop.test che fa al caso nostro.

> prop.test(x = c(30, 65), n = c(74, 103), correct = FALSE)

2-sample test for equality of proportions without continuity
correction

data: c(30, 65) out of c(74, 103)
X-squared = 8.8191, df = 1, p-value = 0.002981
alternative hypothesis: two.sided
95 percent confidence interval:
-0.37125315 -0.08007196
sample estimates:
prop 1 prop 2
0.4054054 0.6310680


La funzione prop.test calcola il valore del chi-quadro, dati i valori dei successi (contenuti nel vettore x) sul numero totale di tentativi (contenuti nel vettore n). I vattori x ed n possono anche essere dichiarati esternamente, e poi venir richiamati come di consueto: prop.test(x, n, correct=FALSE).

Nel caso di campioni piccoli (basso valore di n), si deve specificare correct=TRUE, che nel calcolo del chi-quadro modifica la formula in base alla continuità di Yates:

> prop.test(x = c(30, 65), n = c(74, 103), correct=TRUE)

2-sample test for equality of proportions with continuity correction

data: c(30, 65) out of c(74, 103)
X-squared = 7.9349, df = 1, p-value = 0.004849
alternative hypothesis: two.sided
95 percent confidence interval:
-0.38286428 -0.06846083
sample estimates:
prop 1 prop 2
0.4054054 0.6310680


In entrambi i casi, abbiamo ottenuto p-value minore di 0.05, che ci porta a rifiutare l’ipotesi di uguaglianza delle probabilità. In conclusione, il cliente è un baro. Per conferma confrontiamo il valore chi-quadro-calcolato con il valore chi-quadro-tabulato, che calcoliamo in questo modo:

> qchisq(0.950, 1)
[1] 3.841459


La funzione qchisq calcola il valore del chi-quadro in funzione di alpha e dei gradi di libertà. Poiché chi-quadro-calcolato è maggiore di chi-quadro-tabulato, concludiamo col rifiutare l’ipotesi H0 (come già affermato dal p-value, e dal test parametrico).

NOTA
Da notare è che il valore di Z-calcolato nel primo z-test è esattamente pari alla radice quadrata del valore del chi-quadro calcolato nel prop.test senza la correzione di Yates. Questo non è certo un caso. Per applicare il z test abbiamo supposto che il numero di prove fosse sufficientemente grande per approssimare la binomiale alla Gauss. Nell'esercizio proposto, questo è vero, e proprio per questa ragione abbiamo la uguaglianza dei valori ottenuti (zeta-quadro = chi-quadro). Se il numero delle prove fosse stato più piccolo, avremmo eseguito il test del chi-quadro con la continuità di Yates, ottenendo come risultato un chi-quadro che non è pari al quadrato di zeta.
In altre parole, quando il numero di prove è sufficientemente elevato, non è necessario effettuare il test del chi-quadro (senza la correzione di Yates), perchè possiamo tranquillamente assumere che la binomiale è approssimabile alla Gauss (e difatti i valori sono uguali), procedendo pertanto con un z-test.
Viceversa, in caso di un numero di prove ridotto, bisogna ricorrere al test del chi-quadro con la correzione di Yates, senza effettuare uno zeta-test.

sabato 15 novembre 2008

Test t di verifica della media di due campioni appaiati

Confronto delle medie di due gruppi di campioni appaiati, estratti da due popolazioni a varianza incognita.

Una scuola di atletica ha assunto un nuovo istruttore, e si vuole testare l’efficacia del nuovo tipo di allenamento proposto confrontando le medie dei tempi di 10 centometristi. Vengono riportati di seguito i tempi in secondi prima e dopo l’allenamento di ciascun atleta.

Prima dell’allenamento: 12.9, 13.5, 12.8, 15.6, 17.2, 19.2, 12.6, 15.3, 14.4, 11.3
Dopo l’allenamento: 12.7, 13.6, 12.0, 15.2, 16.8, 20.0, 12.0, 15.9, 16.0, 11.1

Siamo di fronte a due gruppi di campioni appaiati, in quanto le misurazioni sono state effettuate sugli stessi atleti prima e dopo l’allenamento. Per verificare se c’è stato un miglioramento, un peggioramento, o se le medie dei tempi sono rimaste sostanzialmente uguali (ipotesi H0), dobbiamo effettuare un test t di student per campioni appaiati, che in R si richiama così:

> a = c(12.9, 13.5, 12.8, 15.6, 17.2, 19.2, 12.6, 15.3, 14.4, 11.3)
> b = c(12.7, 13.6, 12.0, 15.2, 16.8, 20.0, 12.0, 15.9, 16.0, 11.1)
>
> t.test(a,b, paired=TRUE)

Paired t-test

data: a and b
t = -0.2133, df = 9, p-value = 0.8358
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.5802549 0.4802549
sample estimates:
mean of the differences
-0.05


Il p-value è maggiore di 0.05, quindi possiamo accettare l’ipotesi H0 di uguaglianza delle medie. In conclusione il nuovo allenamento non ha apportato nessun miglioramento (o peggioramento) significativo alla squadra di atleti.
Per conferma, calcoliamo il valore t-tabulato:

> qt(0.975, 9)
[1] 2.262157


Il valore t-calcolato è minore del valore t-tabulato, e questo ci assicura che possiamo tranquillamente accettare l’ipotesi H0.




Supponiamo adesso che il manager della squadra (alla luce dei dati ottenuti) licenzi questo allenatore che non ha apportato alcun miglioramento e ne assuma un altro, più promettente. Riportiamo i tempi degli atleti a seguito del secondo allenamento:

Prima dell’allenamento: 12.9, 13.5, 12.8, 15.6, 17.2, 19.2, 12.6, 15.3, 14.4, 11.3
Dopo il secondo allenamento: 12.0, 12.2, 11.2, 13.0, 15.0, 15.8, 12.2, 13.4, 12.9, 11.0

Adesso verifichiamo se effettivamente c’è stato un miglioramento, ossia eseguiamo un t-test per dati appaiati, specificando in R di verificare l’ipotesi alternativa H1 di miglioramento dei valori. Per far questo è sufficiente aggiungere la sintassi alt=”less” quando richiamiamo il t-test:

> a = c(12.9, 13.5, 12.8, 15.6, 17.2, 19.2, 12.6, 15.3, 14.4, 11.3)
> b = c(12.0, 12.2, 11.2, 13.0, 15.0, 15.8, 12.2, 13.4, 12.9, 11.0)
>
> t.test(a,b, paired=TRUE, alt="less")

Paired t-test

data: a and b
t = 5.2671, df = 9, p-value = 0.9997
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
-Inf 2.170325
sample estimates:
mean of the differences
1.61


Con questa sintassi abbiamo chiesto a R di verificare se la media dei valori contenuti nel vettore a è minore (less, in inglese) della media dei valori contenuti nel vettore b. In risposta abbiamo ottenuto un p-value di gran lunga superiore a 0.05, e questo ci porta a concludere che possiamo accettare l’ipotesi H1: il nuovo allenamento ha apportato sostanziali miglioramenti alla squadra.

Se avessimo scritto: t.test(a,b, paired=TRUE, alt=”greater”), avremmo chiesto a R di verificare se la media dei valori contenuti nel vettore a è maggiore della media dei valori contenuti nel vettore b. Alla luce del risultato precedente, possiamo sospettare che otterremo un p-value molto minore rispetto a 0.05, e difatti:

> a = c(12.9, 13.5, 12.8, 15.6, 17.2, 19.2, 12.6, 15.3, 14.4, 11.3)
> b = c(12.0, 12.2, 11.2, 13.0, 15.0, 15.8, 12.2, 13.4, 12.9, 11.0)
>
> t.test(a,b, paired=TRUE, alt="greater")

Paired t-test

data: a and b
t = 5.2671, df = 9, p-value = 0.0002579
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
1.049675 Inf
sample estimates:
mean of the differences
1.61


Un valore di p-value inferiore a 0.05 ci porta a escludere l'ultima ipotesi H1 formulata, come del resto ci era stato confermato dal calcolo precedente.

Test di Wilcoxon-Mann-Whitney della somma dei ranghi

Confronto delle medie di due gruppi di campioni indipendenti, in cui non si può supporre una distribuzione di tipo gaussiana; è l’equivalente del test U di Mann-Whitney.

Si vuole verificare la media dei goal subito da due squadra calcistiche nel corso degli anni. Vengono riportati di seguito il numero di goal subiti da ciascuna squadra in 6 partite successive, per ciascun anno (dati inventati).
Squadra A: 6, 8, 2, 4, 4, 5
Squadra B: 7, 10, 4, 3, 5, 6

Il test di Wilcoxon-Mann-Whitney (o test di Wilcoxon della somma dei ranghi) si applica nel caso in cui si chiede di confrontare le medie dei valori di due gruppi che non seguono una distribuzione normale. È l’equivalente del test t per campioni indipendenti.
Vediamo come risolvere il problema con R:

> a = c(6, 8, 2, 4, 4, 5)
> b = c(7, 10, 4, 3, 5, 6)
>
> wilcox.test(a,b, correct=FALSE)

Wilcoxon rank sum test

data: a and b
W = 14, p-value = 0.5174
alternative hypothesis: true location shift is not equal to 0


Il p-value è maggiore di 0.05, quindi possiamo accettare l’ipotesi H0 di uguaglianza significativa delle medie dei due gruppi.
Se avessi eseguito il wilcox.test(b,a, correct=FALSE), il p-value sarebbe stato logicamente lo stesso:

> a = c(6, 8, 2, 4, 4, 5)
> b = c(7, 10, 4, 3, 5, 6)
>
> wilcox.test(b,a, correct=FALSE)

Wilcoxon rank sum test

data: b and a
W = 22, p-value = 0.5174
alternative hypothesis: true location shift is not equal to 0


Il valore W (sui libri di testo è solitamente indicato con U) è così calcolato: dopo aver riunito tutti i valori in un unico gruppo, si dispongono in ordine crescente e si assegnano i ranghi. Fatto ciò, si prendono i ranghi del gruppo a e si sommano. A questo valore si sottrae il risultato che si ottiene dalla seguente formula: Na(Na+1)/2 (dove Na è la numerosità campionaria del gruppo a).
Se eseguiamo il wilcox.test(b,a, correct=FALSE), avremmo ottenuto W = 22, che corrisponde al numero che si ottiene sottraendo alla somma dei ranghi del gruppo b il risultato dell’operazione Nb(Nb+1)/2 (dove Nb è la numerosità campionaria del gruppo b).

Verifichiamo quanto detto:

> somma.ranghi.a = sum(rank(c(a,b))[1:6]) #somma dei ranghi assegnati al gruppo a
> W = somma.ranghi.a – (length(a)*(length(a)+1)) / 2 #calcolo il valore di W
> W
[1] 14

> somma.ranghi.b = sum(rank(c(a,b))[7:12]) #somma dei ranghi assegnati al gruppo b
> W = somma.ranghi.b – (length(b)*(length(b)+1)) / 2 #calcolo il valore di W
> W
[1] 22


Per avere una conferma della nostra decisione statistica (di accettare l’ipotesi H0), possiamo fare riferimento ai valori tabulati sulle tavole di Wilcoxon per campioni indipendenti. L’intervallo tabulato per due gruppi da 6 campioni ciascuno è: (26, 52), mentre i valori che otteniamo dal nostro esercizio sono:

> sum(rank(c(a,b))[1:6]) #somma dei ranghi del gruppo a
[1] 35
> sum(rank(c(a,b))[7:12]) #somma dei ranghi del gruppo b
[1] 43


Poiché l’intervallo calcolato (35, 43) è interno all’intervallo tabulato (26, 52), concludiamo accettando l’ipotesi H0 di uguaglianza delle medie (come del resto ci diceva il p-value).

Test t a due campioni #2

Confronto delle medie di due gruppi di campioni indipendenti, estratti da due popolazioni a varianza incognita; varianze campionarie non omogenee.

Si chiede di confrontare le medie delle altezze di due gruppi. I dati vengono riportati qui di seguito (i valori sono completamente inventati).

A: 175, 168, 168, 190, 156, 181, 182, 175, 174, 179
B: 120, 180, 125, 188, 130, 190, 110, 185, 112, 188

Per risolvere questo problema dobbiamo ricorrere ad un test t di student a due campioni (come vedremo un po’ diverso dal precedente), supponendo che i due campioni siano estratti da popolazioni che seguano una distribuzione di tipo gaussiana (nel caso in cui non si possa supporre ciò, si risolve questo problema sfruttando il metodo non parametrico, chiamato test di Wilcoxon-Mann-Whitney). Prima di procedere con il t-test, è necessario valutare le varianze campionarie dei due gruppi, ossia effettuare un test F di Fisher per verificare l’omoschedasticità (omogeneità delle varianze). In R si procede così:

> a = c(175, 168, 168, 190, 156, 181, 182, 175, 174, 179)
> b = c(120, 180, 125, 188, 130, 190, 110, 185, 112, 188)
>
> var.test(b,a)

F test to compare two variances

data: b and a
F = 14.6431, num df = 9, denom df = 9, p-value = 0.0004636
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
3.637133 58.952936
sample estimates:
ratio of variances
14.64308


Abbiamo ottenuto p-value minore di 0.05, quindi le due varianze non sono omogenee. Difatti possiamo confrontare il valore di F ottenuto con il valore di F tabulato per alpha = 0.05, gradi di libertà del numeratore = 9, e gradi di libertà del denominatore = 9, utilizzando la funzione qf(p, df.num, df.den):

> qf(0.95, 9, 9)
[1] 3.178893


Notiamo che il valore di F-calcolato è maggiore del valore di F-tabulato, il che ci porta a rifiutare l’ipotesi di omogeneità delle varianze.
NOTA: come potete vedere, non ho richiamato la funzione var.test(a,b), ma la funzione var.test(b,a), perchè per convenzione per calcolare il valore di F di Fisher si esegue il rapporto tra le varianze campionarie, con al numeratore la varianza maggiore (al fine di ottenere un valore di F sempre maggiore di 1). Se avessi calcolato var.test(a,b) avrei ottenuto un valore di F minore di 1 (0.002 all'incirca): ottenere un risultato del genere ci obbliga a ripetere il test di Fisher invertendo i parametri, così da ottenere il valore di F maggiore di 1, che può essere confrontato con il valore F-tabulato. Se invece basiamo la nostra decisione solo sul valore di p-value, allora non occorre badare a questo dettaglio, perchè il p-value di var.test(a,b) è uguale al p-value di var.test(b,a).

Richiamiamo quindi la funzione t.test per varianze non omogenee (var.equal=FALSE, che si può anche omettere, perché di default la funzione lavora su varianza non omogenee) e campioni indipendenti (paired=FALSE, che si può anche omettere perché di default la funzione lavora su campioni indipendenti) in questo modo:

> t.test(a,b, var.equal=FALSE, paired=FALSE)

Welch Two Sample t-test

data: a and b
t = 1.8827, df = 10.224, p-value = 0.08848
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.95955 47.95955
sample estimates:
mean of x mean of y
174.8 152.8


Come vediamo nella intestazione, è stato eseguito un t-test su due campioni con il calcolo dei gradi di libertà secondo la formula di Welch-Satterthwaite (il risultato della formula applicata lo possiamo leggere in df = 10.224), cui si ricorre nel caso in cui le varianze non siano omogenee (si applica la formula di Welch-Satterwhaite per il calcolo dei gradi di libertà in k gruppi a varianze non omogenee; nel caso k = 2, come questo, la formula prende anche il nome di formula di Dixon-Massey).
Abbiamo ottenuto p-value maggiore di 0.05, quindi possiamo concludere che le medie dei due gruppi sono significativamente simili (seppur di poco, visto che p-value è molto vicino alla soglia limite 0.05). Difatti il valore di t è minore rispetto al valore t-tabulato per 10.224 gradi di libertà, che in R possiamo calcolare così:

> qt(0.975, 10.224)
[1] 2.221539


Questo ci conferma che possiamo accettare l’ipotesi H0 di uguaglianza delle medie.

Per completezza ecco le formule per il calcolo dei gradi di libertà:

Welch-Satterthwaite equation:



Dixon-Massey equation:


Test t a due campioni #1

Confronto delle medie di due gruppi di campioni indipendenti, estratti da due popolazioni a varianza incognita, varianze campionarie omogenee.

Riprendiamo i dati utilizzati in un precedente esercizio
Si chiede di confrontare le medie delle altezze di due gruppi, estratti da due popolazioni a varianza non nota. I dati vengono riportati qui di seguito (i valori sono completamente inventati).

A: 175, 168, 168, 190, 156, 181, 182, 175, 174, 179
B: 185, 169, 173, 173, 188, 186, 175, 174, 179, 180

Per risolvere questo problema dobbiamo ricorrere ad un test t di student a due campioni, supponendo che i due campioni siano estratti da popolazioni che seguano una distribuzione di tipo gaussiana (nel caso in cui non si possa supporre ciò, si risolve questo problema sfruttando il metodo non parametrico, chiamato test di Wilcoxon-Mann-Whitney). Prima di procedere con il test t, è necessario valutare le varianze campionarie dei due gruppi, ossia effettuare un test F di Fisher per verificare l’omoschedasticità (omogeneità delle varianze). In R si procede così:

> a = c(175, 168, 168, 190, 156, 181, 182, 175, 174, 179)
> b = c(185, 169, 173, 173, 188, 186, 175, 174, 179, 180)
>
> var.test(a,b)

F test to compare two variances

data: a and b
F = 2.1028, num df = 9, denom df = 9, p-value = 0.2834
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.5223017 8.4657950
sample estimates:
ratio of variances
2.102784


Abbiamo ottenuto p-value maggiore di 0.05, quindi possiamo supporre che le due varianze siano omogenee. Difatti possiamo confrontare il valore di F ottenuto con il valore di F tabulato per alpha = 0.05, gradi di libertà del numeratore = 9, e gradi di libertà del denominatore = 9, utilizzando la funzione qf(p, df.num, df.den):

> qf(0.95, 9, 9)
[1] 3.178893


Notiamo che il valore di F calcolato è minore del valore di F tabulato, il che ci porta ad accettare l’ipotesi di omogeneità delle varianze.
NOTA: la distribuzione F ha una sola coda, pertanto con un grado di confidenza del 95%, inseriamo nella funzione di R un valore p = 0.95. Viceversa la distribuzione t ha due code, e per questo nella funzione di R qt(p, df) inseriamo un valore p = 0.975.

Richiamiamo quindi la funzione t.test per varianze omogenee (var.equal=TRUE)e campioni indipendenti (paired=FALSE, che si può anche omettere perché di default la funzione lavora su campioni indipendenti) in questo modo:

> t.test(a,b, var.equal=TRUE, paired=FALSE)

Two Sample t-test

data: a and b
t = -0.9474, df = 18, p-value = 0.356
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-10.93994 4.13994
sample estimates:
mean of x mean of y
174.8 178.2


Abbiamo ottenuto p-value maggiore di 0.05, quindi possiamo concludere che le medie dei due gruppi sono significativamente simili. Difatti il valore di t è minore rispetto al valore t-tabulato per 18 gradi di libertà, che in R possiamo calcolare così:

> qt(0.975, 18)
[1] 2.100922


Questo ci conferma che possiamo accettare l’ipotesi H0 di uguaglianza delle medie.

venerdì 14 novembre 2008

Test z a due campioni

Confronto delle medie di due gruppi di campioni indipendenti, estratti da due popolazioni a varianza nota.

Si chiede di confrontare le medie delle altezze di due gruppi. Il primo gruppo (A) è costituito da individui di nazionalità italiana (varianza della popolazione italiana = 5), il secondo gruppo è estratto da individui di nazionalità tedesca (varianza della popolazione tedesca = 8.5). I dati vengono riportati qui di seguito (i valori sono completamente inventati).

A: 175, 168, 168, 190, 156, 181, 182, 175, 174, 179
B: 185, 169, 173, 173, 188, 186, 175, 174, 179, 180

Poiché disponiamo delle varianze della popolazione, dobbiamo procedere con un test z. Anche in questo caso non è disponibile una funzione in R, ma possiamo facilmente crearla noi.

> z.test2 = function(a, b, var.a, var.b){
+ n.a = length(a) #numerosità campionaria del gruppo a
+ n.b = length(b) #numerosità campionaria del gruppo b
+ zeta = (mean(a) - mean(b)) / (sqrt(var.a/n.a + var.b/n.b))
+ return(zeta)
+ }


La funzione z.test2 fornisce in output il valore di zeta, dopo aver ricevuto in input le variabili (a e b), la varianza della prima popolazione (var.a) e la varianza della seconda popolazione (var.b).
Utilizzando questa funzione otteniamo:

> a = c(175, 168, 168, 190, 156, 181, 182, 175, 174, 179)
> b = c(185, 169, 173, 173, 188, 186, 175, 174, 179, 180)
>
> z.test2(a, b, 5, 8.5)
[1] -2.926254


Il valore di zeta è maggiore rispetto al valore di zeta tabulato per alpha pari a 0.05 (z-tabulato = 1.96). Concludiamo che le due medie sono significativamente diverse.

Test t di student a un solo campione

Confronto della media campionaria con un valore noto, varianza della popolazione non nota.

Riprendiamo la traccia di un esercizio precedente.
È stato effettuato un test di intelligenza su 10 soggetti, e di seguito vengono riportati i risultati ottenuti. Il risultato medio della popolazione ottenuto allo stesso test è pari a 75. Si vuole verificare se la media campionaria è significativamente simile (si ponga il livello di significatività pari al 5%) alla media della popolazione, supponendo che la varianza della popolazione non sia nota.

65, 78, 88, 55, 48, 95, 66, 57, 79, 81

Contrariamente a quanto fatto per il test z a un solo campione, per il test t di student a un solo campione esiste una funzione preimpostata in R che possiamo subito applicare.
Si tratta della funzione t.test(a, mu), che possiamo vedere qui di seguito applicata.

> a = c(65, 78, 88, 55, 48, 95, 66, 57, 79, 81)
>
> t.test (a, mu=75)

One Sample t-test

data: a
t = -0.783, df = 9, p-value = 0.4537
alternative hypothesis: true mean is not equal to 75
95 percent confidence interval:
60.22187 82.17813
sample estimates:
mean of x
71.2


La funzione t.test a un solo campione (one sample in inglese) ci fornisce in output il valore di t calcolato; ci dà inoltre i gradi di libertà (df = degree freedom), l’intervallo di confidenza (confidence interval) e la media (mean of x).
Per giungere alla conclusione statistica possiamo procedere in due modi. Possiamo ad esempio confrontare il valore di t con il valore tabulato della t di student con 9 gradi di libertà. Se non disponiamo delle tavole, possiamo calcolare il valore t-tabulato in questo modo:

> qt(0.975, 9)
[1] 2.262157


La funzione qt(p, df) restituisce il valore di t in base alla significatività scelta (abbiamo scelto un livello di significatività pari al 95%, il che significa che in ogni coda abbiamo il 2.5%, che corrisponde al valore di p inserito: 1 – 0.025), e ai gradi di libertà. Confrontando il valore di t calcolato nel t test con il t tabulato risulta minore, il che significa che accettiamo l’ipotesi nulla di uguaglianza delle medie: la mia campionaria è significativamente simile alla media della popolazione.

In alternativa avremmo potuto valutare il valore p-value. Con un livello di significatività del 5%, ricordiamo questa regola: se p-value è maggiore di 0.05 allora accettiamo l’ipotesi H0; se p-value è minore di 0.05 allora rifiutiamo l’ipotesi H0.

Test z a un solo campione

Confronto della media campionaria con un valore noto, varianza nota.

È stato effettuato un test di intelligenza su 10 soggetti, e di seguito vengono riportati i risultati ottenuti. Il risultato medio della popolazione ottenuto allo stesso test è pari a 75. Si vuole verificare se la media campionaria è significativamente simile (si ponga il livello di significatività pari al 5%) alla media della popolazione, supponendo che la varianza sia nota e pari a 18.

65, 78, 88, 55, 48, 95, 66, 57, 79, 81

Per risolvere questo problema è necessario sviluppare un test z a un solo campione. In R non esiste una funzione del genere, per questo motivo possiamo scegliere di eseguire i calcoli singolarmente, oppure di creare una funzione apposita.
Ricordando la formula per calcolare il valore di z, possiamo procedere scrivendo questa funzione:

> z.test = function(a, mu, var){
+ zeta = (mean(a) - mu) / (sqrt(var / length(a)))
+ return(zeta)
+ }


Abbiamo costruito così la funzione z.test che ricevendo in input la sequenza di dati (a), la media della popolazione con cui eseguire il confronto (mu), e la varianza della popolazione (var), restituisce il valore di zeta. Applichiamo ora la funzione al nostro problema.

> a = c(65, 78, 88, 55, 48, 95, 66, 57, 79, 81)
>
> z.test(a, 75, 18)
[1] -2.832353


Il valore di zeta è pari a –2.83, che è maggiore (in valore assoluto) rispetto al valore di zeta tabulato per alfa pari a 0.05 (zeta-tabulato = 1.96). Concludiamo quindi che la media del nostro campione è significativamente diversa dalla media della popolazione.

giovedì 13 novembre 2008

Media geometrica e media armonica in R

Calcolare la media geometrica e la media armonica della sequenza data.

10, 2, 19, 24, 6, 23, 47, 24, 54, 77

Queste funzioni non sono presenti nel pacchetto standard di R, sebbene siano facilmente reperibili in rete pacchetti di funzioni aggiuntive. Ad ogni modo, è facile calcolare questi valori semplicemente riprendendo le formule matematiche, applicandole in R.


Tradotto in R questa formula diventa:

> a = c(10, 2, 19, 24, 6, 23, 47, 24, 54, 77)
>
> 1/mean(1/a) #calcola la media armonica
[1] 10.01109



Tradotto in R questa formula diventa:

> a = c(10, 2, 19, 24, 6, 23, 47, 24, 54, 77)
>
> n = length(a) #assegno alla variabile n un valore pari al numero di elementi in a
>
> prod(a)^(1/n) #calcolo la media geometrica
[1] 18.92809

Esercizio sulla probabilità: distribuzione binomiale negativa

Qual è la probabilità che lanciando una moneta si ottenga la 4 croce prima della 3 testa?


La formula matematica per risolvere questo esercizio, che segue una distribuzione binomiale negativa, è la seguente:

Per risolvere un problema del genere, il numero di lanci non è significativo, in quanto è sufficiente conoscere la probabilità dei singoli lanci che si effettuano e che si vogliono combinare, per ottenere la probabilità richiesta. Per risolvere questo tipo di problema in R si usa la funzione dnbinom(x,y,p). Questa distribuzione permette di calcolare la probabilità che un numero di fallimenti x avvenga prima del successo y in una sequenza di prove bernoulliane per la quali la probabilità del singolo successo è p.

> dnbinom(4,3,0.5)
[1] 0.1171875


La probabilità che esca la quarta croce prima dell’uscita della terza testa (quindi dopo che sono uscite due teste) è pari all’11,72%.

Può sembrare un risultato strano, ma per convincerci dell’esattezza di questa funzione di R, proviamo a considerare quest’altro problema: quali sono le probabilità che esca la prima, la seconda, la terza, ... , la 25esima testa prima della seconda croce?
Possiamo ottenere un istogramma in R, che mostra quanto richiesto, in questa maniera:

> barplot(dnbinom(1:25,2,0.5), col="grey", names.arg=1:25)

Osservate che la probabilità che esca la prima testa, prima della seconda croce è 0.25, ossia il prodotto 0.5 * 0.5 (probabilità che esca T dopo che è uscito C). Man mano che i lanci della moneta vanno avanti, la probabilità scende sempre più.

Esercizio sulla probabilità: distribuzione di Bernoulli

Qual è la probabilità, lanciando una moneta 8 volte, di ottenere la sequenza TTCCCTCC?


La teoria ci insegna che per risolvere questo quesito è sufficiente utilizzare la seguente formula:

Per un problema del genere, in R esiste la funzione dbinom(x,n,p). Ci troviamo nel caso di una distribuzione binomiale, in cui gli eventi possono essere T oppure C. Supponiamo che C sia il successo x (in questo caso 5), mentre n è il numero di tentativi (in questo caso 8). La probabilità del successo è p = 0.5. Immettiamo questi dati in R e otterremo la risposta:

> dbinom(5, 8, 0.5)
[1] 0.21875


La probabilità di ottenere quella specifica sequenza è pari al 21.875%.
Cosa sarebbe accaduto se avessimo scelto come successo l’uscita della testa (ossia imponendo x=3)?

mercoledì 12 novembre 2008

Facciamo pratica con le funzioni di R

Data la seguente serie di dati, calcolare la media aritmetica, la mediana, la varianza, la deviazione standard, il valore più grande, il valore più piccolo, la somma di tutti i valori, il quadrato della somma di tutti i valori, la somma del quadrato di tutti i valori, assegnare i ranghi e sommarli, assegnare i ranghi e sommare quelli relativi ai primi 6 valori.

10, 2, 19, 24, 6, 23, 47, 24, 54, 77

Esercizio estremamente semplice da eseguire in R, volutamente creato per prendere dimestichezza con le funzioni statistiche che ci potranno servire in futuro. Non riporto le formule matematiche, che sono reperibili su qualsiasi libro di matematica, o su Wikipedia.
Vediamo una per una questa funzioni.

> a = c(10, 2, 19, 24, 6, 23, 47, 24, 54, 77)
>
> mean(a) #calcola la media aritmetica
[1] 28.6
>
> median(a) #calcola la mediana
[1] 23.5
>
> var(a) #calcola la varianza
[1] 561.8222
>
> sd(a) #calcola la deviazione standard
[1] 23.70279
>
> max(a) #riporta il valore più grande della sequenza
[1] 77
>
> min(a) #riporta il valore più piccolo della sequenza
[1] 2
>
> sum(a) #somma tutti i valori
[1] 286
>
> sum(a)*sum(a) #calcola il quadrato della somma di tutti i valori
[1] 81796
>
> sum(a*a) #calcola la somma dei quadrati di tutti i valori
[1] 13236
>
> sum(rank(a)) #somma dei ranghi assegnati alle variabili contenute in a
[1] 55
>
> sum(rank(a)[1:6]) #somma dei ranghi assegnati ai primi 6 valori di a
[1] 21.5

sabato 1 novembre 2008

Usefull links & more stuffs

Link utili


Siti web e Blogs che parlano di R


Libri consigliati

  • Handbook of Parametric and Nonparametric Statistical Procedures, Second Edition, by David J. Sheskin


  • Nonparametric Statistics for Non-Statisticians: A Step-by-Step Approach, by Gregory W. Corder


  • R in Action - Data Analysis and Graphics with R, by Robert I. Kabacoff


Indice degli esercizi svolti

Capitolo -1: Come funziona R



Capitolo 0: Esercizi base con R e probabilità



Capitolo 1: Inferenza statistica, verifica di ipotesi su uno o più gruppi



Capitolo 2: Studio della correlazione e della regressione



Capitolo 3: Analisi del trend (trend analysis)



Capitolo 4: Blocking design in R



Capitolo 5: Grafici con R



Capitolo 6: Approfondimenti



Capitolo 7: Cluster Analysis



Capitolo 8: Matrici e algebra lineare e matematica in generale



Capitolo 9: Non solo statistica


Proponi un articolo

I capitoli della statistica sono praticamente infiniti, essendo incredibilmente vasto il campo di applicazione delle indagini.
Ogni giorno quindi sono assalito dal dubbio: "di quale tecnica parlare oggi? Continuo con l'inferenza statistica? Scrivo un post sull'analisi del trend? Un tutorial sui grafici con R?". E i miei dubbi si fanno ancora più atroci non sapendo cosa realmente possa interessare a chi legge il mio blog.

Pur avendo in lista d'attesa numerosi capitoli di cui ancora parlare, vorrei sentire le vostre proposte. Utilizzate lo spazio dei commenti in fondo a questo post per proporre un argomento, ed io sarò ben lieto di approfondire il discorso che voi proporrete (sempre nei limiti delle mie conoscenze, s'intende :-) ).