Tabulare una serie di test di normalità (gaussianità)

Quando con R si esegue una serie di test di normalità (gaussianità) su una variabile, può essere fastidioso consultare e confrontare tra loro risultati che nella Console di R sono riportati separatamente e lontani l'uno dall'altro. Personalmente trovo che gli stessi risultati compattati in una tabella e incolonnati siano molto meglio fruibili. Ma non solo. I valori di probabilità p incolonnati sono poi ancor meglio confrontabili tra loro se riportati in formato fisso anziché in formato esponenziale.    

Vediamo quindi come incolonnare e riportare in formato fisso i risultati dei più comuni test di gaussianità con lo script che segue. Accertatevi di avere installati i pacchetti moments e nortest. Altrimenti dovete scaricarli e installarli dal CRAN. Accertatevi inoltre di avere installato il pacchetto DAAG che contiene i dati impiegati nello script, o in alternativa procedete come indicato [1].

Copiate lo script, incollatelo nella Console di R e premete ↵ Invio.

# TABULARE UNA SERIE DI TEST DI NORMALITA' (GAUSSIANITA')

#

library(DAAG) # carica il pacchetto DAAG incluso il set di dati ais

mydata <- ais[,c(5)] # salva in mydata la colonna 5 con la variabile ferritina

#

# asimmetria, curtosi e quattro test di gaussianità

#

library(moments) # carica il pacchetto

agostino.test(mydata) # test di D'Agostino per il coefficiente di asimmetria

anscombe.test(mydata) # test di Anscombe-Glynn per il coefficiente di curtosi

library(nortest) # carica il pacchetto

ad.test(mydata) # test di Anderson-Darling per la gaussianità

cvm.test(mydata) # test di Cramer-von Mises per la gaussianità

pearson.test(mydata) # test chi-quadrato di Pearson per la gaussianità

sf.test(mydata) # test di Shapiro-Francia per la gaussianità

#

# sintetizza in una tabella i test con i risultati e i valori di p

#

stat <- c(names(agostino.test(mydata)$statistic[2]), names(anscombe.test(mydata)$statistic[2]), names(ad.test(mydata)$statistic), names(cvm.test(mydata)$statistic), names(pearson.test(mydata)$statistic), names(sf.test(mydata)$statistic)) # array con le denominazioni delle statistiche calcolate

#

ris <- c(agostino.test(mydata)$statistic[2], anscombe.test(mydata)$statistic[2], ad.test(mydata)$statistic, cvm.test(mydata)$statistic, pearson.test(mydata)$statistic, sf.test(mydata)$statistic) # array con i risultati delle statistiche calcolate

#

pval <- c(agostino.test(mydata)$p.value, anscombe.test(mydata)$p.value, ad.test(mydata)$p.value, cvm.test(mydata)$p.value, pearson.test(mydata)$p.value, sf.test(mydata)$p.value) # array con i valori di p

#

mydataset <- data.frame(stat, ris, pval) # combina gli array in una tabella

#

rownames(mydataset) <- c(agostino.test(mydata)$method, anscombe.test(mydata)$method, ad.test(mydata)$method, cvm.test(mydata)$method, pearson.test(mydata)$method, sf.test(mydata)$method) # aggiunge i nomi delle righe

colnames(mydataset) <- c("Statistica", "Risultato", "Valore p") # aggiunge i nomi delle colonne

mydataset # mostra la tabella

options(scipen=999) # esprime i numeri in formato fisso

mydataset # mostra nuovamente la tabella

options(scipen=0) # ripristina la notazione scientifica/esponenziale

#

Il test di asimmetria di D'Agostino agostino.test() valuta quanto e come la simmetria della distribuzione campionaria (quella effettivamente osservata) si discosta da quella della distribuzione gaussiana teorica, che è per definizione perfettamente simmetrica. Il test di curtosi di Anscombe-Glynn anscombe.test() valuta se la distribuzione campionaria è troppo appiattita (bassa e larga) o troppo appuntita (alta e stretta) rispetto alla distribuzione gaussiana teorica, in ciascun punto della quale vale un preciso rapporto tra distanza dal centro della distribuzione e altezza della curva.

Mentre i due test precedenti valutano separatamente asimmetria e curtosi, gli altri quattro test qui impiegati, il test di Anderson-Darling ad.test(), il test di Cramer-von Mises cvm.test(), il test chi-quadrato di Pearson pearson.test() e il test di Shapiro-Francia sf.test(), sono test globali di gaussianità (forniscono una valutazione globale del grado di scostamento della distribuzione osservata dalla distribuzione gaussiana teorica, ma non del tipo di scostamento).

Dopo avere calcolato e riportato separatamente i risultati dei sei test, viene generata la tabella che mostra per ciascun test la denominazione, il simbolo della statistica calcolata, il risultato numerico di detta statistica e il corrispondente valore di probabilità p procedendo in questo modo:

→ per ciascun test viene riportato nell'array (o se preferite vettore) stat il simbolo della statistica calcolata names();

→ per ciascun test viene riportato nell'array ris il risultato numerico della statistica calcolata $statistic;

→ per ciascun test viene riportato nell'array pval il valore di probabilità p del risultato numerico $p.value;

→ mediante la funzione data.frame() gli array stat, ris e pval sono combinati nella tabella mydataset;
→ mediante la funzione rownames() a ciascuna riga della tabella mydataset viene assegnata la denominazione del test corrispondente $method;

→ mediante la funzione colnames() sono assegnati i nomi alle colonne della tabella mydataset che contengono, nell'ordine, il simbolo stat della statistica (Statistica), il risultato ris della statistica (Risultato) e il valore pval di probabilità p (Valore p) [2].

La tabella viene quindi mostrata nella Console di R prima con i valori numerici espressi in notazione scientifica

                                  Statistica  Risultato     Valore p

D'Agostino skewness test                   z  6.1376120 8.377114e-10

Anscombe-Glynn kurtosis test               z  2.9448982 3.230609e-03

Anderson-Darling normality test            A  6.0970351 4.903750e-15

Cramer-von Mises normality test            W  0.9447263 2.495032e-09

Pearson chi-square normality test          P 63.7722772 2.530969e-08

Shapiro-Francia normality test             W  0.8913765 1.251172e-09

e poi con i valori numerici espressi in notazione fissa mediante options(scipen=999)

                                  Statistica  Risultato               Valore p

D'Agostino skewness test                   z  6.1376120 0.00000000083771142715

Anscombe-Glynn kurtosis test               z  2.9448982 0.00323060945215406870

Anderson-Darling normality test            A  6.0970351 0.00000000000000490375

Cramer-von Mises normality test            W  0.9447263 0.00000000249503195167

Pearson chi-square normality test          P 63.7722772 0.00000002530968668914

Shapiro-Francia normality test             W  0.8913765 0.00000000125117173937

L'incolonnamento e il formato fisso rendono sicuramente più chiara la sintesi e la valutazione comparativa dei risultati numerici. La probabilità p è la probabilità di osservare per caso il risultato numerico della statistica (z, A, W, P secondo i casi) che misura lo scostamento della distribuzione osservata dalla distribuzione gaussiana teorica. Se tale probabilità è sufficientemente piccola, si conclude che il risultato ottenuto non è imputabile al caso e che pertanto la distribuzione dei dati si discosta significativamente dalla distribuzione gaussiana. Tutti e sei i test confermano che la ferritina non segue una distribuzione gaussiana: infatti abbiamo un valore p < 0.05 - cioè un p inferiore al valore del 5% tradizionalmente assunto come soglia per la significatività - in tutti i casi, e addirittura di molti ordini di grandezza inferiore al 5% in cinque casi su sei. Potete avere una conferma visiva di queste conclusioni con le rappresentazioni grafiche, sempre altamente consigliate come complemento all'analisi statistica, che potete realizzare con lo script riportato nel post Analizzare graficamente la distribuzione di una variabile. 

Come al solito lo script è realizzato in modo da essere facilmente riutilizzabile, essendo sufficiente a questo scopo:

→ sostituire la prima riga di codice con il codice per importare i vostri dati [3];

→ sostituire nella seconda riga di codice ais[,c(5)] con il nome della variabile da analizzare contenuta nei dati che avete importato.

----------

[1] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG.

[2] Potete al bisogno approfondire la documentazione delle funzioni qui impiegate digitando help(nomedellafunzione) nella Console di R.

[3] Per una guida rapida all'importazione dei dati potete consultare i post:

→ importazione dei dati da un file .csv

→ importazione dei dati da un file .xls o .xlsx

→ gestione dei dati mancanti

Test parametrici e non parametrici per dati appaiati

Il confronto tra dati appaiati si applica quando la stessa variabile viene misurata nello stesso caso [1] in due occasioni diverse. In campo medico la situazione tipica è quella di una misura che viene effettuata sullo stesso soggetto prima e dopo uno specifico trattamento.

Il Volume Espiratorio Massimo nel primo Secondo (VEMS), impiegato nella diagnostica della capacità respiratoria, è il volume di aria espirata nel corso del primo secondo di una espirazione massima forzata e indica il grado di pervietà delle grandi vie aeree. Viene anche denominato FEV1 dall'acronimo inglese di Forced Expiratory Volume in the 1st second.

I seguenti dati, ricavati da Campbell [2], riportano i valori di VEMS (FEV1) misurati in 5 soggetti asmatici prima (al tempo t0) e dopo (al tempo t1) l'assunzione di un broncodilatatore, e la loro differenza (t0 - t1).

t0	t1	Differenza
1.5	1.7	-0.2
1.7	1.9	-0.2
2.1	2.2	-0.1
1.6	1.9	-0.3
2.4	2.4	0

Per proseguire ora è necessario:

→ effettuare il download del file FEV1.csv 

→ salvare il file nella cartella C:\Rdati\

Per il file di dati vedere istruzioni e link riportati alla pagina Dati, ma potete anche semplicemente copiare i dati riportati qui sotto aggiungendo un ↵ Invio al termine dell'ultima riga e salvarli in C:\Rdati\ in un file di testo denominato FEV1.csv (assicuratevi che il file sia effettivamente salvato con l'estensione .csv). 

t0;t1;Differenza

1.5;1.7;-0.2

1.7;1.9;-0.2

2.1;2.2;-0.1

1.6;1.9;-0.3

2.4;2.4;0

Copiate e incollate nella Console di R questo script e premete ↵ Invio:

# CONFRONTO TRA DUE CAMPIONI PER DATI APPAIATI

#

mydata <- read.table("c:/Rdati/FEV1.csv", header=TRUE, sep=";") # importa i dati

attach(mydata) # consente di impiegare direttamente i nomi delle variabili

mydata # mostra i dati

#

library(nortest) # carica il pacchetto

lillie.test(Differenza) # test di normalità di Lilliefors (Kolmogorov-Smironv) sulla Differenza (t0 - t1)

#

t.test(t0, t1, paired=TRUE) # test t di Student per dati appaiati (test parametrico)

#

wilcox.test(t0, t1, paired=TRUE, exact=FALSE) # test di Wilcoxon per dati appaiati (test non parametrico)

#

detach(mydata) # termina l'impiego diretto dei nomi delle variabili

#

Dopo avere importato i dati del file nell'oggetto mydata e dopo avere eseguito la funzione attach() che consente nel codice successivo di impiegare direttamente i nomi delle variabili, con mydata sono mostrati i dati importati:

> mydata # mostra i dati 

     t0   t1 Differenza

1   1.5  1.7       -0.2

2   1.7  1.9       -0.2

3   2.1  2.2       -0.1

4   1.6  1.9       -0.3

5   2.4  2.4        0.0

La scelta tra il test parametrico, il test t di Student, e il test non parametrico, il test di Wilcoxon (o Wilcoxon Signed Rank Test) viene effettuata come al solito sulla base dei risultati dei test di normalità (gaussianità).

Qui viene impiegato come test di normalità (gaussianità) il test di Lilliefors (Kolmogorov-Smirnov) che viene applicato alla variabile Differenza (t0 - t1) presa con il segno:

> lillie.test(mydata$Differenza) # test di normalità di Lilliefors (Kolmogorov-Smironv) sulla Differenza (t0 - t1)

        Lilliefors (Kolmogorov-Smirnov) normality test

data:  Differenza

D = 0.23714, p-value = 0.4686

Il test conferma una distribuzione gaussiana dei dati (p = 0.4686). Si procede quindi a impiegare le conclusioni tratte dal test t di Student

> t.test(t0, t1, paired=TRUE) # test t di Student per dati appaiati

        Paired t-test

data:  t0 and t1

t = -3.1379, df = 4, p-value = 0.03492

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

 -0.30157148 -0.01842852 

sample estimates:

mean of the differences 

                  -0.16 

che con un valore di p = 0.03492 indica significatività (anche se abbastanza risicata, di poco inferiore a 0.05) della differenza tra VEMS misurato prima e VEMS misurato dopo la somministrazione del farmaco.

Soggiacente al test t di Student è l'ipotesi che la media delle differenze t0 – t1 sia zero (ipotesi H₀ o ipotesi nulla). Ma il test eseguito con R verifica anche l'ipotesi alternativa, cioè che la media delle differenze non sia uguale a 0 (zero). 

L'intervallo di confidenza al 95% calcolato per questa ipotesi va da -0.30157148 a -0.01842852 e indica che la media delle differenze, uguale a -0.16, è significativamente diversa da 0 (non sarebbe significativamente diversa da 0 se l'intervallo di confidenza della media delle differenze includesse lo 0). Le due soluzioni si confermano l'una con l'altra.

Lo script prevede comunque di eseguire anche il test di Wilcoxon, questo è il risultato:

> wilcox.test(t0, t1, paired=TRUE, exact=FALSE) # test di Wilcoxon  per dati appaiati (test non parametrico)

        Wilcoxon signed rank test with continuity correction

data:  t0 and t1

V = 0, p-value = 0.09751

alternative hypothesis: true location shift is not equal to 0

Dato il ridotto numero di casi l'ipotesi alternativa non viene sviluppata. 

Da notare che assumendo come soglia di significatività il classico p = 0.05 il test di Wilcoxon con un valore di p = 0.09751 indica una differenza non significativa. Una cosa che non deve meravigliare, in quanto è noto che le statistiche non parametriche aumentano il valore di p ovvero rendono la differenza meno significativa rispetto a quella ottenuta con le statistiche parametriche (che hanno evidenziato una significatività risicata).

Il fatto che approcci statistici differenti possono portare a conclusioni opposte, e il fatto che una significatività statistica (risicata) non implica necessariamente una rilevanza pratica, meritano sempre attenzione nella fase di valutazione critica dei risultati di un'analisi statistica.

Una interessante rappresentazione grafica dei dati appaiati qui riportati può essere effettuata mediante un grafico prima-dopo (slopegraph) al quale si rimanda.

----------

[1] Un caso è per definizione un soggetto/oggetto univocamente identificato, come illustrato anche nel post Importazione dei dati da un file .csv.

[2] Campbell MJ, Machin D. Medical Statistics. A Commonsense Approach. John Wiley & Sons, New York, 1993, ISBN 0-471-93764-9, p. 142.

Test parametrici e non parametrici per due campioni indipendenti

Attenzione: nel caso di confronti multipli, cioè se dovete effettuare il confronto tra più campioni indipendenti come ad esempio, nel caso più semplice di tre campioni, confrontare il primo campione con il secondo campione, il primo campione con il terzo e il secondo campione con il terzo, dovete impiegare uno dei test riportati nel post Test parametrici e non parametrici per più campioni indipendenti. 

Qui ci occupiamo invece del confronto tra due campioni indipendenti [1] che viene effettuato verificando se le medie (la tradizionale misura di posizione, parametrica) o le mediane (la misura di posizione non parametrica alternativa alla media) dei due campioni differiscono significativamente tra loro. 

Come nel caso delle statistiche elementari anche in questo caso è necessario effettuare una analisi esplorativa dei dati e i test di normalità (gaussianità) [2] per decidere quale sia il test appropriato da impiegare:

→ il tradizionale test parametrico per il confronto tra medie, rappresentata dal test t di Student, se i dati sono distribuiti in modo gaussiano;

→ uno dei test non parametrici per il confronto tra mediane, che devono essere impiegati quando i dati non sono distribuiti in modo gaussiano, come il test di Wilcoxon per campioni indipendenti (in genere meglio noto come test U di Mann-Whitney) e il test di Kruskal-Wallis. 

Lo script è diviso in cinque parti al fine di sviluppare in modo analitico il percorso logico da seguire, che prevede:

→ verifica della distribuzione dei dati mediante i test di normalità (gaussianità);

→ se la distribuzione è gaussiana in entrambi i campioni messi a confronto, esecuzione del test F di Fisher per verificare se le varianze dei due campioni sono omogenee;

→ se la distribuzione è gaussiana e con il test di Fisher le varianze sono omogenee (non differiscono significativamente) esecuzione dei test t di Student per varianze omogenee;

→ se la distribuzione è gaussiana e con il test di Fisher le varianze non sono omogenee (differiscono significativamente) esecuzione dei test t di Student per varianze non omogenee;

→ se la distribuzione non è gaussiana, esecuzione dei test non parametrici (test di Wilcoxon per campioni indipendenti e/o test di Kruskall-Wallis).

Vediamo un esempio con i dati rilevati in atleti di sesso femminile e di sesso maschile e la domanda: altezza, peso, percentuale di grasso corporeo, emoglobina e gli altri dati rilevati, differiscono significativamente nei due sessi, o sono simili tra loro? I dati sono contenuti nella tabella ais del pacchetto DAAG, accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [3].

Copiate questa prima parte dello script, incollatela nella Console di R e premete ↵ Invio:

# CONFRONTO TRA DUE CAMPIONI INDIPENDENTI (1/5)

#

library(DAAG) # carica il pacchetto DAAG incluso il set di dati ais

str(ais) # mostra la struttura di ais

attach(ais) # consente di impiegare direttamente i nomi delle variabili

#

# test di normalità

#

library(nortest) # carica il pacchetto

lillie.test(wt[sex == "f"]) # test di Lilliefors (Kolmogorov-Smirnov) su peso corporeo donne

lillie.test(wt[sex == "m"]) # idem su peso corporeo uomini

lillie.test(pcBfat[sex == "f"]) # test di Lilliefors (Kolmogorov-Smirnov) su percentuale di grasso corporeo donne

lillie.test(pcBfat[sex == "m"]) # idem su percentuale di grasso corporeo uomini

#

Dopo avere caricato i dati con il pacchetto DAAG, avere mostrato la struttura di ais e avere consentito con attach(ais) di impiegare direttamente nel codice che segue i nomi delle variabili del set di dati ais, viene caricato il pacchetto nortest che consente di effettuare i test di normalità.

Scegliamo uno dei più classici e collaudati, il test di Lilliefors (Kolmogorov-Smirnov) e lo applichiamo al peso corporeo (variabile wt) e alla percentuale di grasso corporeo (variabile pcBfat) di uomini (m) e donne (f).

Nel caso del peso corporeo wt il test di Lilliefors conferma, sia per donne sia per uomini, che i dati non si discostano significativamente da una distribuzione gaussiana (p = 0.6532 e p = 0.6384 rispettivamente).

> lillie.test(wt[sex == "f"]) # test di Lilliefors (Kolmogorov-Smirnov) su peso corporeo donne

        Lilliefors (Kolmogorov-Smirnov) normality test

data:  wt[sex == "f"]

D = 0.05469, p-value = 0.6532

> lillie.test(wt[sex == "m"]) # idem su peso corporeo uomini

        Lilliefors (Kolmogorov-Smirnov) normality test

data:  wt[sex == "m"]

D = 0.054681, p-value = 0.6384

Nel caso della percentuale di grasso corporeo pcBfat invece ci troviamo di fronte a dati che se nelle donne non si discostano significativamente da una distribuzione gaussiana (p = 0.3291), negli uomini invece non sono distribuiti in modo gaussiano (p = 2.254e-08).

> lillie.test(pcBfat[sex == "f"]) # test di Lilliefors (Kolmogorov-Smirnov) su percentuale di grasso corporeo donne

        Lilliefors (Kolmogorov-Smirnov) normality test

data:  pcBfat[sex == "f"]

D = 0.066987, p-value = 0.3291

> lillie.test(pcBfat[sex == "m"]) # idem su percentuale di grasso corporeo uomini

        Lilliefors (Kolmogorov-Smirnov) normality test

data:  pcBfat[sex == "m"]

D = 0.17702, p-value = 2.254e-08

Per poter confrontare due campioni impiegando un test parametrico, è necessario che i dati siano distribuiti in modo gaussiano in entrambi, pertanto:

→ per il confronto del peso corporeo wt in uomini e donne possiamo impiegare un test parametrico, il test t di Student;

→ per il confronto della percentuale di grasso corporeo pcBfat in uomini e donne dobbiamo impiegare un test non parametrico (in realtà ne vediamo due, il test di Wilcoxon per campioni indipendenti e il test di Kruskal-Wallis).

Al peso corporeo, distribuito in modo gaussiano, applichiamo innanzitutto mediante la funzione var.test() il test F di Fisher per il controllo dell'omogeneità tra le varianze dei due campioni da mettere a confronto, che sono ricavati aggregando i valori della variabile peso corporeo (wt) in base al valore della variabile sesso (sex) mediante l'argomento wt~sex:

# test parametrico (confronto tra medie) per il peso corporeo, verifica preliminare (2/5)

#

var.test(wt~sex) # controllo dell'omogeneità delle varianze con il test F di Fisher

#

Il test F di Fisher indica che le varianze osservate in uomini e donne non sono significativamente diverse (la probabilità di osservare per caso una differenza tra le medie pari a quella effettivamente osservata è p = 0.2029) e anche l'ipotesi alternativa conferma questo:

> var.test(wt~sex) # controllo dell'omogeneità delle varianze con il test F di Fisher

        F test to compare two variances

data:  wt by sex

F = 0.77423, num df = 99, denom df = 101, p-value = 0.2029

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

 0.5221795 1.1488587

sample estimates:

ratio of variances 

         0.7742343 

Qualora le varianze dei due campioni fossero state significativamente diverse, sarebbe stato necessario procedere con il test t di Student per varianze non omogenee, impostando nella funzione t.test() l'argomento var.equal=FALSE.

Nel nostro caso invece procediamo con il test t di Student per varianze omogenee impostando nella funzione t.test() l'argomento var.equal=TRUE e aggregando di nuovo i valori della variabile peso corporeo (wt) in base ai valori della variabile sesso (sex) mediante l'argomento wt~sex.

# test parametrico (confronto tra medie) per il peso corporeo (3/5)

#

t.test(wt~sex, var.equal=TRUE) # test t di Student per varianze omogenee

#

Questi sono i risultati:

> t.test(wt~sex, var.equal=TRUE) # test t di Student per varianze omogenee

        Two Sample t-test

data:  wt by sex

t = -9.2272, df = 200, p-value < 2.2e-16

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

 -18.42589 -11.93717

sample estimates:

mean in group f mean in group m 

       67.34200        82.52353 

La media del peso corporeo nelle donne è 67.342, negli uomini è 82.52353 e le due medie sono significativamente diverse, in quanto la probabilità di osservare per caso una tale differenza è inferiore a 0.00000000000000022 (p < 2.2e-16) [1].

Soggiacente al test t di Student è l'ipotesi che le due medie siano uguali ovvero che la differenza tra le due medie sia zero (ipotesi H₀ o ipotesi nulla). Ma il test eseguito con R verifica anche l'ipotesi alternativa, cioè che la differenza tra le due medie non sia uguale a 0 (zero). L'intervallo di confidenza al 95% della differenza tra medie va da -18.42589 a -11.93717 e non include lo 0 indicando che la differenza tra le due medie, pari a -15.18153 (67.342 - 82.52353) è significativamente diversa da 0 (la differenza non sarebbe significativamente diversa se l'intervallo di confidenza della differenza tra le due medie includesse lo 0). Le due soluzioni si confermano l'una con l'altra e ci consentono di affermare che le donne pesano in media significativamente meno degli uomini.

Alla percentuale di grasso corporeo, che non è distribuita in modo gaussiano, applichiamo i test non parametrici per il confronto tra mediane con le funzioni wilcox.test() e kruskal.test(), in ciascuna delle quali i due campioni da mettere a confronto sono ottenuti aggregando i valori della variabile percentuale di grasso corporeo (pcBfat) in base ai valori della variabile sesso (sex) mediante l'argomento pcBfat~sex:

# test non parametrici (confronto tra mediane) per la percentuale di grasso corporeo (4/5)

#

wilcox.test(pcBfat~sex) # test di Wilcoxon per campioni indipendenti

#

kruskal.test(pcBfat~sex) # test di Kruskal-Wallis

#

median(pcBfat[sex == "f"]) # mediana della percentuale di grasso corporeo per sesso = f

median(pcBfat[sex == "m"]) # mediana della percentuale di grasso corporeo per sesso = m

#

Diversamente dalla funzione t.test(), che riporta al termine dei calcoli il valore della media dei due campioni, le funzioni wilcox.test() e kruskal.test() non forniscono la mediana, per cui allo script sono state aggiunte due righe di codice per calcolarla sia nelle donne sia negli uomini.

La mediana della percentuale di grasso corporeo nelle donne è 17.94, negli uomini è 8.625 ed entrambi i test confermano che le due mediane sono significativamente diverse, la probabilità di osservare per caso una tale differenza è per entrambi i test inferiore a 0.00000000000000022 (p < 2.2e-16) [4].

> wilcox.test(pcBfat~sex) # test di Wilcoxon per campioni indipendenti

        Wilcoxon rank sum test with continuity correction

data:  pcBfat by sex

W = 9417.5, p-value < 2.2e-16

alternative hypothesis: true location shift is not equal to 0

> kruskal.test(pcBfat~sex) # test di Kruskal-Wallis

        Kruskal-Wallis rank sum test

data:  pcBfat by sex

Kruskal-Wallis chi-squared = 108.03, df = 1, p-value < 2.2e-16

> median(pcBfat[sex == "f"]) # mediana della percentuale di grasso corporeo per sesso = f

[1] 17.94

> median(pcBfat[sex == "m"]) # mediana della percentuale di grasso corporeo per sesso = m

[1] 8.625

Quindi possiamo affermare che nelle donne la mediana della percentuale di grasso corporeo è significativamente maggiore di quella rilevata uomini.

Un interessante integrazione ai risultati numerici ci viene ancora una volta dalla grafica, e in particolare dai grafici a scatola con i baffi:

# analisi grafica non parametrica delle distribuzioni (5/5)

#

windows() # apre una nuova finestra

par(mfrow=c(1,2)) # suddivide la finestra in due quadranti, uno per grafico

#

boxplot(wt~sex, data=ais, xlab="Sesso", ylab="Peso corporeo in kg", notch=TRUE, col="green") # traccia i boxplot per sesso del peso corporeo

#

boxplot(pcBfat~sex, data=ais, xlab="Sesso", ylab="Grasso corporeo in %", notch=TRUE, col="green") # traccia i boxplot per sesso della percentuale di grasso corporeo

#

Questo è il grafico

nel quale i boxplot per sesso del peso corporeo e della percentuale di grasso corporeo sono tracciati con una incisura (notch=TRUE) che rappresenta i limiti di confidenza al 95% della mediana. Se i limiti di confidenza delle mediane confrontate non si sovrappongono, le mediane sono significativamente diverse.

Qui il fatto interessante è che andiamo ben oltre l'impiego della grafica come complemento dei risultati dell'analisi numerica: abbiamo un esempio di come sia possibile realizzare un test statistico (non parametrico) mediante una rappresentazione grafica dei dati. La conclusione è che sia per il peso corporeo sia per la percentuale di grasso corporeo le incisure dei boxplot di donne e uomini non si sovrappongono, e pertanto le corrispondenti mediane risultano significativamente diverse in entrambi i casi.

----------

[1] Se i due campioni non sono indipendenti impiegare uno dei test riportati nel post Test parametrici e non parametrici per dati appaiati.

[2] Vedere il post Analisi esplorativa dei dati e i post a questo collegati.

[3] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG. 

[4] Viene riportato questo valore quando R non è più in grado di approssimare numericamente il risultato.