Visualizzazione post con etichetta Hmisc. Mostra tutti i post
Visualizzazione post con etichetta Hmisc. Mostra tutti i post

mercoledì 23 gennaio 2019

Coefficienti di correlazione parametrici e non parametrici

Come sottolineato altrove, l'assunzione implicita di un rapporto di causa-effetto in seguito ad un coefficiente di correlazione significativo è una delle trappole mentali più ricorrenti e più subdole della statistica: perché non esiste una connessione logica tra correlazione e causazione [1].

Nonostante da questo derivi l'invito ad impiegare la correlazione con parsimonia, con prudenza e con spirito critico nel trarne le conclusioni, risulta difficile ignorare l'esistenza in R di due strumenti.

Il primo sono i classici coefficienti di correlazione e precisamente:
→ il coefficiente di correlazione (lineare) r di Pearson, il classico test parametrico;
→ il coefficiente di correlazione per ranghi ρ (rho) di Spearman (test non parametrico);
→ il coefficiente di correlazione τ (tau) di Kendall (test non parametrico).

Il secondo sono i correlogrammi [2], una rappresentazione un po' particolare ma interessante della correlazione tra variabili, utile come complemento grafico ai coefficienti di correlazione di cui sopra. 

Va ricordato che il più ampiamente usato (e abbondantemente abusato) coefficiente di correlazione r fornisce una misura del grado di correlazione lineare ed è pertanto valido solamente se le due variabili a confronto hanno una relazione di tipo lineare, mentre i coefficienti di correlazione non parametrici forniscono una misura della correlazione valida anche anche nel caso di relazioni non lineari.

In questo primo script viene impiegato il set di dati galton [3], nel quale sono riportate le altezze (in pollici) di 928 coppie padre-figlio. Sono richiesti i pacchetti aggiuntivi psychTools e Hmisc, che devono essere preventivamente scaricati e installati dal CRAN.

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

# COEFFICIENTI DI CORRELAZIONE parametrici e non parametrici
#
library(psychTools) # carica il pacchetto che include il set di dati galton
str(galton) # mostra la struttura dei dati
#
# calcola i coefficienti di correlazione con i metodi di pearson (il classico r), spearman, kendall
#
round(cor(galton, use="complete.obs", method="pearson"), digits=2) # r di Pearson
round(cor(galton, use="complete.obs", method="spearman"), digits=2) # rho di Spearman
round(cor(galton, use="complete.obs", method="kendall"), digits=2) # tau di Kendall
#
# calcola i coefficienti di correlazione con i livelli di significatività
#
library(Hmisc) # carica il pacchetto
rcorr(as.matrix(galton), type="pearson") # type può essere solo “pearson” (il classico r) o “spearman”
rcorr(as.matrix(galton), type="spearman")
#

Questi sono i risultati ottenuti con i tre metodi per il calcolo della correlazione:

> round(cor(galton, use="complete.obs", method="pearson"), digits=2) # r di Pearson 
       parent child
parent   1.00  0.46
child    0.46  1.00
> round(cor(galton, use="complete.obs", method="spearman"), digits=2) # rho di Spearman 
       parent child
parent   1.00  0.43
child    0.43  1.00
> round(cor(galton, use="complete.obs", method="kendall"), digits=2) # tau di Kendall 
       parent child
parent   1.00  0.33
child    0.33  1.00

La funzione cor() viene impiegata per calcolare i coefficienti di correlazione con gli argomenti:
galton che indica i dati da elaborare;
use=”complete.obs” che comporta l'eliminazione automatica dei casi nei quali manchi l'uno o l'altro dato, argomento che qui non servirebbe impiegare, ma che viene riportato per completezza;
method= che di volta in volta riporta i valori "pearson", "spearman", "kendall", per calcolare le rispettive statistiche.

Da notare che la funzione cor() è stata annidata nella funzione round() che con l'argomento digits=2 provvede ad arrotondare a due cifre decimali i risultati: e che i coefficienti di correlazione non parametrici forniscono, come atteso, valori inferiori a quelli del classico coefficiente di correlazione parametrico r di Pearson (sono più conservativi).

Infine lo script impiega la funzione rcorr() del pacchetto Hmisc per calcolare la significatività del coefficiente r di Pearson e del coefficiente ρ (rho) di Spearman. La probabilità p di osservare per caso i valori di r e di rho qui ottenuti è molto bassa e, come potrete vedere se eseguite lo script, viene riportata come uguale a 0 (zero). Quindi entrambi questi coefficienti di correlazione risultano essere significativi.

In questo secondo script sono impiegati i dati contenuti nella tabella ais del pacchetto DAAG - accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [4].

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

# COEFFICIENTI DI CORRELAZIONE parametrici e non parametrici
#
library(DAAG) # carica il pacchetto incluso il set di dati ais
str(ais) # struttura di ais
mydata <- ais[c(1,2,3,4,5,6,7,8,9,10,11)] # salva sole le colonne con le variabili numeriche in mydata
str(mydata) # mostra la struttura di mydata
#
# method può essere pearson (il classico r), spearman, kendall
#
round(cor(mydata, use="complete.obs", method="pearson"), digits=2) # r di Pearson
round(cor(mydata, use="complete.obs", method="spearman"), digits=2) # rho di Spearman
round(cor(mydata, use="complete.obs", method="kendall"), digits=2) # tau di Kendall
#
# calcola i coefficienti di correlazione con i livelli di significatività
#
library(Hmisc) # carica il pacchetto
rcorr(as.matrix(mydata), type="pearson") # type può essere solo “pearson” (il classico r) o “spearman”
rcorr(as.matrix(mydata), type="spearman")
#
# approfondimenti grafici sulla correlazione lineare
#
library(car) # carica il pacchetto
windows() # apre una finestra grafica
scatterplotMatrix(~rcc+wcc+hc+hg+ferr+bmi+ssf+pcBfat+lbm+ht+wt, regLine = list(method=lm, lty=1, lwd=2, col="red"), smooth=FALSE, diagonal=list(method="density", bw="nrd0", adjust=1, kernel="gaussian", na.rm=TRUE), col = "black", main="Matrice di dispersione", data=ais) # traccia il grafico di dispersione xy che incrocia tutte le variabili
#

Nella parte iniziale dello script dopo avere caricato il pacchetto DAAG con la funzione str(ais) viene mostrata la struttura dell'oggetto/set di dati, che contiene 13 variabili, 11 numeriche e 2 non numeriche. Dato che siamo interessati alle sole variabili numeriche, queste, che sono comprese nelle colonne da 1 a 11 della tabella (o dataframe o dataset) ais, le salviamo in un nuovo oggetto denominato mydata che è quello sul quale lavoreremo. Infine per la conferma dell'avvenuta operazione di selezione delle colonne viene mostrata la struttura del nuovo oggetto con str(mydata).

Il resto del codice è praticamente identico allo script precedente. I risultati sono presentati sotto forma di una matrice che contiene i valori dei coefficienti di correlazione che la funzione cor() calcola per tutte le coppie di variabili. Questi valori sono simmetrici rispetto alla diagonale, nella quale è riportato sempre un valore uguale a 1 per i coefficienti di correlazione di ciascuna variabile con se stessa. I valori evidenziati (manualmente dal sottoscritto) sono quelli che nel successivo grafico mostrano una relazione tra le variabili che più si avvicina ad una retta.

> round(cor(mydata, use="complete.obs", method="pearson"), digits=2) # r di Pearson
         rcc  wcc    hc    hg  ferr  bmi   ssf pcBfat   lbm    ht   wt
rcc     1.00 0.15  0.92  0.89  0.25 0.30 -0.40  -0.49  0.55  0.36 0.40
wcc     0.15 1.00  0.15  0.13  0.13 0.18  0.14   0.11  0.10  0.08 0.16
hc      0.92 0.15  1.00  0.95  0.26 0.32 -0.45  -0.53  0.58  0.37 0.42
hg      0.89 0.13  0.95  1.00  0.31 0.38 -0.44  -0.53  0.61  0.35 0.46
ferr    0.25 0.13  0.26  0.31  1.00 0.30 -0.11  -0.18  0.32  0.12 0.27
bmi     0.30 0.18  0.32  0.38  0.30 1.00  0.32   0.19  0.71  0.34 0.85
ssf    -0.40 0.14 -0.45 -0.44 -0.11 0.32  1.00   0.96 -0.21 -0.07 0.15
pcBfat -0.49 0.11 -0.53 -0.53 -0.18 0.19  0.96   1.00 -0.36 -0.19 0.00
lbm     0.55 0.10  0.58  0.61  0.32 0.71 -0.21  -0.36  1.00  0.80 0.93
ht      0.36 0.08  0.37  0.35  0.12 0.34 -0.07  -0.19  0.80  1.00 0.78
wt      0.40 0.16  0.42  0.46  0.27 0.85  0.15   0.00  0.93  0.78 1.00
> round(cor(mydata, use="complete.obs", method="spearman"), digits=2) # rho di Spearman
         rcc  wcc    hc    hg  ferr  bmi   ssf pcBfat   lbm    ht    wt
rcc     1.00 0.17  0.91  0.89  0.25 0.29 -0.40  -0.50  0.59  0.40  0.42
wcc     0.17 1.00  0.17  0.14  0.05 0.23  0.15   0.14  0.11  0.04  0.17
hc      0.91 0.17  1.00  0.95  0.25 0.32 -0.44  -0.54  0.63  0.43  0.45
hg      0.89 0.14  0.95  1.00  0.31 0.37 -0.43  -0.54  0.67  0.41  0.49
ferr    0.25 0.05  0.25  0.31  1.00 0.30 -0.11  -0.19  0.34  0.17  0.30
bmi     0.29 0.23  0.32  0.37  0.30 1.00  0.29   0.15  0.70  0.35  0.84
ssf    -0.40 0.15 -0.44 -0.43 -0.11 0.29  1.00   0.96 -0.22 -0.10  0.13
pcBfat -0.50 0.14 -0.54 -0.54 -0.19 0.15  0.96   1.00 -0.41 -0.25 -0.05
lbm     0.59 0.11  0.63  0.67  0.34 0.70 -0.22  -0.41  1.00  0.80  0.91
ht      0.40 0.04  0.43  0.41  0.17 0.35 -0.10  -0.25  0.80  1.00  0.78
wt      0.42 0.17  0.45  0.49  0.30 0.84  0.13  -0.05  0.91  0.78  1.00
> round(cor(mydata, use="complete.obs", method="kendall"), digits=2) # tau di Kendall
         rcc  wcc    hc    hg  ferr  bmi   ssf pcBfat   lbm    ht    wt
rcc     1.00 0.11  0.75  0.71  0.17 0.19 -0.26  -0.33  0.41  0.27  0.28
wcc     0.11 1.00  0.11  0.09  0.03 0.15  0.10   0.10  0.08  0.03  0.12
hc      0.75 0.11  1.00  0.83  0.17 0.22 -0.29  -0.35  0.44  0.29  0.30
hg      0.71 0.09  0.83  1.00  0.22 0.25 -0.29  -0.36  0.46  0.28  0.33
ferr    0.17 0.03  0.17  0.22  1.00 0.20 -0.07  -0.13  0.22  0.11  0.20
bmi     0.19 0.15  0.22  0.25  0.20 1.00  0.20   0.11  0.52  0.24  0.65
ssf    -0.26 0.10 -0.29 -0.29 -0.07 0.20  1.00   0.82 -0.13 -0.07  0.09
pcBfat -0.33 0.10 -0.35 -0.36 -0.13 0.11  0.82   1.00 -0.24 -0.16 -0.01
lbm     0.41 0.08  0.44  0.46  0.22 0.52 -0.13  -0.24  1.00  0.61  0.77
ht      0.27 0.03  0.29  0.28  0.11 0.24 -0.07  -0.16  0.61  1.00  0.59
wt      0.28 0.12  0.30  0.33  0.20 0.65  0.09  -0.01  0.77  0.59  1.00

L'ultima parte dello script realizza un grafico che fornisce la necessaria integrazione grafica all'analisi numerica.

Impiegando la funzione scatterplotMatrix() viene realizzato con una sola riga di codice un grafico che riporta sulla diagonale per ciascuna variabile il kernek density plot, e per ciascuna coppia di variabili il grafico di dispersione (xy) con la relativa retta di regressione.


I dettagli delle funzioni qui impiegate, tutte molto semplici, li trovate come al solito digitando help(nomedellafunzione) nella Console di R.


----------

[1] Vedere il post Correlazione e causazione.

[2] Vedere il post Correlogrammi.

[3] Vedere il post Il set di dati Galton.

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


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