lunedì 7 gennaio 2019

Analisi esplorativa dei dati

Con l'espressione “analisi esplorativa dei dati” non si fa riferimento a una tecnica statistica specifica, bensì all'insieme delle valutazioni preliminari che è sempre necessario effettuare nell'ambito di un percorso logico che prevede, per il calcolo delle statistiche elementari di una singola variabile (analisi univariata), i seguenti passi:
→ analisi esplorativa dei dati;
→ esecuzione dei test di normalità (gaussianità), per valutare se i dati seguono una distribuzione gaussiana;
→ calcolo delle statistiche elementari parametriche (media, deviazione standard, varianza, quantili parametrici) se i dati seguono una distribuzione gaussiana [1];
→ calcolo delle statistiche elementari non parametriche (mediana, deviazione assoluta mediana o MAD, quartili e altri quantili non parametrici) se i dati non seguono una distribuzione gaussiana.

In quanto si tratta di una prima fase, l'analisi esplorativa dei dati non prevede test di significatività, che sono riservati alla fasi successive, ma rappresenta piuttosto il momento di valutazione critica preliminare e globale dei dati raccolti, che deve includere quantomeno:
→ l'identificazione di eventuali dati mancanti;
→ un primo confronto orientativo tra i risultati di statistiche parametriche e di statistiche non parametriche;
→ l'individuazione di possibili dati anomali (outliers) cioè di dati che si discostano in modo inatteso dagli altri;
 l'identificazione della possibile origine degli outliers (errori di digitazione? casi inclusi erroneamente? problemi strumentali? altro?) e la valutazione degli interventi correttivi (laddove possibili).

Vediamo ora alcune funzioni che possono aiutarci per le attività previste nei primi tre punti (ovviamente non per il quarto) utilizzando come esempio i dati ematologici e biometrici rilevati in 202 atleti australiani contenuti nella tabella ais del pacchetto DAAG - accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [2] dove trovate anche illustrati i dati contenuti nella tabella.

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

# ANALISI ESPLORATIVA DEI DATI funzioni base
#
library(DAAG) # carica il pacchetto DAAG che include il set di dati ais
#
ais[!complete.cases(ais)] # verifica la presenza di NA (dati mancanti)
colSums(is.na(ais)) # conteggia gli eventuale dati mancanti per colonna
#
summary(ais) # statistiche elementari per tutte le variabili (numeriche e non) di ais
#
media <- mean(ais$ferr) # calcola la media della ferritina
mediana <- median(ais$ferr) # calcola la mediana
data.frame(media, mediana) # le mette a confronto
#
ds <- sd(ais$ferr) # calcola la deviazione standard della ferritina
mad <- mad(ais$ferr) # calcola la Median Absolute Deviation (about the median) o MAD
data.frame(ds, mad) # le mette a confronto
#
qpar <- round(qnorm(c(seq(0.025, 0.975, 0.025)), mean=mean(ais$ferr), sd=sd(ais$ferr)), digits=2) # calcola i quantili parametrici della ferritina
qnon <- round(quantile(ais$ferr, probs=seq (0.025, 0.975, 0.025)), digits=2) # calcola i quantili non parametrici
data.frame(qpar, qnon) # li mette a confronto
#
windows() # apre e inizializza una nuova finestra grafica
boxplot(ais$ferr, range=1.5, horizontal=FALSE, main="Valori oltre la mediana ± 1.5 · IQR", ylab="Ferritina in µg/L", notch=FALSE, col="yellow") # rappresenta graficamente i valori della ferritina
boxplot.stats(ais$ferr, coef=1.5)$stats # mostra i 5 punti notevoli dei baffi e della scatola
boxplot.stats(ais$ferr, coef=1.5)$out # mostra i valori che si trovano oltre la mediana ± 1.5 volte il range interquartile (IQR)
#

Con la prima riga di codice viene caricato il pacchetto che contiene la tabella ais con i dati che ci servono. 

Quindi con la funzione complete.cases() viene fatta la ricerca dei casi con dati non (!) completi, nei quali dovremmo avere i dati mancanti sostituiti con la sigla NA (Not Available). 

> ais[!complete.cases(ais)] # verifica la presenza di NA (dati mancanti)
data frame con 0 colonne e 202 righe

Con la funzione colSums() sono ricercati i dati mancanti is.na() per ciascuna delle colonne/variabili della tabella:

> colSums(is.na(ais)) # conteggia gli eventuale dati mancanti per colonna
   rcc    wcc     hc     hg   ferr    bmi    ssf pcBfat    lbm     ht     wt 
     0      0      0      0      0      0      0      0      0      0      0 
   sex  sport 
     0      0 

Ora sappiamo che nella tabella i dati sono completi, non abbiamo dati mancanti.

Segue la funzione summary() che per ciascuna variabile della tabella riporta valore minimo (Min.), primo quartile (1st Qu.), mediana o secondo quartile (Median), media (Mean), terzo quartile (3rd Qu.), valore massimo (Max.). 

> summary(ais) # statistiche elementari per tutte le variabili (numeriche e non) di ais
      rcc             wcc               hc              hg       
 Min.   :3.800   Min.   : 3.300   Min.   :35.90   Min.   :11.60  
 1st Qu.:4.372   1st Qu.: 5.900   1st Qu.:40.60   1st Qu.:13.50  
 Median :4.755   Median : 6.850   Median :43.50   Median :14.70  
 Mean   :4.719   Mean   : 7.109   Mean   :43.09   Mean   :14.57  
 3rd Qu.:5.030   3rd Qu.: 8.275   3rd Qu.:45.58   3rd Qu.:15.57  
 Max.   :6.720   Max.   :14.300   Max.   :59.70   Max.   :19.20  
                                                                 
      ferr             bmi             ssf             pcBfat      
 Min.   :  8.00   Min.   :16.75   Min.   : 28.00   Min.   : 5.630  
 1st Qu.: 41.25   1st Qu.:21.08   1st Qu.: 43.85   1st Qu.: 8.545  
 Median : 65.50   Median :22.72   Median : 58.60   Median :11.650      
 Mean   : 76.88   Mean   :22.96   Mean   : 69.02   Mean   :13.507         
 3rd Qu.: 97.00   3rd Qu.:24.46   3rd Qu.: 90.35   3rd Qu.:18.080  
 Max.   :234.00   Max.   :34.42   Max.   :200.80   Max.   :35.520  
                                                                   
      lbm               ht              wt         sex         sport   
 Min.   : 34.36   Min.   :148.9   Min.   : 37.80   f:100   Row    :37  
 1st Qu.: 54.67   1st Qu.:174.0   1st Qu.: 66.53   m:102   T_400m :29  
 Median : 63.03   Median :179.7   Median : 74.40           B_Ball :25  
 Mean   : 64.87   Mean   :180.1   Mean   : 75.01           Netball:23  
 3rd Qu.: 74.75   3rd Qu.:186.2   3rd Qu.: 84.12           Swim   :22  
 Max.   :106.00   Max.   :209.4   Max.   :123.20           Field  :19  
                                                           (Other):47  

In una distribuzione gaussiana media e mediana sono identiche [1]. Qui si osservano tra le due in genere valori abbastanza simili, tranne che nel caso della ferritina (ferr), della somma dello spessore delle pliche cutanee (ssf) e della percentuale di grasso corporeo (pcBfat).

Continuiamo quindi con l'analisi esplorativa dei dati, che per semplicità è qui limitata alla ferritina, calcolando e confrontando media e mediana

> media <- mean(ais$ferr) # calcola la media della ferritina
> mediana <- median(ais$ferr) # calcola la mediana
> data.frame(media, mediana) # le mette a confronto
     media mediana
1 76.87624    65.5

calcolando e confrontando deviazione standard e MAD

> ds <- sd(ais$ferr) # calcola la deviazione standard della ferritina
> mad <- mad(ais$ferr) # calcola la Median Absolute Deviation (about the median) o MAD
> data.frame(ds, mad) # le mette a confronto
        ds     mad
1 47.50124 37.8063

e infine calcolando e confrontando quantili parametrici e quantili non parametrici

> qpar <- round(qnorm(c(seq (0.025, 0.975, 0.025)), mean = mean(ais$ferr), sd = sd(ais$ferr)), digits=2) # calcola i quantili parametrici della ferritina
> qnon <- round(quantile(ais$ferr, probs = seq (0.025, 0.975, 0.025)), digits=2) # calcola i quantili non parametrici
> data.frame(qpar, qnon) # li mette a confronto
        qpar   qnon
2.5%  -16.22  20.00
5%     -1.26  22.00
7.5%    8.50  27.15
10%    16.00  30.00
12.5%  22.23  33.12
15%    27.64  35.15
17.5%  32.48  37.18
20%    36.90  39.20
22.5%  40.99  41.00
25%    44.84  41.25
27.5%  48.48  43.00
30%    51.97  44.00
32.5%  55.32  48.00
35%    58.57  50.00
37.5%  61.74  53.00
40%    64.84  55.00
42.5%  67.89  58.00
45%    70.91  59.45
47.5%  73.90  62.48
50%    76.88  65.50
52.5%  79.85  68.53
55%    82.85  71.00
57.5%  85.86  73.00
60%    88.91  76.00
62.5%  92.01  80.00
65%    95.18  85.65
67.5%  98.43  87.68
70%   101.79  90.70
72.5% 105.27  93.73
75%   108.92  97.00
77.5% 112.76 102.00
80%   116.85 107.00
82.5% 121.27 114.13
85%   126.11 122.00
87.5% 131.52 125.88
90%   137.75 138.40
92.5% 145.26 155.93
95%   155.01 182.95
97.5% 169.98 212.00

Nel caso di una distribuzione gaussiana tutti questi valori devono essere identici [1] e qui non lo sono.

Per l'identificazione degli outliers impieghiamo la funzione boxplot() realizzando un grafico a scatola con i baffi (boxplot) nel quale sono riportati come singoli punti separati i dati che si trovano oltre la mediana ± 1.5 volte il range interquartile (argomento range=1.5) [4].


Con le due ultime righe dello script sono quindi presentati due risultati della funzione boxplot.stats()

Con $stats sono riportati i valori delle ferritina corrispondenti:
→ al baffo inferiore (8.0)
→ al margine inferiore della scatola (41.0)
→ alla mediana (65.5)
→ al margine superiore della scatola (97.0)
→ al baffo superiore (177.0)

> boxplot.stats(ais$ferr, coef=1.5)$stats # mostra i 5 punti notevoli dei baffi e della scatola
[1]   8.0  41.0  65.5  97.0 177.0

Con $out sono riportati i dati anomali (182 183 212 213 184 220 191 189 212 234 214 233), o outliers. 

> boxplot.stats(ais$ferr, coef=1.5)$out # mostra i valori che si trovano oltre la mediana ± 1.5 volte il range interquartile (IQR)
 [1] 182 183 212 213 184 220 191 189 212 234 214 233

Si tratta dei dati che nel grafico si trovano oltre la mediana ± 1.5 volte il range interquartile (argomento coef=1.5), che non per questo devono essere esclusi dalle statistiche, ma sui quali sarebbe opportuno, dato il loro eccessivo scostamento dai dati rimanenti, effettuare una rivalutazione (cosa qui ovviamente impraticabile non avendo accesso alla complessa catena di eventi che ha portato alla produzione dei dati).

Oltre alla funzione summary() inclusa nel pacchetto base di R [5] altre funzioni [6] per l'analisi esplorativa dei dati sono disponibili nei pacchetti Hmiscpastecspsych.

Potete scaricare e installare i pacchetti aggiuntivi dal CRAN e trovate la loro documentazione completa, incluso il manuale di riferimento, sul repository della documentazione di R [7].

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

ANALISI ESPLORATIVA DEI DATI funzioni per una analisi globale
#
library(DAAG) # carica il pacchetto DAAG incluso il set di dati ais
mydata <- ais[c(1,2,3,4,5,6,7,8,9,10,11)] # salva le colonne con le sole variabili numeriche in mydata
#
summary(mydata) # statistiche elementari per tutte le variabili 
#
library(Hmisc) # carica il pacchetto
describe(mydata) # statistiche del pacchetto Hmisc
#
library(pastecs) # carica il pacchetto
stat.desc(mydata) # statistiche del pacchetto pastecs
#
library(psych) # carica il pacchetto
describe(mydata) # statistiche del pacchetto psych
describeBy(mydata, ais$sex) # statistiche del pacchetto psych separate per sesso
describeBy(mydata, ais$sport) # statistiche del pacchetto psych separate per sport
#

Le funzioni espandono, ciascuna a modo suo,  il quadro dei dati fornito dalla funzione summary() del pacchetto base di R, e sono state riportate perché ciascuno ne possa valutare i risultati e l'utilità per gli scopi che si propone.

L'aspetto forse più interessante è rappresentato dalla possibilità offerta dalla funzione describeBy() del pacchetto psych di riportare statistiche riepilogative riaggregando i dati in base ai fattori presenti nei record, e quindi, nel caso specifico, di riportare statistiche separate per sesso (m,f), come pure statistiche separate per ciascuno degli sport praticati dagli atleti (B_Ball, Field, Gym, Netball, Row, Swim, T_400m, T_Sprnt, Tennis, W_Polo).

Conclusione: l'analisi esplorativa dei dati fornisce informazioni utili a evidenziare per la ferritina alcuni problemi che meritano di essere approfonditi, valutati e opportunamente gestiti continuando a seguire le fasi del percorso logico riportato all'inizio. 

Potete riutilizzare facilmente lo script sostituendo all'oggetto ais l'oggetto contenente i vostri dati, opportunamente strutturati. Per una guida rapida all'importazione dei dati potete consultare i link:
 importare i dati di un file .csv
 importare i dati di un file .xls o .xlsx
 gestione dei dati mancanti


----------

[1] Una tipica distribuzione gaussiana è riportata nel post: La distribuzione gaussiana.

[2] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG. Per manuale di riferimento del pacchetto vedere nel repository della documentazione: Package 'DAAG'. URL consultato il 20/01/2023: https://goo.gl/gzKabK

[3] La “Median Absolute Deviation (about median)” o MAD ovvero la deviazione assoluta mediana (dalla mediana) è l'equivalente non parametrico della deviazione standard e in una distribuzione gaussiana ha lo stesso valore di questa. Vedere: Rousseeuw PJ, Croux C. Alternatives to the Median Absolute Deviation. Journal of the American Statistical Association 88 (424), 1273-1283, 1993. URL consultato il 04/01/2019: https://goo.gl/4Rh53b

[4] Per i dettagli sulla funzione boxplot() digitare help(boxplot) nella Console di R e vedere il post Grafici a scatola con i baffi.

[5] Vedere il manuale di riferimento del pacchetto base R: A Language and Environment for Statistical Computing, Reference Index. URL consultato il 28/11/2018: https://goo.gl/kes91X

[6] Per la loro documentazione digitate come al solito help(nomedellafunzione) nella Console di R.

[7] Available CRAN Packages By NameURL consultato il 20/10/2018: https://goo.gl/hLC9BB

Nessun commento:

Posta un commento