Test parametrici e non parametrici per più campioni indipendenti

Per confrontare due campioni indipendenti abbiamo visto [1] che è necessario impiegare:

→ metodi parametrici quando i dati sono distribuiti in modo gaussiano;

→ metodi non parametrici quando i dati non sono distribuiti in modo gaussiano.

Lo stesso vale per il confronto tra più di due campioni indipendenti. Ma in questo caso però vale anche un altro principio: non è corretto applicare un test per due campioni indipendenti ripetitivamente a tutte le diverse possibili coppie di campioni. Così, ad esempio, nel caso più semplice, quello di tre campioni, non è corretto confrontare il primo campione con il secondo campione, il primo campione con il terzo e il secondo campione con il terzo. Questo perché la probabilità di commettere un errore di primo tipo [2], cioè la probabilità di considerare la differenza tra le medie di due campioni significativa quando invece non è significativa (probabilità di un falso positivo), nel caso di confronti multipli viene moltiplicata per il numero di confronti effettuati.

In altre parole un valore che nel caso di due campioni (per i quali un solo confronto è possibile) corrisponde a una probabilità di commettere un errore di primo tipo del 5%:
→ nel caso di tre campioni (essendo 3 il numero di confronti possibili) corrisponde a una probabilità di commettere un errore di primo tipo del 3 x 5% = 15%;
→ nel caso di quattro campioni (essendo 6 il numero di confronti possibili) corrisponde a una probabilità di commettere un errore di primo tipo del 6 x 5% = 30%;
→ nel caso di cinque campioni (essendo 10 il numero di confronti possibili) corrisponde a una probabilità di commettere un errore di primo tipo del 10 x 5% = 50%;
e così via.

Reciprocamente se nel confronto tra due campioni si assume come valore soglia per la significatività il classico p = 0.05 è indispensabile ricordarsi che questo valore:
→ quando si confrontano tre campioni deve essere sostituito con p = 0.017 (cioè 0.05/3);
→ quando si confrontano quattro campioni deve essere sostituito con p = 0.008 (cioè 0.05/6);
→ quando si confrontano cinque campioni deve essere sostituito con p = 0.005 (cioè 0.05/10);
e così via.

Anche se si potrebbe procedere seguendo queste semplici regole, di fatto sono disponibili test che, nel caso di confronti multipli, applicano automaticamente le correzioni necessarie per controbilanciare l'aumento dei falsi positivi, cioè l'aumento della probabilità di considerare la differenza tra le medie di due campioni significativa quando invece non è significativa.

Vediamo quindi i test parametrici e non parametrici per più campioni indipendenti impiegando come esempio il confronto tra i valori di concentrazione dell'emoglobina nel sangue rilevata in 202 atleti australiani suddivisi per sport e riportata nella variabile hg della tabella ais inclusa nel pacchetto DAAG. Accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [3].

Test parametrici

Copiate lo script che segue, incollatelo nella Console di R e premete ↵ Invio:

# CONFRONTARE TRA DI LORO PIU' CAMPIONI con un metodo parametrico

# test t di Student con la correzione di Bonferroni per confronti multipli

#

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

#

boxplot(hg~sport, data=ais, horizontal=FALSE, notch=FALSE, col="green", las=2, xlab="Sport praticato", ylab="Emoglobina nel sangue (g/dL)") # boxplot della concentrazione di emoglobina per sport praticato

#

pairwise.t.test(ais$hg, ais$sport, p.adjust.method="bonferroni") # test t di Student per campioni indipendenti con la correzione di Bonferroni per confronti multipli

#

Per prima cosa visualizziamo i dati realizzando un grafico a scatola con i baffi della concentrazione hg dell'emoglobina nel sangue (espressa in g/dL) aggregando i valori in base allo sport praticato (hg~sport).

Ora per avere un riscontro quantitativo e oggettivo delle differenze tra le medie riportate nel grafico utilizziamo il più classico test parametrico, il test t di Student (un altro test parametrico lo vedremo con lo script successivo). 

Lo facciamo con la funzione pairwise.t.test() che consente di aggiustare i valori di p (p.adjust.method) con la correzione per confronti multipli di Bonferroni (method=”bonferroni”):

       Pairwise comparisons using t tests with pooled SD 

data:  ais$hg and ais$sport 

        B_Ball Field   Gym    Netball Row    Swim   T_400m T_Sprnt Tennis

Field   0.0023 -       -      -       -      -      -      -       -     

Gym     1.0000 0.1012  -      -       -      -      -      -       -     

Netball 0.0050 2.4e-11 1.0000 -       -      -      -      -       -     

Row     1.0000 0.1676  1.0000 6.5e-07 -      -      -      -       -     

Swim    1.0000 1.0000  1.0000 2.3e-06 1.0000 -      -      -       -     

T_400m  1.0000 0.8321  1.0000 2.7e-07 1.0000 1.0000 -      -       -     

T_Sprnt 0.0015 1.0000  0.0606 5.1e-11 0.0928 0.5925 0.4455 -       -     

Tennis  1.0000 0.2188  1.0000 0.0178  1.0000 1.0000 1.0000 0.1213  -     

W_Polo  0.0036 1.0000  0.1078 8.8e-11 0.2180 1.0000 0.9715 1.0000  0.2490

P value adjustment method: bonferroni 

Il test t di Student con la correzione di Bonferroni conferma che le atlete (sono solo donne) che praticano Netball hanno in comune una bassa concentrazione dell'emoglobina, e che la concentrazione media di questo gruppo si discosta significativamente dalle concentrazioni medie dei rimanenti gruppi di atleti, con la sola eccezione (p = 1.0000) del gruppo delle atlete che pratica la ginnastica (Gym). Questo è il dato certamente più rilevante, con poche altre differenze significative che sono evidenziate nella tabella.

Da notare che la funzione pairwise.t.test() prevede le seguenti correzioni alternative meno conservative [4] di quella di Bonferroni: secondo Holm (1979) (method="holm"), Hochberg (1988) (method="hochberg"), Hommel (1988) (method="hommel"), Benjamini & Hochberg (1995) (method="BH" o method="fdr"), Benjamini & Yekutieli (2001) (method="BY"), e method="none" per nessuna correzione (ovviamente sconsigliato).

Test non parametrici

Vediamo ora di applicare ai dati alcuni test non parametrici, contenuti nel pacchetto PMCMRplus [5], che è necessario scaricare e installare dal CRAN. Copiate lo script che segue, incollatelo nella Console di R e premete ↵ Invio:

# CONFRONTARE TRA DI LORO PIU' CAMPIONI con metodi non parametrici

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

#

boxplot(hg~sport, data=ais, horizontal=FALSE, notch=FALSE, col="green", las=2, xlab="Sport praticato", ylab="Emoglobina nel sangue (g/dL)") # boxplot della concentrazione di emoglobina per sport praticato

#

library(PMCMRplus) # carica il pacchetto

spearmanTest(hg~sport, data=ais) # rho di Spearman (test omnibus, non parametrico)

kwAllPairsConoverTest(hg~sport, data=ais) # test di Conover (non parametrico)

kwAllPairsDunnTest(hg~sport, data=ais) # test di Dunn (non parametrico)

kwAllPairsNemenyiTest(hg~sport, data=ais) # test di Nemenyi (non parametrico)

#

Quelli che qui impieghiamo sono tutti test non parametrici che quindi non prevedono che i dati siano distribuiti in modo gaussiano e sono eseguiti mediante:

→ la funzione spearmanTest() per il test ρ (rho) di Spearman [6],  un test non parametrico che riporta un risultato unico per tutte le coppie di campioni (test omnibus);

→ la funzione kwAllPairsConoverTest() per il test di Conover, un test non parametrico basato sui ranghi di tipo Kruskal, che riporta un risultato per ogni coppia di campioni;

→ la funzione kwAllPairsDunnTest() per il test di Dunn, un test non parametrico basato sui ranghi di tipo Kruskal, che riporta un risultato per ogni coppia di campioni;

→ la funzione kwAllPairsNemenyiTest() per il test di Nemenyi, un test non parametrico basato sul confronto tra i ranghi, che riporta un risultato per ogni coppia di campioni.

I test di Conover, di Dunn e di Nemenyi, che riportano un risultato per ogni coppia di campioni, forniscono risultati sostanzialmente sovrapponibili, qui riporto per semplicità solamente i risultati dell'ultimo di questi.

> kwAllPairsNemenyiTest(hg~sport, data=ais) # test di Nemenyi (non parametrico)

        Pairwise comparisons using Tukey-Kramer-Nemenyi all-pairs test with Tukey-Dist approximation

data: hg by sport

        B_Ball Field   Gym   Netball Row   Swim  T_400m T_Sprnt Tennis

Field   0.023  -       -     -       -     -     -      -       -     

Gym     0.997  0.152   -     -       -     -     -      -       -     

Netball 0.028  4.3e-09 0.996 -       -     -     -      -       -     

Row     0.981  0.201   0.905 8.6e-05 -     -     -      -       -     

Swim    0.867  0.697   0.785 7.8e-05 1.000 -     -      -       -     

T_400m  0.844  0.578   0.785 2.0e-05 1.000 1.000 -      -       -     

T_Sprnt 0.133  1.000   0.279 5.9e-07 0.573 0.934 0.891  -       -     

Tennis  1.000  0.196   0.997 0.158   0.999 0.981 0.981  0.441   -     

W_Polo  0.018  1.000   0.126 5.0e-09 0.160 0.614 0.494  1.000   0.158 

P value adjustment method: single-step

alternative hypothesis: two.sided

Messaggio di avvertimento:

In kwAllPairsNemenyiTest.default(c(12.3, 12.7, 11.6, 12.6, 14, 12.5,  :

  Ties are present, p-values are not corrected.

Il pacchetto PMCMRplus include tra le altre anche la funzione tukeyTest() per il test di Tukey, un test parametrico per il confronto di campioni con distribuzione gaussiana e con uguale varianza, che riporta un risultato per ogni coppia di campioni, un test sovrapponibile per applicazione e conclusioni al test t di Student con la correzione di Bonferroni, anche se meno conservativo. La sintassi è sempre la stessa, per visualizzarne i risultati copiate e incollate nella Console di R questa riga di codice e premete ↵ Invio:

tukeyTest(hg~sport, data=ais) # test di Tukey (parametrico)

Da notare che i due test parametrici, il test di Tukey e il test t di Student con la correzione di Bonferroni, forniscono risultati molto simili:

        Pairwise comparisons using Tukey's test

data: hg by sport

        B_Ball Field   Gym    Netball Row    Swim   T_400m T_Sprnt Tennis

Field   0.0020 -       -      -       -      -      -      -       -     

Gym     0.9982 0.0671  -      -       -      -      -      -       -     

Netball 0.0043 2.4e-11 0.9553 -       -      -      -      -       -     

Row     0.8104 0.1028  0.8151 6.4e-07 -      -      -      -       -     

Swim    0.6798 0.4141  0.7173 2.3e-06 1.0000 -      -      -       -     

T_400m  0.5633 0.3463  0.6829 2.7e-07 1.0000 1.0000 -      -       -     

T_Sprnt 0.0013 1.0000  0.0429 5.1e-11 0.0623 0.2739 0.2224 -       -     

Tennis  1.0000 0.1278  0.9871 0.0142  0.9993 0.9920 0.9872 0.0784  -     

W_Polo  0.0031 1.0000  0.0709 8.8e-11 0.1275 0.4467 0.3833 1.0000  0.1418

P value adjustment method: single-step

alternative hypothesis: two.sided

In conclusione: i risultati di due test parametrici per il confronto tra medie (il test t di Student con la correzione di Bonferroni e il test di Tukey) sono simili a quelli dei principali test non parametrici per il confronto tra mediane (test ρ di Spearman, test di Conover, di Dunn, di Nemenyi). Questo è determinato dal fatto che in questo caso le distribuzioni dei dati non si allontanano molto da una distribuzione gaussiana: e per definizione in una distribuzione gaussiana i risultati dei test parametrici e dei test non parametrici coincidono, e coincidono medie e mediane. Tuttavia va sottolineato che un approccio rigoroso e puntuale deve prevedere, al fine di consentire la corretta applicazione dei test, la verifica preliminare della gaussianità dei dati – che qui per semplicità non è stata riportata – con uno dei metodi che trovate descritti [7]. 

----------

[1] Vedere il post Test parametrici e non parametrici per due campioni indipendenti.

[2] Wikipedia. Type I and type II errors. URL consultato il 10/01/2020: http://bit.ly/2R1VdyF

[3] Vedere il post Il set di dati ais. La concentrazione dell'emoglobina viene espressa in grammi per decilitro di sangue (g/dL). Nel post trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG. 

[4] Un metodo statistico “meno conservativo” a parità di condizioni riporta più frequentemente differenze significative. Personalmente preferisco in ogni caso impiegare i metodi più conservativi, in quanto riducono la probabilità di considerare significativa una differenza che invece non è significativa (ovvero riducono la probabilità di un falso positivo).

[5] PMCMR è l'acronimo di Calculate Pairwise Multiple Comparisons of Mean Rank Sums Extended. Si tratta di un pacchetto che include una serie molto ampia sia di test omnibus (che forniscono un solo risultato per tutti i confronti effettuati), sia per il confronto tra tutte le coppie di campioni, sia parametrici, sia non parametrici. Vedere sul CRAN il reference manual PMCMRplus: Calculate Pairwise Multiple Comparisons of Mean Rank Sums Extended.

https://cran.r-project.org/web/packages/PMCMRplus/index.html

[6] Per il test ρ (rho) di Spearman vedere anche il post Coefficienti di correlazione parametrici e non parametrici.

[7] Vedere ad esempio:

→ il post Test di normalità (gaussianità)

→ il post Tabulare una serie di test di normalità (gaussianità)

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.