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 (13 variables) 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: prima di procedere con l'analisi statistica o grafica, 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')] 

#

quindi dovremmo convertire i valori estratti in valori numerici

#

BasketBall <- as.numeric(BasketBall) 

#

e infine potremmo applicare il test di Shapiro-Wilk

#

shapiro.test(BasketBall) 

#

ottenendo questo risultato 

        Shapiro-Wilk normality test

data:  BasketBall

W = 0.95839, p-value = 0.3833

E quanto sopra dovrebbe essere ripetuto per altre nove volte, cambiando manualmente il nome dello sport (e possibilmente per chiarezza anche quello della variabile nella quale i dati vengono estratti). 

Fortunatamente esiste una strada molto più semplice: innanzitutto si preparano i dati impiegando la funzione split(), che consente di estrarre e separare automaticamente i valori della variabile hg della tabella ais (ais$hg) in base a tutti i valori assunti dal fattore sport (ais$sport), per farlo copiate e incollate nella Console di R questa riga di 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 ancora tutti i valori di concentrazione dell'emoglobina del set di dati ais, ma ora separati per ciascuno dei dieci sport praticati, con l'indicazione tra parentesi quadra del numero di dati per ciascuno di essi.

> 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 rapido anche applicare il test di normalità 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 valori di concentrazione dell'emoglobina 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 dei valori di concentrazione dell'emoglobina.

#

# 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).

Ma non basta, perché al bisogno con la funzione split() è possibile separare i dati per più fattori: nel nostro caso possiamo separare i valori di concentrazione dell'emoglobina 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 ricavati dalla scomposizione dei valori di concentrazione dell'emoglobina 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. 

Come impiegare l'operatore %>% e l'argomento FUN

Anche se l'obiettivo di questi post rimane sempre quello di fornire gli strumenti più semplici, quelli utili nella fase iniziale di apprendimento – notoriamente tutta in salita e nella quale quello che conta è "rompere il ghiaccio" con R – occasionalmente può essere interessante affrontare argomenti un poco più tecnici, come i due che vediamo ora.

L'operatore %>%

Il simbolo %>% viene impiegato per indicare l'operatore pipe: il termine sta per il sostantivo "tubo" ma anche per l'espressione verbale "trasportare con una tubatura", ed è proprio in questo senso che viene impiegato. L'operatore pipe è stato introdotto la prima volta con il pacchetto magrittr [1] che ne riporta la seguente descrizione: 

"[il pacchetto] fornisce un meccanismo per concatenare i comandi con il nuovo operatore di inoltro pipe, %>%. Questo operatore inoltrerà un valore, o il risultato di un'espressione, alla successiva chiamata/espressione di funzione. [l'operatore] fornisce un supporto flessibile per il tipo di espressioni sul lato destro. Per ulteriori informazioni, consultare la descrizione del pacchetto. Per citare René Magritte, 'Ceci n'est pas une pipe'. " [2]

In altre parole: l'operatore %>% fornisce la "tubatura" che consente – in modo completamente trasparente per il programmatore – di inoltrare i dati in uscita da una funzione all'ingresso della funzione che segue. Vediamo questo breve esempio. Abbiamo il set di dati ais [3] 

> 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

…...

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

 

e vogliamo calcolare la media, la deviazione standard e l'errore standard della concentrazione di emoglobina (hg) nel sangue dei 202 atleti australiani ivi inclusi suddividendoli per sport praticato (sport) e per sesso (sex).

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

# COME IMPIEGARE l'operatore %>% (pipe)

#

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

library(dplyr) # carica il pacchetto con l'operatore pipe e altre funzioni

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

# 

# calcola media, deviazione standard ed errore standard separatamente per sport e per sesso

# 

pipe_hg <- ais %>% group_by(sport, sex) %>% summarise(media=mean(hg), ds=sd(hg), es=std.error(hg)) %>% ungroup()

#

pipe_hg # mostra i risultati

#

Dopo avere installato e caricato i pacchetti necessari, per realizzare la tabella che contiene i risultati desiderati impieghiamo più volte 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() – che raggruppa i dati in sottoinsiemi per sport praticato (sport) e per sesso dell'atleta (sex) – trasferisce il risultato di questa operazione (%>%) alla funzione summarise();

→ la funzione summarise() calcola la media mean(), la deviazione standard sd() e l'errore standard std.error() della variabile hg – la concentrazione nel sangue dell'emoglobina (espressa in g/dL, grammi per decilitro di sangue) dei 202 atleti australiani suddivisi per sport e per sesso – e trasferisce il risultato di questa operazione (%>%) alla funzione ungroup();

→ la funzione ungroup() rimuove una serie di valori impiegati provvisoriamente per generare il risultato finale, che viene infine salvato (<-) nella tabella pipe_hg.

Ed ecco la tabella risultante, ottenuta quindi con una sola riga di codice:

> pipe_hg # mostra i risultati

# A tibble: 17 × 5

   sport   sex   media    ds    es

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

 1 B_Ball  f      13.1 0.878 0.243

 2 B_Ball  m      15.1 0.922 0.266

 3 Field   f      14.6 0.682 0.258

 4 Field   m      16.0 0.805 0.232

 5 Gym     f      13.6 0.860 0.430

 6 Netball f      12.8 0.567 0.118

 7 Row     f      14.0 0.740 0.158

 8 Row     m      15.4 0.711 0.184

 9 Swim    f      13.6 0.583 0.194

10 Swim    m      15.5 0.655 0.182

11 T_400m  f      13.8 1.04  0.312

12 T_400m  m      15.3 0.824 0.194

13 T_Sprnt f      14.2 0.556 0.278

14 T_Sprnt m      16.2 1.49  0.450

15 Tennis  f      13.5 1.10  0.414

16 Tennis  m      15.6 1.48  0.741

17 W_Polo  m      15.5 0.718 0.174

L'argomento FUN

L'argomento FUN consente di specificare una o più funzioni e di applicarle in sequenza a una variabile. Per illustrarlo impieghiamo gli stessi dati dell'esempio precedente allo stesso scopo: calcolare la media, la deviazione standard e l'errore standard della concentrazione di emoglobina nel sangue di 202 atleti australiani suddividendoli per sport praticato e per sesso. 

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

# COME IMPIEGARE l'argomento FUN

#

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

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

#

# calcola media, deviazione standard ed errore standard separatamente per sport e per sesso

#

FUN_hg <- aggregate(hg~sport+sex, data=ais, FUN=function(x) c(media=mean(x), ds=sd(x), es=std.error(x))) 

#

FUN_hg # mostra i risultati 

#

Dopo avere installato e caricato i pacchetti necessari, la tabella che contiene i risultati desiderati viene realizzata con la funzione aggregate() che nell'ordine:

→ nel nostro set di dati (data=ais) raggruppa (~) i valori di concentrazione dell'emoglobina (hg) in sottoinsiemi per sport praticato (sport) e per sesso dell'atleta (sex);

→ con l'argomento FUN e con il valore =function(x) specifica che sui dati devono essere effettuati i calcoli di seguito elencati;

→ con la funzione c() elenca le funzioni da impiegare per i calcoli;

→ con le funzioni mean(), sd() e std.error() effettua i calcoli di media, deviazione standard ed errore standard sui dati raggruppati per sport e per sesso, salvando i risultati nelle variabili media, ds ed es. 

I risultati sono infine riportati (<-) nella tabella FUN_hg e, a parte l'ordinamento e gli arrotondamenti, sono identici a quelli ottenuti con lo script precedente, impiegando anche in questo caso una sola riga di codice.

> FUN_hg # mostra i risultati 

     sport sex   hg.media      hg.ds      hg.es

1   B_Ball   f 13.1307692  0.8778616  0.2434750

2    Field   f 14.6285714  0.6824326  0.2579353

3      Gym   f 13.6000000  0.8602325  0.4301163

4  Netball   f 12.8173913  0.5670114  0.1182301

5      Row   f 14.0318182  0.7396179  0.1576871

6     Swim   f 13.5666667  0.5830952  0.1943651

7   T_400m   f 13.7727273  1.0354621  0.3122036

8  T_Sprnt   f 14.1750000  0.5560276  0.2780138

9   Tennis   f 13.5285714  1.0965313  0.4144499

10  B_Ball   m 15.1416667  0.9219134  0.2661335

11   Field   m 16.0250000  0.8046738  0.2322893

12     Row   m 15.3866667  0.7110020  0.1835799

13    Swim   m 15.5076923  0.6550592  0.1816807

14  T_400m   m 15.3055556  0.8235306  0.1941080

15 T_Sprnt   m 16.1909091  1.4936228  0.4503442

16  Tennis   m 15.6500000  1.4821156  0.7410578

17  W_Polo   m 15.5176471  0.7178399  0.1741017

Conclusione: nel corso dell'elaborazione dei dati può essere utile tenere a portata di mano qualche esempio di codice, come questi due, che ricorda come impiegare l'operatore %>% e l'argomento FUN. Non trovate molto altro sul web, ma potete provare ricercando ad esempio "R %>% operator examples" o "R FUN function examples".

-----------

[1] Vedere il Reference manual del pacchetto: Package ‘magrittr’.

https://cran.r-project.org/web/packages/magrittr/magrittr.pdf

[2] "Provides a mechanism for chaining commands with a new forward-pipe operator, %>%. This operator will forward a value, or the result of an expression, into the next function call/expression. There is flexible support for the type of right-hand side expressions. For more information, see package vignette. To quote Rene Magritte, 'Ceci n'est pas un pipe.' ".

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