Statistica e grafica con R

L'obiettivo è riprendere gli aspetti della correlazione che sono stati trattati separatamente per la parte statistica (coefficienti di correlazione parametrici e non parametrici) e per la parte grafica (correlogrammi), aggiungere test di normalità e regressione lineare, quindi riunire il tutto in uno strumento unico di utilizzo pratico che consente di analizzare i dati per ricercare una possibile correlazione.

Lo facciamo con sei script, tre (uno per ciascuno dei coefficienti di correlazione che vedremo) per la ricerca di una possibile correlazione nel caso di dati multivariati e tre (per altrettanti coefficienti di correlazione e due tipi di regressione lineare) per la ricerca di una possibile correlazione tra due variabili, nello spirito che muove tutti i post del blog: costituire e arricchire sempre più una collezione personale di script da tenere a portata di mano per un impiego immediato al bisogno.

Gli script che vediamo consentono per la parte statistica di calcolare:

→ il coefficiente di correlazione lineare r di Pearson (test parametrico)

→ il coefficiente di correlazione per ranghi ρ (rho) di Spearman (test non parametrico)

→ il coefficiente di correlazione τ (tau) di Kendall (test non parametrico)

Quale dei tre coefficienti di correlazione impiegare per trarre le conclusioni dai propri dati deve essere valutato di volta in volta alla luce delle seguenti considerazioni:

1) come suggerisce anche la sua denominazione completa e corretta, il coefficiente di correlazione lineare r di Pearson fornisce una misura accurata del grado di correlazione quando sono soddisfatti i seguenti due requisiti:

→ i dati hanno una distribuzione normale (o, meglio, gaussiana) [1] in quanto si tratta di un test parametrico;

→ tra le due variabili a confronto esiste una relazione lineare;

2) la differenza sostanziale con il coefficiente di correlazione lineare r consiste nel fatto che i due coefficienti di correlazione non parametrici (ρ e τ) forniscono una misura della correlazione valida anche quando i dati non soddisfano i requisiti di normalità (gaussianità) e di linearità previsti per r;

3) i test non parametrici sono più conservativi cioè tendono a fornire valori di p superiori (quindi valori di p meno significativi) rispetto al test parametrico; se si preferisce essere più prudenti nelle conclusioni, si possono applicare i test non parametrici anche quando i dati soddisfano i requisiti di gaussianità e di linearità che consentirebbero l'applicazione del test parametrico. Il coefficiente di correlazione τ di Kendall è il più conservativo.

I primi tre script, uno per ciascuno dei coefficienti di correlazione di cui ci occupiamo, sono applicati al contemporaneo confronto tra tutte le possibili coppie di variabili nel caso di dati multivariati e prevedono:

a) il calcolo dello specifico coefficiente di correlazione;

b) il calcolo della significatività del coefficiente di correlazione (con l'eccezione del test τ di Kendall per il quale nel caso di dati multivariati il codice sarebbe troppo complesso - il problema si risolve ricorrendo di volta in volta al confronto tra le due variabili che interessano);

c) la rappresentazione dei dati sotto forma di correlogrammi, che consentono di effettuare una valutazione grafica della linearità dei dati messi a confronto;

d) i test di gaussianità (normalità) dei dati.

Per eseguire gli script è necessario disporre dei pacchetti corrgram, gclus, Hmisc, mblm, moments [2], se non li avete dovete scaricarli e installarli dal CRAN. Installate infine il pacchetto DAAG che contiene nella tabella ais i dati impiegati come esempio.

Vediamo per primo lo script per il calcolo del coefficiente di correlazione lineare r di Pearson, copiatelo e incollatelo nella Console di R e premete ↵ Invio:

# COEFFICIENTE DI CORRELAZIONE LINEARE r (erre) di Pearson (parametrico)

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

mydata <- ais[c(1:11)] # salva in mydata le colonne contenenti variabili numeriche

# coefficiente di correlazione lineare r di Pearson con valori p di significatività

library(Hmisc) # carica il pacchetto

rcorr(as.matrix(mydata), type="pearson") # mostra i valori di r di Pearson e i valori di p corrispondenti

# grafico xy per tutte le coppie di variabili

library(car) # carica il pacchetto

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="Grafici di dispersione", data=mydata) # traccia i grafici di dispersione xy che incrociano tutte le variabili

# correlogrammi

library(corrgram) # carica il pacchetto

windows() # apre una nuova finestra

corrgram(mydata, cor.method = "pearson", order=FALSE, lower.panel=panel.pts, upper.panel=panel.pie, text.panel=panel.txt, main="Correlogramma (r di Pearson)") # correlogramma con grafici xy e grafici a torta

windows() # apre una nuova finestra

corrgram(mydata, cor.method = "pearson", order=FALSE, lower.panel=panel.pts, upper.panel=panel.conf, text.panel=panel.txt, main="Correlogramma (r di Pearson)") # correlogramma con grafici xy e r di Pearson

windows() # apre una nuova finestra

corrgram(mydata, cor.method="pearson", order=FALSE, lower.panel=panel.pts, upper.panel=panel.ellipse, text.panel=panel.txt, diag.panel=panel.density, main="Correlogramma (r di Pearson)") # correlogramma con grafici xy ed ellissi delle tendenze con curva loess

library(gclus) # carica il pacchetto

windows() # apre una nuova finestra

prs <- cor(mydata, method="pearson") # r di Pearson

mydata.col <- dmat.color(abs(prs)) # applica i colori

mydata.o <- order.single(abs(prs)) # dispone le variabili meglio correlate più vicine alla diagonale

cpairs(mydata, mydata.o, panel.colors=mydata.col, gap=.5, main="Gradiente di correlazione (r di Pearson)") # mostra il grafico

# test di gaussianità (asimmetria e curtosi)

library(moments) # carica il pacchetto

print(paste("numero dei dati =", n <- nrow(mydata)), quote=FALSE)

print(paste("deviazione standard del coefficiente di asimmetria (sa) =", sa <- sqrt(6/n)), quote=FALSE)

print(paste("deviazione standard del coefficiente di curtosi(sc) =", sc <- sqrt(24/n)), quote=FALSE)

ca <- skewness(mydata) # coefficiente di asimmetria (ca)

ca_sa <- round(abs(ca/sa), digits=2) # rapporto (ca/sa)

cc <- kurtosis(mydata)-3 # coefficiente di curtosi (cc)

cc_sc <- round(abs(cc/sc), digits=2) # rapporto (cc/sc)

ascurt <- data.frame(ca, ca_sa, cc, cc_sc) # costruisce la tabella

names(ascurt) <- c("ca", "ca/sa", "cc", "cc/sc") # assegna i nomi alle colonne

ascurt # mostra la tabella

Nel primo blocco di codice viene innanzitutto caricato con library(DAAG) il pacchetto che contiene nella tabella ais i dati utilizzati come esempio. Quindi all'oggetto mydata sono assegnate (<-) le variabili numeriche della tabella ais che si vogliono analizzare indicando le colonne che le contengono. Da notare che qui è possibile impiegare la forma ais[c(1:11)] in quanto le colonne sono contigue: in caso contrario è necessario specificare le colonne una per una nella forma ais[c(1,2,3,4,5,6,8,9,10,11)] nella quale, a titolo di esempio, è stata esclusa la colonna 7.

Potete in alternativa scaricare il file ais.csv che contiene gli stessi dati seguendo le indicazioni fornite alla pagina Dati e sostituire

library(DAAG)

con

ais <- read.table("c:/Rdati/ais.csv", header=TRUE, sep=";", dec=",")

come riportato nel post Il set di dati ais. Dovete fare la stessa cosa per importare i dati da un vostro file: in questo caso sostituite ais con il nome che volete assegnare ai dati, sostituite c:/Rdati/ais.csv con il nome e il percorso completo del file che contiene i vostri dati, al bisogno sostituite sep=";" con il separatore dei campi utilizzato nel file e sostituite dec="," con il separatore dei decimali utilizzato nel file. Infine nella seconda riga indicate le colonne con le variabili da analizzare, in una delle due forme riportate sopra.

Nel secondo blocco di codice, dopo avere caricato il pacchetto Hmisc, viene eseguita l'analisi statistica sull'oggetto mydata impiegando la funzione rcorr() e specificando l'argomento type="pearson": in questo modo sono calcolati e riportati in forma di tabella i risultati per il coefficiente di correlazione lineare r di Pearson e i relativi valori di probabilità p.

Questi sono i valori di r calcolati

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

n= 202

e questi i rispettivi valori di p tabulati (i valori 0.0000 sono da intendere come p < 0.0001) che consentono di stabilirne la significatività. Potete impiegare come valore soglia il classico p = 0.05 ovvero essere un poco più esigenti (e prudenti) e porre la soglia per la significatività di r in corrispondenza del valore p = 0.01.

rcc wcc hc hg ferr bmi ssf pcBfat lbm ht wt

rcc 0.0367 0.0000 0.0000 0.0003 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

wcc 0.0367 0.0294 0.0559 0.0610 0.0118 0.0519 0.1262 0.1460 0.2773 0.0270

hc 0.0000 0.0294 0.0000 0.0002 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

hg 0.0000 0.0559 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

ferr 0.0003 0.0610 0.0002 0.0000 0.0000 0.1252 0.0090 0.0000 0.0805 0.0000

bmi 0.0000 0.0118 0.0000 0.0000 0.0000 0.0000 0.0075 0.0000 0.0000 0.0000

ssf 0.0000 0.0519 0.0000 0.0000 0.1252 0.0000 0.0000 0.0030 0.3136 0.0284

pcBfat 0.0000 0.1262 0.0000 0.0000 0.0090 0.0075 0.0000 0.0000 0.0074 0.9978

lbm 0.0000 0.1460 0.0000 0.0000 0.0000 0.0000 0.0030 0.0000 0.0000 0.0000

ht 0.0000 0.2773 0.0000 0.0000 0.0805 0.0000 0.3136 0.0074 0.0000 0.0000

wt 0.0000 0.0270 0.0000 0.0000 0.0000 0.0000 0.0284 0.9978 0.0000 0.0000

Con il terzo blocco di codice, dopo avere caricato il pacchetto car, con la funzione scatterplotMatrix() sono riportati in un'unica figura i grafici xy (grafici di dispersione) che incrociano tutte le variabili: questo consente una valutazione grafica preliminare dei dati, delle loro distribuzioni e delle loro possibili relazioni.

Sulla diagonale sono riportati i kernel density plot che consentono di valutare la distribuzione dei dati. In ciascun grafico xy è riportata a scopo orientativo la retta di regressione, la cui significatività peraltro deve essere valutata con gli altri strumenti che trovate più avanti. L'argomento ~rcc+wcc+hc+hg+ferr+bmi+ssf+pcBfat+lbm+ht+wt qui specifica tutte le variabili del set di dati ais allo scopo di non perdere nessuna potenziale correlazione. Se avete importato i dati da un vostro file ricordatevi di modificare questo argomento, riportando dopo la tilde (~) l'elenco delle variabili che volete rappresentare e collegando l'una all'altra con il segno +.

Nel quarto blocco di codice, sono realizzati innanzitutto mediante le funzioni del pacchetto corrgram tre correlogrammi che riportano come riferimento in basso a sinistra sotto alla diagonale i grafici xy dei dati (lower.panel=panel.pts).

Il primo correlogramma riporta in alto a destra i grafici a torta (upper.panel=panel.pie) che rappresentano i valori di r (cor.method="pearson") con superficie colorata della torta e intensità del colore proporzionali al valore di r e con le gradazioni di blu e di rosso riferite rispettivamente a valori positivi e a valori negativi di r.

Il secondo correlogramma riporta in alto a destra i valori del coefficiente di correlazione lineare r di Pearson con i limiti di confidenza al 95% (upper.panel=panel.conf). Questi forniscono un test di significatività: il valore di r è significativo se i limiti di confidenza non includono lo 0 (zero), mentre non è significativo se includono lo 0.

Il terzo correlogramma riporta in alto a destra le ellissi (upper.panel=panel.ellipse) i cui assi maggiori delineano la tendenza delle due variabili a variare congiuntamente, con le regressioni loess che aiutano a valutare la linearità della loro relazione. Inoltre riporta sulla diagonale i kernel density plot delle variabili rappresentate (diag.panel=panel.density) già visti nel primo grafico, che ci consentono di evidenziare variabili distribuite in modo spiccatamente asimmetrico, come la variabile ferr, o che hanno una evidente distribuzione bimodale, come la variabile rcc: fatti dei quali si deve tenere conto.

Un quarto e ultimo correlogramma viene quindi realizzato mediante le funzioni del pacchetto gclus. Innanzitutto con la funzione corr() e con l'argomento method="pearson" sono calcolati e salvati (<-) in prs i valori del coefficiente di regressione lineare r di Pearson. Quindi sono predisposti i colori (mydata.cor) per evidenziare le variabili meglio correlate e viene stabilito l'ordine della variabili (mydata.o) in modo che le meglio correlate siano più vicine alla diagonale. Ordine e colore delle variabili sono quindi impiegati nella funzione cpairs() per organizzare i grafici all'interno del correlogramma, con i valori di r più significativi le relazioni lineari più evidenti disposti attorno alla diagonale.

Ai correlogrammi seguono i test di normalità (gaussianità). Viene mostrato il numero dei dati

> print(paste("numero dei dati =", n <- nrow(mydata)), quote=FALSE)

[1] numero dei dati = 202

viene calcolata la deviazione standard sa del coefficiente di asimmetria

> print(paste("deviazione standard del coefficiente di asimmetria (sa) =", sa <- sqrt(6/n)), quote=FALSE)

[1] deviazione standard del coefficiente di asimmetria (sa) = 0.172345496886428

e viene calcolata la deviazione standard sc del coefficiente di curtosi

> print(paste("deviazione standard del coefficiente di curtosi(sc) =", sc <- sqrt(24/n)), quote=FALSE)

[1] deviazione standard del coefficiente di curtosi(sc) = 0.344690993772856

Per ognuna della variabili con la funzione skewness() viene calcolato il coefficiente di asimmetria ca dal quale si ricava il rapporto ca/sa, sempre per ognuna delle variabili con la funzione kurtosis() viene calcolato il coefficiente di curtosi cc dal quale si ricava il rapporto cc/sc. Infine i valori calcolati sono assemblati mediante la funzione data.frame() in una tabella:

ca ca/sa cc cc/sc

rcc 0.4160046 2.41 0.662946166 1.92

wcc 0.8352336 4.85 1.449462762 4.21

hc 0.2752244 1.60 0.891229892 2.59

hg 0.1759351 1.02 0.008064219 0.02

ferr 1.2805839 7.43 1.420173468 4.12

bmi 0.9465155 5.49 2.183474974 6.33

ssf 1.1746685 6.82 1.365135223 3.96

pcBfat 0.7595480 4.41 -0.173299598 0.50

lbm 0.3585088 2.08 -0.240119844 0.70

ht -0.1993028 1.16 0.528116566 1.53

wt 0.2406053 1.40 0.385610737 1.12

Se il rapporto ca/sa è superiore a 2.6 la asimmetria è significativa. Se il rapporto cc/sc è superiore a 2.6 la curtosi è significativa. Questo modo semplificato di effettuare il test di significatività è stato adottato per alleggerire il codice dello script, ma nulla impedisce di approfondire al bisogno lo studio delle variabili cui si è maggiormente interessati con altri test, come quelli contenuti nello script riportato nel post test di normalità (gaussianità).

Integrando i risultati:

→ del coefficiente di correlazione e della sua significatività

→ della valutazione grafica della linearità

→ dei test di normalità (gaussianità) dei dati

abbiamo le basi concrete per decidere se sia opportuno basare le nostre conclusioni sul più classico e tradizionale test r di Pearson, test parametrico basato sull'assunto di gaussianità e di linearità dei dati, o ricorrere in alternativa a un test non parametrico, il coefficiente di correlazione per ranghi ρ (rho) di Spearman oppure il coefficiente di correlazione τ (tau) di Kendall, che fornisce una misura della correlazione valida anche quando i dati non soddisfano i requisiti di normalità (gaussianità) e di linearità necessari per applicare r.

Ad esempio per le variabili hc e hg (per il significato delle variabili vedere il post Il set di dati ais) abbiamo che:

→ il valore del coefficiente di correlazione lineare r di Pearson è elevato ed è significativo;

→ i grafici mostrano una evidente relazione lineare tra le due variabili;

→ i test di normalità (gaussianità) confermano per entrambe le variabili una distribuzione gaussiana;

pertanto possiamo basare le nostre conclusioni in merito alla correlazione sul coefficiente r di Pearson.

Nel caso della correlazione tra le variabili bmi e ssf abbiamo che:

→ i valori del coefficiente di correlazione ρ (rho) di Spearman e del coefficiente di correlazione τ (tau) di Kendall sono elevati e sono significativi;

→ i grafici ci mostrano tra le due variabili una relazione praticamente casuale (difficile poterla considerare "lineare");

→ i test di normalità (gaussianità) confermano per entrambe le variabili una distribuzione non gaussiana;

pertanto r non è applicabile e per le nostre conclusioni in merito alla correlazione dobbiamo ricorrere a uno dei due test non parametrici.

Non ho consigli pratici in merito ai criteri da impiegare per la scelta tra l'uno o l'altro dei due test non parametrici, a parte il fatto che rispetto al ρ di Spearman il coefficiente di correlazione τ di Kendall è più robusto (meno influenzato da valori estremi) ed è più conservativo (a parità di condizioni fornisce valori di p maggiori quindi significatività inferiori).

Lo script per l'analisi della correlazione basata sul coefficiente di correlazione per ranghi ρ (rho) di Spearman ricalca il precedente, copiatelo e incollatelo nella Console di R e premete ↵ Invio:

# COEFFICIENTE DI CORRELAZIONE PER RANGHI ρ (rho) di Spearman (non parametrico)

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

mydata <- ais[c(1:11)] # salva in mydata solamente le colonne contenenti variabili numeriche

# coefficiente di correlazione per ranghi ρ di Spearman con valori p di significatività

library(Hmisc) # carica il pacchetto

rcorr(as.matrix(mydata), type="spearman") # mostra i valori di ρ di Sperman e i valori di p corrispondenti

# grafico xy per tutte le coppie di variabili

library(car) # carica il pacchetto

# correlogrammi

library(corrgram) # carica il pacchetto

windows() # apre una nuova finestra

corrgram(mydata, cor.method="spearman", order=FALSE, lower.panel=panel.pts, upper.panel=panel.pie, text.panel=panel.txt, main="Correlogramma (ρ di Spearman)") # correlogramma con grafici xy e grafici a torta

windows() # apre una nuova finestra

corrgram(mydata, cor.method="spearman", order=FALSE, lower.panel=panel.pts, upper.panel=panel.cor, text.panel=panel.txt, main="Correlogramma (ρ di Spearman)") # correlogramma con grafici xy e ρ di Spearman

windows() # apre una nuova finestra

corrgram(mydata, cor.method="spearman", order=FALSE, lower.panel=panel.pts, upper.panel=panel.ellipse, text.panel=panel.txt, diag.panel=panel.density, main="Correlogramma (ρ di Spearman)") # correlogramma con grafici xy ed ellissi delle tendenze con curva loess

library(gclus) # carica il pacchetto

windows() # apre una nuova finestra

prs <- cor(mydata, method="spearman") # ρ di Spearman

mydata.col <- dmat.color(abs(prs)) # applica i colori

mydata.o <- order.single(abs(prs)) # dispone le variabili meglio correlate più vicine alla diagonale

cpairs(mydata, mydata.o, panel.colors=mydata.col, gap=.5, main="Gradiente di correlazione (ρ di Spearman)") # mostra il grafico

# test di gaussianità (asimmetria e curtosi)

library(moments) # carica il pacchetto

print(paste("numero dei dati =", n <- nrow(mydata)), quote=FALSE)

print(paste("deviazione standard del coefficiente di asimmetria (sa) =", sa <- sqrt(6/n)), quote=FALSE)

print(paste("deviazione standard del coefficiente di curtosi(sc) =", sc <- sqrt(24/n)), quote=FALSE)

ca <- skewness(mydata) # coefficiente di asimmetria (ca)

ca_sa <- round(abs(ca/sa), digits=2) # rapporto (ca/sa)

cc <- kurtosis(mydata)-3 # coefficiente di curtosi (cc)

cc_sc <- round(abs(cc/sc), digits=2) # rapporto (cc/sc)

ascurt <- data.frame(ca, ca_sa, cc, cc_sc) # costruisce la tabella

names(ascurt) <- c("ca", "ca/sa", "cc", "cc/sc") # assegna i nomi alle colonne

ascurt # mostra la tabella

Da notare che l'argomento upper.panel=panel.conf che riporta nel correlogramma i limiti di confidenza del coefficiente di correlazione può essere impiegato solamente con il test r di Pearson ed è stato qui sostituito con l'argomento upper.panel=panel.cor

L'analisi della correlazione per dati multivariati si conclude con lo script per l'analisi della correlazione basata sul coefficiente di correlazione τ (tau) di Kendall, che ricalca i due precedenti. Copiatelo e incollatelo nella Console di R e premete ↵ Invio:

# COEFFICIENTE DI CORRELAZIONE PER RANGHI τ (tau) di Kendall (non parametrico)

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

mydata <- ais[c(1:11)] # salva in mydata solamente le colonne contenenti variabili numeriche

# coefficiente di correlazione τ di Kendall

cor(mydata, method="kendall") # mostra i valori di τ di Kendall

# grafico xy per tutte le coppie di variabili

library(car) # carica il pacchetto

# correlogrammi

library(corrgram) # carica il pacchetto

windows() # apre una nuova finestra

corrgram(mydata, cor.method="kendall", order=FALSE, lower.panel=panel.pts, upper.panel=panel.pie, text.panel=panel.txt, main="Correlogramma (τ di Kendall)") # correlogramma con grafici xy e grafici a torta

windows() # apre una nuova finestra

corrgram(mydata, cor.method="kendall", order=FALSE, lower.panel=panel.pts, upper.panel=panel.cor, text.panel=panel.txt, main="Correlogramma (τ di Kendall)") # correlogramma con grafici xy e τ di Kendall

windows() # apre una nuova finestra

corrgram(mydata, cor.method="kendall", order=FALSE, lower.panel=panel.pts, upper.panel=panel.ellipse, text.panel=panel.txt, diag.panel=panel.density, main="Correlogramma (τ di Kendall)") # correlogramma con grafici xy ed ellissi delle tendenze con curva loess

library(gclus) # carica il pacchetto

windows() # apre una nuova finestra

prs <- cor(mydata, method="kendall") # τ di Kendall

mydata.col <- dmat.color(abs(prs)) # applica i colori

mydata.o <- order.single(abs(prs)) # dispone le variabili meglio correlate più vicine alla diagonale

cpairs(mydata, mydata.o, panel.colors=mydata.col, gap=.5, main="Gradiente di correlazione (τ di Kendall)") # mostra il grafico

# test di gaussianità (asimmetria e curtosi)

library(moments) # carica il pacchetto

print(paste("numero dei dati =", n <- nrow(mydata)), quote=FALSE)

print(paste("deviazione standard del coefficiente di asimmetria (sa) =", sa <- sqrt(6/n)), quote=FALSE)

print(paste("deviazione standard del coefficiente di curtosi(sc) =", sc <- sqrt(24/n)), quote=FALSE)

ca <- skewness(mydata) # coefficiente di asimmetria (ca)

ca_sa <- round(abs(ca/sa), digits=2) # rapporto (ca/sa)

cc <- kurtosis(mydata)-3 # coefficiente di curtosi (cc)

cc_sc <- round(abs(cc/sc), digits=2) # rapporto (cc/sc)

ascurt <- data.frame(ca, ca_sa, cc, cc_sc) # costruisce la tabella

names(ascurt) <- c("ca", "ca/sa", "cc", "cc/sc") # assegna i nomi alle colonne

ascurt # mostra la tabella

Anche in questo caso l'argomento upper.panel=panel.conf che riporta nel correlogramma i limiti di confidenza del coefficiente di correlazione e che può essere impiegato solamente con il test r di Pearson viene sostituito con l'argomento upper.panel=panel.cor

Inoltre la funzione rcorr() che tabula anche i limiti di confidenza è impiegata nel test r di Pearson e nel test ρ di Sperman, ma non nel test τ di Kendall, per il quale il calcolo dei limiti di confidenza non è previsto, ed è stata sostituita con la funzione cor()

Il fatto che per i valori di τ di Kendall non siano riportati i valori di p per la significatività (nel caso dei dati multivariati il codice dello script diventerebbe troppo complesso) si risolve semplicemente ricorrendo di volta in volta al confronto tra le due variabili che interessano con uno degli script che seguono.

Vediamo quindi il calcolo del coefficiente di correlazione nel caso del confronto tra due variabili con il primo di tre script che impiegano la stessa sequenza logica descritta per i dati multivariati applicandola al confronto tra due variabili, un fatto che semplifica le cose in quanto:

a) consente il calcolo di tutti e tre i coefficienti di correlazione con un unico script;

b) consente il calcolo della significatività dei coefficienti di correlazione, inclusa quella per il test τ di Kendall;

c) prevede di sostituire i correlogrammi con un semplice grafico xy;

d) rende più ampia la scelta e più semplice l'applicazione dei test per la valutazione della normalità (gaussianità) dei dati.

Il primo script calcola il coefficiente di correlazione lineare r di Pearson e riporta la regressione lineare parametrica. Copiate e incollate nella Console di R lo script e premete ↵ Invio:

# COEFFICIENTE DI CORRELAZIONE nel caso di due variabili

# coefficiente di correlazione lineare r di Pearson e regressione lineare parametrica

x <- c(18, 46, 89, 72, 96, 63, 23, 58, 36, 82, 88) # variabile in ascisse

y <- c(21, 41, 94, 64, 91, 71, 19, 66, 42, 89, 79) # variabile in ordinate

# coefficienti di correlazione

options(scipen=999) # predispone espressione dei numeri in formato fisso

cor.test(x, y, method="pearson") # r di Pearson

cor.test(x, y, method="spearman") # ρ (rho) di Spearman

cor.test(x, y, method="kendall") # τ (tau) di Kendall

# test di gaussianità

library(moments) # carica il pacchetto

agostino.test(x) # test di D'Agostino per il coefficiente di asimmetria di x

anscombe.test(x) # test di Anscombe-Glynn per il coefficiente di curtosi

agostino.test(y) # test di D'Agostino per il coefficiente di asimmetria di y

anscombe.test(y) # test di Anscombe-Glynn per il coefficiente di curtosi

shapiro.test(x) # test di Shapiro-Wilk

shapiro.test(y)

# grafico xy

plot(x, y, xlim = c(min(x), max(x)), ylim = c(min(y), max(y)), main="Regressione lineare (parametrica)")

# aggiunge al grafico la regressione lineare (parametrica) con r di Pearson

reglin <- lm(y~x) # calcola la regressione lineare

a <- reglin$coefficients[1] # intercetta a

b <- reglin$coefficients[2] # coefficiente angolare b

r <- cor(x, y, method="pearson") # coefficiente di correlazione r

lines(c(min(x), max(x)), c(a+b*min(x), a+b*max(x)), col="black", lty=2, lwd=1) # traccia la retta di regressione

legend(min(x), max(y), legend=c(paste("Intercetta a =", round(a, digits=0)), paste("Coefficiente angolare b =", round(b, digits=3)), paste("Equazione: y =", round(a, digits=0), ifelse(sign(b) == 1, "+", "-"), round(b, digits=3), "· x"), paste("r di Pearson =", round(r, digits=3)))) # aggiunge al grafico la legenda

options(scipen=0) # ripristina espressione dei numeri in formato esponenziale

Dopo avere definito nelle prime due righe di codice i dati da analizzare come x e come y, il che consente di sostituirli agevolmente con i propri, e dopo avere predisposto l'espressione dei numeri in formato fisso con options(scipen=999), mediante la funzione cor.test() sono subito calcolati i tre coefficienti di correlazione.

> cor.test(x, y, method="pearson") # r di Pearson

Pearson's product-moment correlation

data: x and y

t = 11.83, df = 9, p-value = 0.0000008695

alternative hypothesis: true correlation is not equal to 0

95 percent confidence interval:

0.8827108 0.9922369

sample estimates:

cor

0.9693169

> cor.test(x, y, method="spearman") # ρ (rho) di Spearman

Spearman's rank correlation rho

data: x and y

S = 14, p-value < 0.00000000000000022

alternative hypothesis: true rho is not equal to 0

sample estimates:

rho

0.9363636

> cor.test(x, y, method="kendall") # τ (tau) di Kendall

Kendall's rank correlation tau

data: x and y

T = 49, p-value = 0.0003334

alternative hypothesis: true tau is not equal to 0

sample estimates:

tau

0.7818182

I coefficienti r di Pearson (p = 0.0000008695), ρ (rho) di Spearman (p < 0.00000000000000022) e τ (tau) di Kendall (p = 0.0003334) sono sempre significativi. Quest'ultimo con suo il valore inferiore di p conferma di essere il più conservativo dei tre.

Il test di asimmetria e il test di curtosi

> agostino.test(x) # test di D'Agostino per il coefficiente di asimmetria

D'Agostino skewness test

data: x

skew = -0.30031, z = -0.55045, p-value = 0.582

alternative hypothesis: data have a skewness

> anscombe.test(x) # test di Anscombe-Glynn per il coefficiente di curtosi

Anscombe-Glynn kurtosis test

data: x

kurt = 1.7415, z = -1.1572, p-value = 0.2472

alternative hypothesis: kurtosis is not equal to 3

confermano una distribuzione normale della variabile x.

Ancora il test di asimmetria e il test di curtosi

> agostino.test(y) # test di D'Agostino per il coefficiente di asimmetria

D'Agostino skewness test

data: y

skew = -0.38458, z = -0.70302, p-value = 0.482

alternative hypothesis: data have a skewness

> anscombe.test(y) # test di Anscombe-Glynn per il coefficiente di curtosi

Anscombe-Glynn kurtosis test

data: y

kurt = 1.7976, z = -1.0283, p-value = 0.3038

alternative hypothesis: kurtosis is not equal to 3

confermano una distribuzione normale della variabile y.

Anche un test globale per la normalità

> shapiro.test(x) # test di Shapiro-Wilk

Shapiro-Wilk normality test

data: x

W = 0.93397, p-value = 0.4524

> shapiro.test(y)

Shapiro-Wilk normality test

data: y

W = 0.91105, p-value = 0.2509

conferma una distribuzione normale delle due variabili.

Integriamo i test statistici con un grafico xy dei dati realizzato mediante la funzione plot(). Da notare che gli argomenti che definiscono limite inferiore e limite superiore dell'asse delle x (xlim) e dell'asse delle y (ylim) sono parametrizzati, nel senso che ricavano valore minimo (min(x) e min(y)) e valore massimo (max(x) e max(y)) impiegati nel grafico dai dati in ingresso, in modo che cambiando i dati le scale degli assi siano automaticamente riadattate. Ma nulla impedisce di sostituire tali valori di volta in volta con valori numerici fissi che possono fornire una miglior rappresentazione del grafico (in questo caso ad esempio sarebbe più appropriato dimensionare entrambi gli assi tra da 0 a 100).

A questo punto sappiamo di avere:

→ due variabili che i test statistici ci indicano avere una distribuzione gaussiana;

→ un grafico che ci indica che tra le due sussiste, nell'ambito dei valori osservati, una relazione lineare.

Pertanto abbiamo quanto ci serve per essere autorizzati a:

→ riportare la correlazione tra le due variabili sotto forma del coefficiente di correlazione lineare r di Pearson;

→ descrivere la relazione di funzione che collega le due variabili mediante l'equazione della retta di regressione.

L'equazione della retta viene calcolata con il tradizionale metodo parametrico (dei minimi quadrati) come lm(y~x) [3] e il coefficiente di correlazione r viene calcolato come cor(x, y, method="pearson"). L'intercetta a e il coefficiente angolare b della retta di regressione sono impiegati per riportarla mediante la funzione lines() nel grafico, che viene infine completato aggiungendo mediante la funzione legend() una legenda che riporta l'equazione della retta e il coefficiente di correlazione r.

Lo script si chiude ripristinando l'espressione dei numeri in formato esponenziale con options(scipen=0).

Vediamo ora cosa accade fare quando i dati non sono distribuiti normalmente (in modo gaussiano). Copiate e incollate nella Console di R lo script, che ricalca il precedente, e premete ↵ Invio:

# COEFFICIENTE DI CORRELAZIONE nel caso di due variabili

# coefficiente di correlazione per ranghi ρ (rho) di Spearman e regressione lineare non parametrica di Siegel

x <- c(18, 16, 19, 12, 17, 13, 15, 58, 64, 36, 82, 88) # variabile in ascisse

y <- c(21, 14, 15, 11, 20, 18, 19, 46, 38, 42, 89, 79) # variabile in ordinate

# coefficienti di correlazione

options(scipen=999) # predispone espressione dei numeri in formato fisso

cor.test(x, y, method="pearson") # r di Pearson

cor.test(x, y, method="spearman") # ρ (rho) di Spearman

cor.test(x, y, method="kendall") # τ (tau) di Kendall

# test di gaussianità

library(moments) # carica il pacchetto

agostino.test(x) # test di D'Agostino per il coefficiente di asimmetria di x

anscombe.test(x) # test di Anscombe-Glynn per il coefficiente di curtosi

agostino.test(y) # test di D'Agostino per il coefficiente di asimmetria di y

anscombe.test(y) # test di Anscombe-Glynn per il coefficiente di curtosi

shapiro.test(x) # test di Shapiro-Wilk

shapiro.test(y)

# grafico xy

plot(x, y, xlim = c(min(x), max(x)), ylim = c(min(y), max(y)), main="Regressione lineare di Siegel (non parametrica)")

# aggiunge al grafico la regressione lineare non parametrica con ρ (rho) di Spearman

library(mblm) # carica il pacchetto per la regressione di Siegel

reglin <- mblm(y ~ x, repeated=TRUE) # calcola la regressione lineare non parametrica di Siegel

a <- reglin$coefficients[1] # intercetta a

b <- reglin$coefficients[2] # coefficiente angolare b

rho <- cor(x, y, method="spearman") # rho di Spearman

lines(c(min(x), max(x)), c(a+b*min(x), a+b*max(x)), col="black", lty=2, lwd=1) # traccia la retta di regressione

legend(min(x), max(y), legend=c(paste("Intercetta a =", round(a, digits=0)), paste("Coefficiente angolare b =", round(b, digits=3)), paste("Equazione: y =", round(a, digits=0), ifelse(sign(b) == 1, "+", "-"), round(b, digits=3), "· x"), paste("rho di Spearman =", round(rho, digits=3)))) # aggiunge al grafico la legenda

options(scipen=0) # ripristina espressione dei numeri in formato esponenziale

In questo caso uno dei test indica una significativa non normalità / non gaussianità nei dati (p < 0.05):

> shapiro.test(x) # test di Shapiro-Wilk

Shapiro-Wilk normality test

data: x

W = 0.79299, p-value = 0.007807

> shapiro.test(y)

Shapiro-Wilk normality test

data: y

W = 0.80417, p-value = 0.01045

A questo punto sappiamo di avere:

→ due variabili che non sono distribuite in modo gaussiano;

→ un grafico che ci indica che tra le due sussiste, nell'ambito dei valori osservati, una relazione lineare.

Le indicazioni che dobbiamo seguire ora sono:

→ esprimere la correlazione mediante un coefficiente di correlazione non parametrico - in questo caso impieghiamo il coefficiente di correlazione per ranghi ρ (rho) di Spearman;

→ riportare una regressione lineare non parametrica - in questo caso impieghiamo la regressione lineare non parametrica di Siegel [4].

Questo è il grafico generato per presentare in forma sintetica i risultati:

Se preferite in alternativa potete:

→ esprimere la correlazione mediante il coefficiente di correlazione τ (tau) di Kendall;

→ riportare (anche in questo caso) la regressione lineare non parametrica di Siegel [5].

Per farlo copiate e incollate nella Console di R quest'ultimo script e premete ↵ Invio:

# COEFFICIENTE DI CORRELAZIONE nel caso di due variabili

# coefficiente di correlazione τ (tau) di Kendall e regressione lineare non parametrica di Siegel

x <- c(18, 16, 19, 12, 17, 13, 15, 58, 64, 36, 82, 88) # variabile in ascisse

y <- c(21, 14, 15, 11, 20, 18, 19, 46, 38, 42, 89, 79) # variabile in ordinate

# coefficienti di correlazione

options(scipen=999) # predispone espressione dei numeri in formato fisso

cor.test(x, y, method="pearson") # r di Pearson

cor.test(x, y, method="spearman") # ρ (rho) di Spearman

cor.test(x, y, method="kendall") # τ (tau) di Kendall

# test di gaussianità

library(moments) # carica il pacchetto

agostino.test(x) # test di D'Agostino per il coefficiente di asimmetria di x

anscombe.test(x) # test di Anscombe-Glynn per il coefficiente di curtosi

agostino.test(y) # test di D'Agostino per il coefficiente di asimmetria di y

anscombe.test(y) # test di Anscombe-Glynn per il coefficiente di curtosi

shapiro.test(x) # test di Shapiro-Wilk

shapiro.test(y)

# grafico xy

plot(x, y, xlim = c(min(x), max(x)), ylim = c(min(y), max(y)), main="Regressione lineare di Siegel (non parametrica)")

# aggiunge al grafico la regressione lineare non parametrica e τ (tau) di Kendall

library(mblm) # carica il pacchetto per la regressione di Siegel

reglin <- mblm(y ~ x, repeated=TRUE) # calcola la regressione lineare non parametrica di Siegel

a <- reglin$coefficients[1] # intercetta a

b <- reglin$coefficients[2] # coefficiente angolare b

rho <- cor(x, y, method="kendall") # tau di Kendall

lines(c(min(x), max(x)), c(a+b*min(x), a+b*max(x)), col="black", lty=2, lwd=1) # traccia la retta di regressione

legend(min(x), max(y), legend=c(paste("Intercetta a =", round(a, digits=0)), paste("Coefficiente angolare b =", round(b, digits=3)), paste("Equazione: y =", round(a, digits=0), ifelse(sign(b) == 1, "+", "-"), round(b, digits=3), "· x"), paste("tau di Kendall =", round(rho, digits=3)))) # aggiunge al grafico la legenda

options(scipen=0) # ripristina espressione dei numeri in formato esponenziale

Questo è il grafico che illustra in forma sintetica i risultati:

La conclusione?

Non è vero che ricercare una correlazione equivale semplicemente a calcolare il coefficiente di correlazione r e la sua significatività.

Nel confronto che abbiamo fatto tra le variabili rcc e wcc

abbiamo un esempio tipico di una situazione nella quale un r significativo (p = 0.0367) è fuorviante e sarebbe sbagliato basarsi solamente su di esso senza considerare il contesto nel quale è stato generato, infatti:

→ una delle due condizioni per applicare il coefficiente di correlazione lineare r di Pearson è che i dati siano distribuiti in modo gaussiano mentre i test mostrano che la variabile rcc ha una asimmetria significativa e la variabile wcc ha asimmeria e curtosi significative, quindi i dati non sono distribuiti in modo gaussiano;

→ l'altra delle due condizioni per applicare il coefficiente di correlazione lineare r di Pearson è che i dati seguano un andamento lineare. Essendo r = 0.15 abbiamo che R² - il coefficiente di determinazione, che rappresenta la variabilità spiegata dalla regressione lineare - è uguale a 0.0225 il che significa che solamente il 2.2% della variabilità osservata è spiegata dalla regressione lineare, troppo poco per assumere che la relazione tra rcc e wcc possa essere rappresentata mediante una retta;

→ il grafico xy di rcc e wcc mostra una distribuzione dei punti in un'area sostanzialmente circolare, nella quale non vi è alcuna tendenza che possa essere ragionevolmente approssimata mediante una retta;

→ queste osservazioni confermano che il coefficiente di correlazione r non è applicabile e indicano la necessità di ricorrere a un coefficiente di correlazione non parametrico - il coefficiente di correlazione per ranghi ρ (rho) di Spearman o il coefficiente di correlazione τ (tau) di Kendall - che è in grado di fornire una misura della correlazione valida anche quando i dati non soddisfano i requisiti di normalità (gaussianità) e di linearità previsti per r.

Ma c'è dell'altro che coinvolge tutti i coefficienti di correlazione in quanto:

→ per valori identici del coefficiente di correlazione, tanto maggiore è il numero di dati sui quali un valore è stato calcolato, tanto maggiore è la sua significatività (cioè tanto minore è il valore di p) - in altre parole è possibile rendere significativo un coefficiente di correlazione aumentando opportunamente il numero dei dati senza che il coefficiente cambi valore;

→ il coefficiente di correlazione misura la forza del legame tra due variabili ponendola tra 0 e 1 e l'esperienza insegna che valori inferiori a 0.5 corrispondono a correlazioni che non hanno rilevanza pratica - qui abbiamo 0.15 per r, 0.17 per ρ, 0.11 per τ, cioè valori ormai prossimi allo 0, in sostanza valori inadeguati;

→ a prescindere dal fatto che r in questo caso non deve essere impiegato per le ragioni riportate sopra, questi valori inadeguati del coefficiente di correlazione risultano paradossalmente significativi sia per r (p = 0.0367), sia per ρ (p = 0.01697), sia per τ (p = 0.02104) in virtù del numero elevato (202) di dati sui quali è stato calcolato;

→ è ragionevole concludere che i coefficienti di correlazione r, ρ e τ in questo caso sono falsamente significativi:

a) a causa della distorsione indotta sui risultati dalla numerosità campionaria;

b) perché p è appena al di sotto della soglia del 5% (p = 0.05) e potremmo optare per una più prudente soglia dell'1% (p = 0.01) che li renderebbe tutti e tre non significativi;

c) perché dal grafico non si riesce proprio a immaginare un qualsivoglia tipo di relazione che colleghi le due variabili;

d) in quanto (soprattutto) sappiamo che nei soggetti sani (come sono gli atleti australiani i cui valori ematologici e biometrici sono riportati nel set di dati ais qui impiegato) i meccanismi di produzione dei globuli rossi (rcc) e dei globuli bianchi (wcc) del sangue sono tra loro indipendenti e quindi non in relazione l'uno con l'altro.

Quest'ultima osservazione ci riporta ai princìpi generali soggiacenti alla correlazione e ci ricorda che:

→ una correlazione in definitiva ha un senso solamente quando sia supportata da un modello di un evento - fisico, chimico-fisico, biologico, fisiologico, econometrico o quant'altro - che si assume spieghi e descriva adeguatamente la relazione di funzione che intercorre tra due variabili;

→ correlazione statistica non significa causazione [6], come ho già sottolineato nel post Correlazione e causazione e come ci ricorda Tyler Vigen nel suo Spurious correlations - l'assunzione implicita di un rapporto di causa-effetto in seguito ad un coefficiente di correlazione significativo è una delle trappole più subdole e più ricorrenti nella statistica.

Nonostante questi limiti della correlazione - anzi, proprio per illustrarli e poterli gestire - mi è parso interessante presentare alcuni strumenti contenuti nei pacchetti aggiuntivi di R e che, impiegati con consapevolezza e spirito critico, possono fornire un valido aiuto quando si vogliono analizzare i dati per ricercare una possibile correlazione.

----------

[1] Cerco per quanto possibile di evitare l'aggettivo "normale" in quanto viene impiegato con troppi differenti significati: la distribuzione cui ci si riferisce è a rigore la "distribuzione gaussiana".

[2] Trovate i relativi manuali in: Available CRAN Packages By Name. URL consultato il 20/04/2022: https://goo.gl/hLC9BB

[3] Per i dettagli sull'impiego della funzione e sulla rappresentazione della retta vedere i post La regressione lineare: assunti e modelli e il post Regressione lineare semplice parametrica.

[4] Due altri metodi per il calcolo della regressione lineare non parametrica, alternativi a quello di Siegel, sono riportati nel post Regressione lineare semplice non parametrica.

[5] In realtà i modelli di regressione qui impiegati, parametrici e non parametrici, prevedono tutti che x sia la variabile indipendente mentre i dati reali che si intende analizzare potrebbero non rispettare questo assunto. Per questo tema, e per i modelli di regressione lineare alternativi che devono essere presi in considerazione, vedere il post La regressione lineare: assunti e modelli.

[5] “Correlation is not causation”. In: Campbell MJ, Machin D. Medical Statistics. A Commonsense Approach. John Wiley & Sons, New York, 1993, ISBN 0-471-93764-9, p. 103.

Statistica e grafica con R

martedì 3 maggio 2022

Esprimere i numeri in formato fisso

mercoledì 20 aprile 2022

Ricercare una possibile correlazione