Calcolare la media mobile e la mediana mobile

Data una sequenza di valori numerici la media mobile è la media calcolata in corrispondenza di ciascun dato su un numero fisso di dati, determinato da una finestra di ampiezza prestabilita che scorre i dati e nella quale il dato più vecchio viene ogni volta sostituito con un nuovo dato.

Così ad esempio per una finestra con l'ampiezza di 5 dati, in corrispondenza dei primi 4 dati la media mobile non può essere calcolata, in corrispondenza del dato numero 5 viene calcolata la media dei dati da 1 a 5, in corrispondenza del dato numero 6 viene calcolata la media dei dati da 2 a 6, in corrispondenza del dato numero 7 viene calcolata la media dei dati da 3 a 7 e così via. La mediana mobile a sua volta è la mediana calcolata in corrispondenza di ciascun dato nello stesso modo.

Per calcolare la media mobile e la mediana mobile facciamo ricorso al pacchetto zoo, che deve essere scaricato dal CRAN.

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

# CALCOLARE LA MEDIA MOBILE E LA MEDIANA MOBILE 

#

library(zoo) # carica il pacchetto

rileva_x <- seq(1:50) # numero progressivo della rilevazione in ascisse

valore_y <- c(5, 9, 12, 7, 8, 13, 26, 7, 8, 19, 10, 7, 5, 36, 12, 18, 31, 17, 16, 43, 7, 12, 9, 27, 15, 18, 24, 29, 32, 23, 29, 12, 7, 48, 34, 45, 19, 12, 16, 18, 22, 31, 38, 7, 11, 4, 38, 36 ,43, 15) # valore osservato corrispondente in ordinate

mydata <- data.frame(rileva_x, valore_y) # combina i vettori in una matrice

#

mydata$mediamob <- rollmean(valore_y, k=7, fill=NA, align='right') # calcola la media mobile (dato corrente + 6 precedenti) e la aggiunge in una nuova colonna

mydata$medianamob <- rollmedian(valore_y, k=7, fill=NA, align='right') # calcola la mediana mobile (dato corrente + 6 precedenti) e la aggiunge in una nuova colonna

mydata # mostra la matrice con i dati definitivi

#

plot(mydata$rileva_x, mydata$valore_y, xlim=c(0,50), ylim=c(0,60), type="l", lty=1, lwd=1, col="black", main="Valori osservati, media mobile e mediana mobile", xlab="Rilevazioni", ylab="Valori osservati") # grafico dei dati originali

lines(mydata$rileva_x, mydata$mediamob, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="red") # sovrappone il grafico della media mobile

lines(mydata$rileva_x, mydata$medianamob, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="blue") # sovrappone il grafico della mediana mobile

#

legend(0,60, legend=c("Valori osservati nelle rilevazioni", "Media mobile dei valori osservati", "Mediana mobile dei valori osservati"), col=c("black", "red", "blue"), lty=c(1,2,2), lwd=c(1,1,1), cex=0.8) # aggiunge la legenda

#

Dopo avere caricato il pacchetto zoo i dati impiegati come esempio sono generati in questo modo:

→ nel vettore rileva_x viene riportato il numero progressivo della rilevazione;

→ nel vettore valore_y viene riportato il valore numerico osservato in corrispondenza di ciascuna rilevazione, valore numerico sul quale calcoleremo la media mobile e la mediana mobile;

→ con la funzione dataframe() i due vettori sono combinati nella matrice mydata [1].

A questo punto alla matrice mydata sono aggiunte due colonne:

→ la colonna/vettore mydata$mediamob nella quale è riportata la media mobile calcolata in corrispondenza di ogni dato;

→ la colonna/vettore mydata$medianamob nella quale è riportata la mediana mobile calcolata in corrispondenza di ogni dato.

La media mobile e la mediana mobile sono calcolate rispettivamente con le funzioni rollmean() e rollmedian() che impiegano i seguenti argomenti:

→ la variabile (valore_y) sulla quale calcolare la media mobile;

→ l'ampiezza della finestra (k=7) che specifica che la media (o mediana) mobile va calcolata, in corrispondenza di ogni dato, su quel dato e sui sei dati che lo precedono;

→ fill=NA che specifica di riportare NA in corrispondenza dei valori non computabili;

→ align='right' che - in quanto con k=7 la media e la mediana possono essere calcolate solamente a partire dal settimo dato - specifica di riportare i valori non computabili NA nelle prime sei posizioni.

Per i dettagli delle funzioni impiegate digitare help(nomedellafunzione) nella Console di R. Al bisogno potete consultare il manuale del pacchetto zoo [2].

Non resta quindi che mostrare il contenuto della matrice mydata nella versione definitiva contenente il numero della rilevazione (rileva_x), il valore osservato (valore_y) e le rispettive media mobilea (mediamob) e mediana mobile (medianamob)

> mydata # mostra la matrice con i dati definitivi

   rileva_x valore_y mediamob medianamob

1         1        5       NA         NA

2         2        9       NA         NA

3         3       12       NA         NA

4         4        7       NA         NA

5         5        8       NA         NA

6         6       13       NA         NA

7         7       26 11.42857          9

8         8        7 11.71429          9

9         9        8 11.57143          8

10       10       19 12.57143          8

11       11       10 13.00000         10

12       12        7 12.85714         10

13       13        5 11.71429          8

14       14       36 13.14286          8

15       15       12 13.85714         10

16       16       18 15.28571         12

17       17       31 17.00000         12

18       18       17 18.00000         17

19       19       16 19.28571         17

20       20       43 24.71429         18

21       21        7 20.57143         17

22       22       12 20.57143         17

23       23        9 19.28571         16

24       24       27 18.71429         16

25       25       15 18.42857         15

26       26       18 18.71429         15

27       27       24 16.00000         15

28       28       29 19.14286         18

29       29       32 22.00000         24

30       30       23 24.00000         24

31       31       29 24.28571         24

32       32       12 23.85714         24

33       33        7 22.28571         24

34       34       48 25.71429         29

35       35       34 26.42857         29

36       36       45 28.28571         29

37       37       19 27.71429         29

38       38       12 25.28571         19

39       39       16 25.85714         19

40       40       18 27.42857         19

41       41       22 23.71429         19

42       42       31 23.28571         19

43       43       38 22.28571         19

44       44        7 20.57143         18

45       45       11 20.42857         18

46       46        4 18.71429         18

47       47       38 21.57143         22

48       48       36 23.57143         31

49       49       43 25.28571         36

50       50       15 22.00000         15

che sono impiegati per la successiva rappresentazione grafica

che viene realizzata con tre sole funzioni:

→ la funzione plot() che traccia il grafico del valore osservato mydata$valore_y in corrispondenza di ciascuna rilevazione mydata$rileva_x impiegando una linea (type="l") continua (lty=1) di larghezza unitaria (lwd=1) e di colore nero (col="black");

→ la funzione lines() che viene impiegata una prima volta per sovrapporre il grafico della media mobile con una linea tratteggiata (lty=2) di colore rosso (col="red") e una seconda volta per sovrapporre il grafico della mediana mobile con una linea tratteggiata (lty=2) di colore blu (col="blue") 

→ la funzione legend() che riporta la legenda esplicativa dei grafici [3].

Ovviamente lo script può essere riadattato per rappresentare dati differenti ma anche per riportare solamente l'uno o l'altro dei grafici.

----------

[1] Parliamo di array o vettore nel caso di dati numerici monodimensionali, disposti su una sola riga, 

8

6

11

7

di matrice nel caso di dati numerici disposti su più righe e più colonne

8	9	15	14
6	7	18	12
11	8	17	13
7	4	19	17

e di tabella nei casi in cui il contenuto, disposto su più righe e più colonne, è rappresentato oltre che da dati numerici, anche da testo e/o operatori logici

M	7	9	VERO
F	3	12	VERO
F	5	10	FALSO

[2] Vedere la voce Reference manual del pacchetto zoo: S3 Infrastructure for Regular and Irregular Time Series (Z's Ordered Observations).

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

[3] Per un approfondimento sull'impiego della funzione legend() rimando al post Aggiungere la legenda a un grafico.

Calcolare la media cumulativa e la mediana cumulativa

Data una sequenza di valori numerici la media cumulativa è la media, calcolata in corrispondenza di ciascun dato, di tutti i dati dal primo fino a quello corrente.

Così ad esempio in corrispondenza del dato numero 5 la media cumulativa è la media dei dati da 1 a 5, in corrispondenza del dato numero 6 la media cumulativa è la media dei dati da 1 a 6, in corrispondenza del dato numero 7 la media cumulativa è la media dei dati da 1 a 7 e così via. Lo stesso principio vale per la mediana cumulativa che a sua volta è la mediana, calcolata in corrispondenza di ciascun dato, di tutti i dati dal primo fino a quello corrente.

Per il calcolo della media cumulativa e della mediana cumulativa impieghiamo le funzioni contenute nel pacchetto cumstats, che deve essere scaricato dal CRAN.

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

# CALCOLARE LA MEDIA CUMULATIVA E LA MEDIANA CUMULATIVA 

#

library(cumstats) # carica il pacchetto

rileva_x <- seq(1:50) # numero progressivo della rilevazione in ascisse

valore_y <- c(5, 9, 12, 7, 8, 13, 26, 7, 8, 19, 10, 7, 5, 36, 12, 18, 31, 17, 16, 43, 7, 12, 9, 27, 15, 18, 24, 29, 32, 23, 29, 12, 7, 48, 34, 45, 19, 12, 16, 18, 22, 31, 38, 7, 11, 4, 38, 36 ,43, 15) # valore osservato corrispondente in ordinate

mydata <- data.frame(rileva_x, valore_y) # combina i vettori in una matrice

#

mydata$mediacum <- cummean(mydata$valore_y) # calcola la media cumulativa e la aggiunge in una nuova colonna

mydata$medianacum <- cummedian(mydata$valore_y) # calcola la mediana cumulativa e la aggiunge in una nuova colonna

mydata # mostra la matrice con i dati definitivi

#

plot(mydata$rileva_x, mydata$valore_y, xlim=c(0,50), ylim=c(0,60), type="l", lty=1, lwd=1, col="black", main="Valori osservati, media cumulativa e mediana cumulativa", xlab="Rilevazioni", ylab="Valori osservati") # grafico dei dati originali

lines(mydata$rileva_x, mydata$mediacum, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="red") # sovrappone grafico della media cumulativa

lines(mydata$rileva_x, mydata$medianacum, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="blue") # sovrappone grafico della mediana cumulativa

#

legend(0,60, legend=c("Valori osservati nelle rilevazioni", "Media cumulativa dei valori osservati", "Mediana cumulativa dei valori osservati"), col=c("black", "red", "blue"), lty=c(1,2,2), lwd=c(1,1,1), cex=0.8) # aggiunge la legenda

#

Dopo avere caricato il pacchetto cumstats i dati impiegati come esempio sono generati in questo modo:

→ nel vettore rileva_x viene riportato il numero progressivo della rilevazione;

→ nel vettore valore_y viene riportato il valore numerico osservato in corrispondenza di ciascuna rilevazione, valore numerico sul quale calcoleremo la media cumulativa e la mediana cumulativa.

→ con la funzione dataframe() i due vettori sono combinati nella matrice mydata [1].

A questo punto alla matrice mydata sono aggiunte due colonne:

→ la colonna (vettore) mediacum nella quale è riportata la media cumulativa calcolata in corrispondenza di ogni dato con la funzione cummmean();

→ la colonna (vettore) medianacum nella quale è riportata la mediana cumulativa calcolata in corrispondenza di ogni dato con la funzione cummedian().

Per i dettagli delle funzioni impiegate digitare help(nomedellafunzione) nella Console di R. Al bisogno potete consultare il manuale del pacchetto cumstats [2].

Non resta quindi che mostrare il contenuto della matrice mydata nella versione definitiva contenente il numero della rilevazione (rileva_x), il valore osservato (valore_y) e le rispettive media cumulativa (mediacum) e mediana cumulativa (medianacum)

> mydata # mostra la matrice con i dati definitivi

   rileva_x valore_y  mediacum medianacum

1         1        5  5.000000        5.0

2         2        9  7.000000        7.0

3         3       12  8.666667        9.0

4         4        7  8.250000        8.0

5         5        8  8.200000        8.0

6         6       13  9.000000        8.5

7         7       26 11.428571        9.0

8         8        7 10.875000        8.5

9         9        8 10.555556        8.0

10       10       19 11.400000        8.5

11       11       10 11.272727        9.0

12       12        7 10.916667        8.5

13       13        5 10.461538        8.0

14       14       36 12.285714        8.5

15       15       12 12.266667        9.0

16       16       18 12.625000        9.5

17       17       31 13.705882       10.0

18       18       17 13.888889       11.0

19       19       16 14.000000       12.0

20       20       43 15.450000       12.0

21       21        7 15.047619       12.0

22       22       12 14.909091       12.0

23       23        9 14.652174       12.0

24       24       27 15.166667       12.0

25       25       15 15.160000       12.0

26       26       18 15.269231       12.0

27       27       24 15.592593       12.0

28       28       29 16.071429       12.5

29       29       32 16.620690       13.0

30       30       23 16.833333       14.0

31       31       29 17.225806       15.0

32       32       12 17.062500       14.0

33       33        7 16.757576       13.0

34       34       48 17.676471       14.0

35       35       34 18.142857       15.0

36       36       45 18.888889       15.5

37       37       19 18.891892       16.0

38       38       12 18.710526       15.5

39       39       16 18.641026       16.0

40       40       18 18.625000       16.0

41       41       22 18.707317       16.0

42       42       31 19.000000       16.5

43       43       38 19.441860       17.0

44       44        7 19.159091       16.5

45       45       11 18.977778       16.0

46       46        4 18.652174       16.0

47       47       38 19.063830       16.0

48       48       36 19.416667       16.5

49       49       43 19.897959       17.0

50       50       15 19.800000       16.5

che sono impiegati per la successiva rappresentazione grafica

che viene realizzata con tre sole funzioni:

→ la funzione plot() che traccia il grafico del valore osservato mydata$valore_y in corrispondenza di ciascuna rilevazione mydata$rileva_x impiegando una linea (type="l") continua (lty=1) di larghezza unitaria (lwd=1) e di colore nero (col="black");

→ la funzione lines() che viene impiegata una prima volta per sovrapporre il grafico della media cumulativa con una linea tratteggiata (lty=2) di colore rosso (col="red") e una seconda volta per sovrapporre il grafico della mediana cumulativa con una linea tratteggiata (lty=2) di colore blu (col="blue") 

→ la funzione legend() che riporta la legenda esplocativa dei grafici [3].

Ovviamente lo script può essere riadattato per rappresentare dati differenti ma anche per riportare solamente l'uno o l'altro dei grafici.

----------

[1] Parliamo di array o vettore nel caso di dati numerici monodimensionali, disposti su una sola riga, 

8

6

11

7

di matrice nel caso di dati numerici disposti su più righe e più colonne

8	9	15	14
6	7	18	12
11	8	17	13
7	4	19	17

e di tabella nei casi in cui il contenuto, disposto su più righe e più colonne, è rappresentato oltre che da dati numerici, anche da testo e/o operatori logici

M	7	9	VERO
F	3	12	VERO
F	5	10	FALSO

[2] Vedere la voce Reference manual del pacchetto cumstats: Cumulative Descriptive Statistics.

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

[3] Per un approfondimento sull'impiego della funzione legend() rimando al post Aggiungere la legenda a un grafico.

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.