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.