Gestione dei dati mancanti

Per illustrare il tema, che è strettamente collegato all'ordinamento dei dati [1], vediamo tre cose:

→ come individuare i dati mancanti;

→ come i dati mancanti - identificati in R con la sigla NA (Not Available) - possono bloccare l'esecuzione di calcoli sulle variabili numeriche che li contengono e come in questo caso sia possibile omettere selettivamente dai calcoli i dati NA; 

→ come al bisogno sia possibile eliminare per intero da un vettore o matrice o tabella [2] i casi con dati NA.

Per proseguire dovete, seguendo le istruzioni fornite alla pagina Dati:

→ effettuare il download del file importa_csv.csv

→ salvare il file nella cartella C:\Rdati\

In alternativa copiate le otto righe riportate qui sotto, incollatele in un editor di file di testo aggiungendo un ↵ Invio al termine dell'ultima riga, e salvate il tutto in C:\Rdati\ in un file di testo denominato importa_csv.csv (attenzione: assicuratevi che il file sia effettivamente salvato con l'estensione .csv):

id;sesso;anni;peso_kg;altezza_m

MT;M;69;76;1,78

GF;F;56;63;

MC;F;53;71;1,60

SB;M;28;73;1,78

FE;F;61;54;1,54

AB;M;46;92;1,84

RF;F;31;81;1,56

Vediamo come individuare i dati mancanti. Copiate e incollate nella Console di R questo script e premete ↵ Invio:

# IDENTIFICA E CONTEGGIA I DATI MANCANTI 

#

# importa i dati, notare / invece di \ su windows

mydata <- read.table("C:/Rdati/importa_csv.csv", header=TRUE, sep=";", dec=",", row.names="id")

#

mydata # mostra i dati importati

#

is.na(mydata) # identifica i dati mancanti

#

colSums(is.na(mydata)) # conteggia i dati mancanti per colonna/variabile

#

which(is.na(mydata$altezza_m)) # identifica la posizione dei dati mancanti

#

Dopo avere importato i dati, richiamando il nome dell'oggetto mydata che li contiene questi sono mostrati con la sigla NA riportata automaticamente da R in corrispondenza di ciascuno dei dati mancanti.

> mydata # mostra i dati importati

   sesso anni peso_kg altezza_m

MT     M   69      76      1.78

GF     F   56      63        NA

MC     F   53      71      1.60

SB     M   28      73      1.78

FE     F   61      54      1.54

AB     M   46      92      1.84

RF     F   31      81      1.56

Nel caso di piccole tabelle come questa non serve altro: ma nel caso di estesi database ci servono funzioni che forniscano qualche soluzione più pratica e immediata. 

Impieghiamo quindi la funzione is.na() che identifica nella tabella mydata i dati mancanti riportando TRUE in corrispondenza di ciascuno di essi 

> is.na(mydata) # identifica i dati mancanti

   sesso  anni peso_kg altezza_m

MT FALSE FALSE   FALSE     FALSE

GF FALSE FALSE   FALSE      TRUE

MC FALSE FALSE   FALSE     FALSE

SB FALSE FALSE   FALSE     FALSE

FE FALSE FALSE   FALSE     FALSE

AB FALSE FALSE   FALSE     FALSE

RF FALSE FALSE   FALSE     FALSE

I dati della tabella mydata trasformati nei corrispettivi valori logici con la funzione is.na() possono allora essere impiegati come argomento della funzione colSums() che effettua il conteggio dei dati mancanti per ciascuna colonna/variabile

colSums(is.na(mydata)) # conteggia i dati mancanti per colonna/variabile

    sesso      anni   peso_kg altezza_m 

        0         0         0         1 

A questo punto abbiamo individuato in altezza_m la variabile con i dati mancanti e abbiamo conteggiato il loro numero. 

Non ci resta quindi che impiegare i dati della variabile mydata$altezza_m trasformati nei corrispettivi valori logici con is.na() come argomento della funzione which() per identificare la posizione dei  dati mancanti

> which(is.na(mydata$altezza_m)) # identifica la posizione dei dati mancanti

[1] 2

che nel nostro caso risultano essere il dato della riga numero 2 della variabile mydata$altezza_m.

Per l'effetto che i dati mancanti possono determinare sui calcoli copiate e incollate nella Console di R questo script e premete ↵ Invio:

# ESEMPIO DI MEDIA ARITMETICA NON COMPUTABILE A CAUSA DI DATI MANCANTI 

#

# importa i dati, notare / invece di \ su windows

mydata <- read.table("C:/Rdati/importa_csv.csv", header=TRUE, sep=";", dec=",", row.names="id")

#

mydata # mostra i dati importati

#

mean(mydata$peso_kg) # calcola la media aritmetica dei valori di peso

mean(mydata$altezza_m) # calcola la media aritmetica dei valori di altezza

#

Il risultato che compare nella Console di R è il seguente:

> mydata # mostra i dati importati 

   sesso anni peso_kg altezza_m

MT     M   69      76      1.78

GF     F   56      63        NA

MC     F   53      71      1.60

SB     M   28      73      1.78

FE     F   61      54      1.54

AB     M   46      92      1.84

RF     F   31      81      1.56

> #

> mean(mydata$peso_kg) # calcola la media aritmetica dei valori di peso 

[1] 72.85714

> mean(mydata$altezza_m) # calcola la media aritmetica dei valori di altezza 

[1] NA

Come potete vedere il calcolo della media effettuato con la funzione mean() ha avuto successo solamente per il peso. Per l'altezza in luogo del valore della media è stata riportata la sigla NA, evidente conseguenza del valore mancante dell'altezza nel caso GF.

Copiate e incollate nella Console di R queste due ulteriori righe di codice e premete ↵ Invio:

# la corretta gestione dei dati mancanti consente di calcolare la media dei valori di altezza 

#

mean(mydata$peso_kg) # ricalcola la media aritmetica dei valori di peso

mean(mydata$altezza_m, na.rm=TRUE) # ricalcola la media aritmetica dei valori di altezza

#

Il risultato dei calcoli della media che compare nella Console di R ora è questo:

> mean(mydata$peso_kg) # ricalcola la media aritmetica dei valori di peso

[1] 72.85714

> mean(mydata$altezza_m, na.rm=TRUE) # ricalcola la media aritmetica dei valori di altezza

[1] 1.683333

Quindi con l'aggiunta nella funzione mean() dell'argomento na.rm=TRUE, che rimuove dal calcolo i dati mancanti (in questo caso uno solo, ma potrebbero essere anche molti), è stato reso possibile il calcolo della media anche per l'altezza.

Con la funzione na.omit() è infine possibile eliminare definitivamente i casi con dati mancanti: anche un solo dato mancante/campo vuoto comporta la cancellazione dell'intera riga/caso. Copiate e incollate nella Console di R questo script e premete ↵ Invio:

# ELIMINA I CASI CON DATI MANCANTI 

#

# importa i dati, notare / invece di \ su windows

mydata <- read.table("C:/Rdati/importa_csv.csv", header=TRUE, sep=";", dec=",", row.names="id")

#

mydata # mostra i dati importati

#

newdata <- na.omit(mydata) # elimina da mydata i casi con dati mancanti

newdata # mostra i dati dopo eliminazione dei casi con dati mancanti

#

mean(newdata$peso_kg) # calcola la media aritmetica dei valori di peso

mean(newdata$altezza_m) # calcola la media aritmetica dei valori di altezza

#

Dai dati importati nella tabella mydata mediante la funzione na.omit() viene generata una nuova tabella denominata newdata dalla quale sono esclusi i casi con dati mancanti. Con l'eliminazione dalla tabella del caso GF i calcoli vanno subito a buon fine sia per il peso sia per l'altezza.

> mydata # mostra i dati importati 

   sesso anni peso_kg altezza_m

MT     M   69      76      1.78

GF     F   56      63        NA

MC     F   53      71      1.60

SB     M   28      73      1.78

FE     F   61      54      1.54

AB     M   46      92      1.84

RF     F   31      81      1.56

> #

> newdata <- na.omit(mydata) # elimina da mydata i casi con dati mancanti

> newdata # mostra i dati dopo eliminazione dei casi con dati mancanti

   sesso anni peso_kg altezza_m

MT     M   69      76      1.78

MC     F   53      71      1.60

SB     M   28      73      1.78

FE     F   61      54      1.54

AB     M   46      92      1.84

RF     F   31      81      1.56

> #  

> mean(newdata$peso_kg) # calcola la media aritmetica dei valori di peso 

[1] 74.5

> mean(newdata$altezza_m) # calcola la media aritmetica dei valori di altezza 

[1] 1.683333

Da notare che impiegando la funzione na.omit() la media dei valori di peso, che con lo script precedente era 72.85714 kg, ora è cambiata ed è diventata 74.5 kg in quanto in seguito all'eliminazione dell'intero caso GF è stato eliminato anche il suo valore di peso 63 kg che contribuiva a determinare la media di quest'ultimo.

----------

[1] Vedere il post Ordinamento dei dati.

[2] 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

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 quantili parametrici e quantili non parametrici 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 Hmisc, pastecs, psych.

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 proprio,  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'.

https://cran.stat.unipd.it/web/packages/DAAG/DAAG.pdf

[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.

https://www.jstor.org/stable/2291267?seq=1#page_scan_tab_contents

[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.

https://cran.r-project.org/doc/manuals/r-release/fullrefman.pdf

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

[7] Available CRAN Packages By Name. 

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