Estrazione dei dati per fattore

Sul set di dati ais contenuto nel pacchetto DAAG – accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [1] – vogliamo eseguire una analisi statistica e grafica selezionando alcuni dati specifici: vediamo cosa contiene quindi lo impieghiamo in un esempio di estrazione dei dati per fattore.

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

# ESTRAZIONE DEI DATI PER FATTORE

#

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

head(ais, 3) # mostra le prime tre righe del set di dati ais

tail(ais,3) # mostra le ultime tre righe del set di dati ais

str(ais) # mostra la struttura dei dati 

#

Dopo avere caricato il pacchetto, con head(ais, 3) vediamo le prime tre righe di ais

> head(ais, 3) # mostra le prime tre righe del set di dati ais

   rcc wcc   hc   hg ferr   bmi   ssf pcBfat   lbm    ht   wt sex  sport

1 3.96 7.5 37.5 12.3   60 20.56 109.1  19.75 63.32 195.9 78.9   f B_Ball

2 4.41 8.3 38.2 12.7   68 20.67 102.8  21.30 58.55 189.7 74.4   f B_Ball

3 4.14 5.0 36.4 11.6   21 21.86 104.6  19.88 55.36 177.8 69.1   f B_Ball

con tail(ais, 3) vediamo le ultime tre righe 

> tail(ais,3) # mostra le ultime tre righe del set di dati ais

     rcc wcc   hc   hg ferr   bmi  ssf pcBfat lbm    ht   wt sex  sport

200 5.03 6.4 42.7 14.3  122 22.01 47.6   8.51  68 183.1 73.8   m Tennis

201 4.97 8.8 43.0 14.9  233 22.34 60.4  11.50  63 178.4 71.1   m Tennis

202 5.38 6.3 46.0 15.7   32 21.07 34.9   6.26  72 190.8 76.7   m Tennis

e con str(ais) vediamo la struttura dei dati contenuti in ais

> str(ais) # mostra la struttura dei dati 

'data.frame':   202 obs. of  13 variables:

 $ rcc   : num  3.96 4.41 4.14 4.11 4.45 4.1 4.31 4.42 4.3 4.51 ...

 $ wcc   : num  7.5 8.3 5 5.3 6.8 4.4 5.3 5.7 8.9 4.4 ...

 $ hc    : num  37.5 38.2 36.4 37.3 41.5 37.4 39.6 39.9 41.1 41.6 ...

 $ hg    : num  12.3 12.7 11.6 12.6 14 12.5 12.8 13.2 13.5 12.7 ...

 $ ferr  : num  60 68 21 69 29 42 73 44 41 44 ...

 $ bmi   : num  20.6 20.7 21.9 21.9 19 ...

 $ ssf   : num  109.1 102.8 104.6 126.4 80.3 ...

 $ pcBfat: num  19.8 21.3 19.9 23.7 17.6 ...

 $ lbm   : num  63.3 58.5 55.4 57.2 53.2 ...

 $ ht    : num  196 190 178 185 185 ...

 $ wt    : num  78.9 74.4 69.1 74.9 64.6 63.7 75.2 62.3 66.5 62.9 ...

 $ sex   : Factor w/ 2 levels "f","m": 1 1 1 1 1 1 1 1 1 1 ...

 $ sport : Factor w/ 10 levels "B_Ball","Field",..: 1 1 1 1 1 1 1 1 1 1 ...

Tutto molto chiaro e molto semplice, abbiamo 202 osservazioni (202 obs) di 13 variabili di cui: 

→ 11 variabili numeriche (num), risultato di misure di dati ematologici e di dati biometrici effettuate in altrettanti atleti australiani [1]; 

→ due variabili "fattore" (Factor) cioè variabili qualitative che sono state impiegate per classificare gli atleti in base al sesso (sex) in due classi (2 levels) e in base allo sport praticato (sport) in dieci classi (10 levels). 

Eccoci quindi al punto: vogliamo verificare in via preliminare se, nel caso della variabile emoglobina (hg), all'interno di ciascuno dei dieci sport praticati (sport) i dati sono distribuiti in modo normale (gaussiano), impiegando il test di Shapiro-Wilk. 

Per farlo potremmo innanzitutto estrarre i valori di concentrazione dell'emoglobina (ais$hg) con la funzione which() per il primo sport (ais$sport=='B_Ball') nell'elenco, impiegando questo codice

#

BasketBall <- ais$hg[which(ais$sport=='B_Ball')] 

#

e ripetere poi l'operazione per altre nove volte, cambiando manualmente il nome dello sport (e della variabile nella quale memorizzare i dati estratti). Ottenuti così i dieci sottoinsiemi di dati, dovremmo poi applicare a ciascuno di essi, quindi per dieci volte, il test di normalità. 

Fortunatamente esiste una soluzione molto più semplice, che innanzitutto prevede di impiegare la funzione split() che applicata alla variabile hg della tabella ais (ais$hg) estrae e separa automaticamente i dati in base a tutti i valori assunti dal fattore sport (ais$sport), per farlo copiate e incollate nella Console di R questo codice e premete ↵ Invio.

#

# estrae e separa i valori di hg in base al fattore sport

#

hg_sport <- split(ais$hg, ais$sport) 

#

Per capire cosa è accaduto vediamo con la funzione str() la struttura dell'oggetto hg_sport così creato

#

# struttura dell'oggetto con i valori di hg separati in base al fattore sport 

#

str(hg_sport) 

#

che contiene quindi i valori di concentrazione dell'emoglobina, ma ora separati per ciascuno dei dieci sport praticati.

> str(hg_sport) 

List of 10

 $ B_Ball : num [1:25] 12.3 12.7 11.6 12.6 14 12.5 12.8 13.2 13.5 12.7 ...

 $ Field  : num [1:19] 13.3 14.7 15.3 14.3 14.6 15 15.2 15.4 16.7 15.5 ...

 $ Gym    : num [1:4] 13.4 14 12.5 14.5

 $ Netball: num [1:23] 13.6 12.7 12.3 12.3 12.8 11.8 12.7 12.4 12.4 14.1 ...

 $ Row    : num [1:37] 13.9 14.7 13.3 12.9 14.7 14.3 14.5 13 12.5 14.5 ...

 $ Swim   : num [1:22] 13.8 13.3 12.9 12.7 14.4 14 13.9 14 13.1 15.9 ...

 $ T_400m : num [1:29] 15.9 12.4 14.7 13.9 13 13.8 13.2 12.6 13.5 13.7 ...

 $ T_Sprnt: num [1:15] 13.4 14.4 14.7 14.2 15.6 15.2 15.5 15.5 14.9 19.2 ...

 $ Tennis : num [1:11] 12 13.9 13.5 12.1 14.8 14.5 13.9 17.7 14.3 14.9 ...

 $ W_Polo : num [1:17] 14.4 15.2 15 16.2 15.6 16.2 14.8 15.8 14.4 15.8 ...

Il lavoro di preparazione dei dati è stato rapido, ma la cosa più interessante è che ora diventa rapidissimo applicare l'analisi statistica a ciascuno dei dieci sport, è sufficiente infatti una sola riga di codice, copiatela e incollatela nella Console di R e premete ↵ Invio.

#

# test di normalità (gaussianità) di Shapiro-Wilk sui valori di hg per i dieci sport

#

sapply(hg_sport, shapiro.test) 

#

Come si vede con la funzione sapply() il test di Shapiro-Wilk (shapiro.test) per la normalità (gaussianità) viene automaticamente applicato a ciascuno dei dieci sottoinsiemi ricavati dalla scomposizione dei dati in base al fattore sport e contenuti nell'oggetto hg_sport. Da notare per inciso che la distribuzione dei valori di emoglobina non si allontana mai in modo significativo da una distribuzione gaussiana (in nessun caso abbiamo p <0.05)

> sapply(hg_sport, shapiro.test)

          B_Ball                        Field                        

statistic 0.9583881                     0.9740128                    

p.value   0.3833057                     0.8528473                    

method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"

data.name "X[[i]]"                      "X[[i]]"                     

          Gym                           Netball                      

statistic 0.9791217                     0.95727                      

p.value   0.896848                      0.4105184                    

method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"

data.name "X[[i]]"                      "X[[i]]"                     

          Row                           Swim                         

statistic 0.9782721                     0.9457531                    

p.value   0.6714365                     0.2595768                    

method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"

data.name "X[[i]]"                      "X[[i]]"                     

          T_400m                        T_Sprnt                      

statistic 0.9780738                     0.8952224                    

p.value   0.7872988                     0.08047871                   

method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"

data.name "X[[i]]"                      "X[[i]]"                     

          Tennis                        W_Polo                       

statistic 0.9399817                     0.9481069                    

p.value   0.520283                      0.4272774                    

method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"

data.name "X[[i]]"                      "X[[i]]"                

Gli stessi dati contenuti nell'oggetto hg_sport possono essere impiegati anche per generare (ad esempio) i grafici a scatola con i baffi delle dieci distribuzioni.

#

# boxplot dei valori di hg per i dieci sport

#

boxplot(hg_sport) 

#

Ed ecco il grafico, in una versione base ma adeguata (per vedere i nomi di tutte le variabili, dopo avere generato il grafico potete "afferrare" con il mouse il lato sinistro della finestra e allargarla).

L'estrazione e la separazione dei valori di concentrazione dell'emoglobina (ais$hg) in base al sesso (ais$sex) viene lasciata come esercizio. 

Ma non basta, perché al bisogno con la funzione split() è possibile separare i dati per più fattori: nel nostro caso possiamo separarli contemporaneamente per i due fattori sport e sesso.

#

# estrae e separa i valori di hg in base al fattore sport e al fattore sesso

#

hg_sport_sesso <- split(ais$hg, list(ais$sport, ais$sex), drop=TRUE) 

#

Le cose da fare in questo caso sono:

→ impiegare la funzione list() per elencare i fattori in base ai quali i dati devono essere suddivisi, nel nostro caso sia lo sport praticato (ais$sport) sia il sesso (ais$sex);

→ aggiungere drop=TRUE per scartare i sottoinsiemi vuoti, nel nostro caso invece di 20 sottoinsiemi (10 sport x 2 sessi) avremo solamente 17 sottoinsiemi perché dei dieci sport le donne ne praticano 9 e gli uomini ne praticano 8.

Ora con la funzione str() si può avere la conferma che nell'oggetto hg_sport_sesso si trovano i 17 sottoinsiemi separati sia per sport sia per sesso

> str(hg_sport_sesso)

List of 17

 $ B_Ball.f : num [1:13] 12.3 12.7 11.6 12.6 14 12.5 12.8 13.2 13.5 12.7 ...

 $ Field.f  : num [1:7] 13.3 14.7 15.3 14.3 14.6 15 15.2

 $ Gym.f    : num [1:4] 13.4 14 12.5 14.5

 $ Netball.f: num [1:23] 13.6 12.7 12.3 12.3 12.8 11.8 12.7 12.4 12.4 14.1 ...

 $ Row.f    : num [1:22] 13.9 14.7 13.3 12.9 14.7 14.3 14.5 13 12.5 14.5 ...

 $ Swim.f   : num [1:9] 13.8 13.3 12.9 12.7 14.4 14 13.9 14 13.1

 $ T_400m.f : num [1:11] 15.9 12.4 14.7 13.9 13 13.8 13.2 12.6 13.5 13.7 ...

 $ T_Sprnt.f: num [1:4] 13.4 14.4 14.7 14.2

 $ Tennis.f : num [1:7] 12 13.9 13.5 12.1 14.8 14.5 13.9

 $ B_Ball.m : num [1:12] 15.9 15.6 15.9 15.7 16.4 14.4 13.7 15.9 13.5 15 ...

 $ Field.m  : num [1:12] 15.4 16.7 15.5 14.7 18 16.1 15.9 16.3 15.8 15.7 ...

 $ Row.m    : num [1:15] 15 14.8 14.5 15.5 14.7 15.4 16.2 14.9 17.3 15.8 ...

 $ Swim.m   : num [1:13] 15.9 15.2 15.9 15 15.6 15.2 16.3 15.2 16.5 15.4 ...

 $ T_400m.m : num [1:18] 14.4 14 15.1 16.5 16.1 15 15.4 15.4 15.8 14.3 ...

 $ T_Sprnt.m: num [1:11] 15.6 15.2 15.5 15.5 14.9 19.2 15.1 14.9 17.2 16.5 ...

 $ Tennis.m : num [1:4] 17.7 14.3 14.9 15.7

 $ W_Polo.m : num [1:17] 14.4 15.2 15 16.2 15.6 16.2 14.8 15.8 14.4 15.8 ...

ed è quindi possibile, automaticamente per tutti i sottoinsiemi, effettuare il test di normalità e rappresentare i boxplot con le funzioni di cui sopra, semplicemente sostituendo hg_sport con hg_sport_sesso. 

La conclusione? Quando si inizia a lavorare con R molte volte ci si arena su problemi banali per i quali esistono soluzioni delle quali non si è a conoscenza data la vastità del linguaggio. Un esempio ne è l'estrazione dei dati per fattore da un set di dati mediante la funzione split(), che consente – indicando la variabile da estrarre e il fattore (o congiuntamente con list() i fattori) in base al quale separare i dati – di preparare i dati organizzandoli in un formato utile ai fini della loro successiva analisi statistica e della loro rappresentazione grafica.

----------

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

Barre dell'errore e limiti di confidenza

Questa volta vediamo come con ggplot2 sia possibile, impiegando una sola riga di codice, ottenere rappresentazioni multiple di una statistica (come ad esempio la media) con le corrispondenti barre dell'errore, che tecnicamente definiscono i limiti di confidenza della statistica.

Nel primo dei due esempi che seguono sottraiamo e sommiamo alla media di una variabile 1.96 volte la sua deviazione standard ottenendo i limiti di confidenza al 95% della distribuzione campionaria (la distribuzione dei dati del campione in esame).

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

# BARRE DELL'ERRORE E LIMITI DI CONFIDENZA della distribuzione campionaria

#

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

library(dplyr) # carica il pacchetto per i calcoli

library(ggplot2) # carica il pacchetto per la grafica

# 

# calcola media e deviazione standard separatamente per sport e per sesso

# 

ais_stat <- ais %>% group_by(sport, sex) %>% summarise(mean_hg=mean(hg), sd_hg=sd(hg)) %>% ungroup()

#

# crea il grafico ma alcuni dati risultano sovrapposti

#

ggplot(ais_stat, aes(x=sport, y=mean_hg, color=sex)) + geom_point(size=3) + geom_errorbar(aes(ymin=mean_hg-1.96*sd_hg, ymax=mean_hg+1.96*sd_hg), width=0.2) + labs(x="Sport praticato", y="Emoglobina in g/dL", color="Sesso dell'atleta") + theme_minimal() + theme(legend.position="top") 

#

Le prime tre righe servono a caricare i pacchetti necessari, che devono essere preventivamente scaricati dal CRAN, il primo dei quali contiene il set di dati ais impiegato come esempio [1]. 

La quarta riga di codice serve per preparare i dati impiegando l'operatore pipe (%>%) che trasferisce i dati da una funzione alla successiva in questo modo:

→ i dati (ais) sono trasferiti (%>%) alla funzione group_by();

→ la funzione group_by() raggruppa i dati in sottoinsiemi per sport praticato (sport) e per sesso dell'atleta (sex) e il risultato di questa suddivisione viene trasferito (%>%) alla funzione summarise();

→ la funzione summarise() calcola la media mean() e deviazione standard sd() della variabile hg che è la concentrazione nel sangue dell'emoglobina (espressa in g/dL, grammi per decilitro di sangue) dei 202 atleti australiani [1] suddivisi per sport e per sesso, e il risultato viene salvato (<-) nella tabella ais_stat;

→ i risultati sono trasferiti (%>%) alla funzione ungroup() che rimuove dalla tabella ais_stat una serie di valori impiegati provvisoriamente per generarla.

L'ultima riga di codice realizza un grafico:

→ impiegando le statistiche salvate nella variabile ais_stat;

→ riportando in ascisse lo sport praticato (sport), in ordinate la media dell'emoglobina (mean_hg), assegnando colori diversi in base al sesso (sex);

→ riportando un grafico a punti (geom_point) dei dati;

→ aggiungendo con geom_errorbar() le barre che indicano l'errore che nel nostro caso è rappresentato da 1.96 volte la deviazione standard (sd_hg);

→ aggiungendo con labs() le opportune etichette; 

→ impiegando con theme_minimal() una rappresentazione semplice;

→ posizionando infine la legenda in alto (legend.position="top").

Questo è il grafico.

Poiché le barre dell'errore nei due sessi risultano sovrapposte, riducendo la chiarezza del grafico, riprendiamo lo stesso codice con una piccola modifica. Copiate e incollate nella Console di R questo script e premete ↵ Invio.

# BARRE DELL'ERRORE E LIMITI DI CONFIDENZA della distribuzione campionaria

#

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

library(dplyr) # carica il pacchetto per i calcoli

library(ggplot2) # carica il pacchetto per la grafica

# 

# calcola media e deviazione standard separatamente per sport e per sesso

# 

ais_stat <- ais %>% group_by(sport, sex) %>% summarise(mean_hg=mean(hg), sd_hg=sd(hg)) %>% ungroup()

#

# crea il grafico evitando la sovrapposizione dei dati

ggplot(ais_stat, aes(x=sport, y=mean_hg, color=sex)) + geom_point(size=3, position=position_dodge(width=0.5)) + geom_errorbar(aes(ymin=mean_hg-1.96*sd_hg, ymax=mean_hg+1.96*sd_hg), width=0.2, position=position_dodge(width=0.5)) + labs(x="Sport praticato", y="Emoglobina in g/dL", color="Sesso dell'atleta") + theme_minimal() + theme(legend.position="top") 

#

L'unica differenza è che aggiungiamo alla funzione geom_point() e alla funzione geom_errorbar() l'argomento position=position_dodge(width=0.5): questo consente di eliminare la sovrapposizione delle barre dell'errore. 

Se ora digitate ais_stat nella Console di R potete visualizzare i dati impiegati per realizzare i due grafici precedenti: le medie (mean_hg) e le deviazioni standard (sd_hg), suddivise per ciascuno sport (sport) e per ciascun sesso (sex).

> ais_stat

# A tibble: 17 × 4

   sport   sex   mean_hg sd_hg

   <fct>   <fct>   <dbl> <dbl>

 1 B_Ball  f        13.1 0.878

 2 B_Ball  m        15.1 0.922

 3 Field   f        14.6 0.682

 4 Field   m        16.0 0.805

 5 Gym     f        13.6 0.860

 6 Netball f        12.8 0.567

 7 Row     f        14.0 0.740

 8 Row     m        15.4 0.711

 9 Swim    f        13.6 0.583

10 Swim    m        15.5 0.655

11 T_400m  f        13.8 1.04 

12 T_400m  m        15.3 0.824

13 T_Sprnt f        14.2 0.556

14 T_Sprnt m        16.2 1.49 

15 Tennis  f        13.5 1.10 

16 Tennis  m        15.6 1.48 

17 W_Polo  m        15.5 0.718

Ripetiamo l'analisi dei dati: ma questa volta sottraiamo e sommiamo alla media della variabile 1.96 volte il suo errore standard ottenendo i limiti di confidenza al 95% della media campionaria.

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

# BARRE DELL'ERRORE E LIMITI DI CONFIDENZA della media campionaria

#

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

library(dplyr) # carica il pacchetto per i calcoli

library(ggplot2) # carica il pacchetto per la grafica

library(plotrix) # carica il pacchetto per il calcolo dell'errore standard

# 

# calcola media e deviazione standard separatamente per sport e per sesso

# 

ais_stat <- ais %>% group_by(sport, sex) %>% summarise(mean_hg=mean(hg), es_hg=std.error(hg)) %>% ungroup()

#

# crea il grafico evitando la sovrapposizione dei dati

#

ggplot(ais_stat, aes(x=sport, y=mean_hg, color=sex)) + geom_point(size=3, position=position_dodge(width=0.5)) + geom_errorbar(aes(ymin=mean_hg-1.96*es_hg, ymax=mean_hg+1.96*es_hg), width=0.2, position=position_dodge(width=0.5)) + labs(x="Sport praticato", y="Emoglobina in g/dL", color="Sesso dell'atleta") + theme_minimal() + theme(legend.position="top") 

#

Da notare che deve essere installato anche il pacchetto plotrix che viene impiegato per calcolare l'errore standard della media. Per il resto viene impiegato esattamente lo stesso codice dello script precedente, semplicemente sostituendo la deviazione standard dell'emoglobina con l'errore standard della media (es_hg) calcolato con la funzione std.error(). Questo è il risultato. 

Come si rileva dal grafico, a parte il caso del tennis (Tennis), all'interno di ciascuno sport i limiti di confidenza al 95% delle medie dell'emoglobina negli atleti dei due sessi non si sovrappongono (la considerazione che se ne può trarre è riportata qui sotto nella Conclusione).

Se digitate ais_stat nella Console di R potete ora visualizzare la tabella che contiene media (mean_hg) ed errore standard (es_hg) calcolati per ciascuno sport (sport) e per ciascun sesso (sex).

> ais_stat

# A tibble: 17 × 4

   sport   sex   mean_hg es_hg

   <fct>   <fct>   <dbl> <dbl>

 1 B_Ball  f        13.1 0.243

 2 B_Ball  m        15.1 0.266

 3 Field   f        14.6 0.258

 4 Field   m        16.0 0.232

 5 Gym     f        13.6 0.430

 6 Netball f        12.8 0.118

 7 Row     f        14.0 0.158

 8 Row     m        15.4 0.184

 9 Swim    f        13.6 0.194

10 Swim    m        15.5 0.182

11 T_400m  f        13.8 0.312

12 T_400m  m        15.3 0.194

13 T_Sprnt f        14.2 0.278

14 T_Sprnt m        16.2 0.450

15 Tennis  f        13.5 0.414

16 Tennis  m        15.6 0.741

17 W_Polo  m        15.5 0.174

Conclusione: 

→ il termine barra dell'errore o errorbar indica semplicemente l'espediente grafico impiegato per rappresentare attorno al valore di una statistica i limiti di confidenza della statistica al livello di probabilità desiderato. In genere si impiega il livello di probabilità del 95% in modo che la probabilità che il valore della statistica – stimato mediante i dati campionari – possa uscire dai limiti di confidenza così calcolati sia "solamente" del 5%;

→ i limiti di confidenza al 95% della distribuzione campionaria – calcolati come media±1.96 ∙ ds (deviazione standard) – sono i limiti all'interno dei quali cade il 95% dei valori osservati nel campione esaminato, quindi la probabilità che un nuovo dato risulti al di fuori di questi limiti è del 5%;

→ i limiti di confidenza al 95% della media campionaria – calcolati come media±1.96 ∙ es (errore standard) – sono i limiti all'interno dei quali abbiamo il 95% di probabilità che cada la media, stante l'informazione contenuta nel nostro campione. Questi limiti forniscono un test di significatività della differenza tra le medie: quando i limiti di confidenza al 95% di due medie non si sovrappongono, possiamo affermare che tra le due medie esiste una differenza statisticamente significativa. Quindi nel nostro caso possiamo concludere che all'interno di ciascuno degli sport praticati da entrambi i sessi, con una sola eccezione, la concentrazione media dell'emoglobina nelle atlete è significativamente inferiore a quella degli atleti e che la probabilità che questa conclusione sia errata è del 5% o meno;

→ il risultato di una statistica dovrebbe essere sempre riportato con la sua incertezza, espressa in termini dei suoi limiti di confidenza opportunamente specificati: limiti di confidenza al 95%, limiti di confidenza al 99%, media±n∙ds, media±n∙es, ove n è l'opportuno moltiplicatore, eccetera.

----------

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