Regressione polinomiale di terzo grado

Supponiamo di avere, per ciascuno di una serie di valori prefissati x (riportati in ascisse), il risultato di altrettante misure y (riportate in ordinate) e con la funzione plot() visualizziamo i dati sotto forma di un grafico xy.

Copiate il codice che segue nella Console di R e premete ↵ Invio:

#

x <- c(20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300) # valori in ascisse

y <- c(0.695, 0.699, 0.707, 0.722, 0.740, 0.765, 0.790, 0.815, 0.840, 0.865, 0.885, 0.902, 0.920, 0.932, 0.940) # valori in ordinate

plot(y~x, col="black", pch=1, main="Grafico xy", xlab="Valori prefissati (x)", ylab="Risultati delle misure effettuate (y)") # traccia il grafico xy

#

Dalla ispezione del grafico risulta evidente una relazione non lineare.

Questo ci fa escludere la possibilità di descrivere la relazione di funzione che intercorre tra le due variabili mediante una retta. E pone il tema di quale sia la curva che potrebbe meglio descriverla.

La prima cosa da valutare è se esiste un modello deterministico di un evento fisico, chimico-fisico, biologico, econometrico o quant'altro, che si assume descriva adeguatamente la relazione di funzione che intercorre tra le due variabili: in tal caso non resta che applicare l'equazione fornita dal modello. In caso contrario è necessario ricorrere ad un modello empirico: ed è questo che consideriamo essere il caso nostro.

Qualora si debba ricorrere ad un modello empirico, la regola fondamentale è impiegare l'equazione più semplice compatibile con una adeguata descrizione dei dati. Nel nostro caso visto che è evidente dall'ispezione del grafico che l'equazione di una retta 

y = a + b· ​⁠⁠x

calcolata mediante la classica regressione lineare con il metodo dei minimi quadrati sarebbe inadeguata a descrivere la relazione di funzione che lega la y alla x, possiamo pensare di ricorrere ad una equazione appena più complessa, come quella fornita da una regressione polinomiale di secondo grado nella forma

y = a + b· x + c· x²

e calcolata anch'essa con il metodo dei minimi quadrati [1].

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

# REGRESSIONE POLINOMIALE di secondo grado

#

x <- c(20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300) # valori in ascisse

y <- c(0.695, 0.699, 0.707, 0.722, 0.740, 0.765, 0.790, 0.815, 0.840, 0.865, 0.885, 0.902, 0.920, 0.932, 0.940) # valori in ordinate

#

options(scipen=999) # esprime i numeri in formato fisso

#

# calcola la regressione 

#

polisec <- lm(y~poly(x, degree=2, raw=TRUE)) # calcola la regressione

summary(polisec) # riporta le statistiche della regressione

#

# traccia il grafico

#

newx <- seq(min(x), max(x), (max(x)-min(x))/1000) # valori x in corrispondenza dei quali calcolare il valore del polinomio

newy <- predict(polisec, list(x=newx)) # calcola il valore corrispondente del polinomio

plot(y~x, col="black", pch=1, main="Regressione polinomiale di secondo grado", xlab="Valori prefissati (x)", ylab="Risultati delle misure effettuate (y)") # predispone il grafico e riporta i dati

lines(newx, newy, col="black", lty=2, lwd=1) # traccia il polinomio

#

options(scipen=0) # ripristina la notazione scientifica

#

La prima cosa che notiamo è che il polinomio di secondo grado descrive in pratica una retta.

Inoltre se andate alla quinta riga di codice tra le statistiche generate con la funzione summary(polisec) trovate:

Coefficients:

                                     Estimate   Std. Error t value             Pr(>|t|)    

(Intercept)                      0.6573076923 0.0100102431  65.664 < 0.0000000000000002 ***

poly(x, degree = 2, raw = TRUE)1 0.0009554000 0.0001439488   6.637             0.000024 ***

poly(x, degree = 2, raw = TRUE)2 0.0000001299 0.0000004374   0.297                0.772    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Mentre l'intercetta a con un p=0.0000000000000002 e il coefficiente di primo grado b del polinomio con un p=0.000024 sono significativi, il coefficiente di secondo grado c con un p=0.772 non lo è affatto. Pertanto l'equazione del polinomio di secondo grado

y = a + b· x + c· x²  

privata del coefficiente c· x² si riduce all'equazione di una retta

y = a + b · x   

cosa che conferma il grafico ottenuto.

Dati che descrivono una curva con un raggio di curvatura omogeneo, senza massimi, senza minimi e senza punti di flesso - cioè senza punti di passaggio tra due curvature che vanno in senso opposto - sono compatibili con il ramo di una parabola descritta appunto da un polinomio di secondo grado. Ma in questo caso i dati presentano un punto di flesso tra una curvatura iniziale rivolta verso l'alto e una curvatura finale rivolta verso il basso. Il polinomio di secondo grado media tra le due, finendo con il fornire una curvatura talmente ridotta da potere essere assimilata ad una retta.

Dobbiamo pertanto vedere cosa accade applicando il modello immediatamente successivo, rappresentato da una regressione polinomiale di terzo grado nella forma

 y = a + b· x + c· x² + d· x³     

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

# REGRESSIONE POLINOMIALE di terzo grado

#

x <- c(20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300) # valori in ascisse

y <- c(0.695, 0.699, 0.707, 0.722, 0.740, 0.765, 0.790, 0.815, 0.840, 0.865, 0.885, 0.902, 0.920, 0.932, 0.940) # valori in ordinate

#

options(scipen=999) # esprime i numeri in formato fisso

#

# calcola la regressione 

#

politer <- lm(y~poly(x, degree=3, raw=TRUE)) # calcola la regressione

summary(politer) # riporta le statistiche della regressione

#

# traccia il grafico

#

newx <- seq(min(x), max(x), (max(x)-min(x))/1000) # valori x in corrispondenza dei quali calcolare il valore del polinomio

newy <- predict(politer, list(x=newx)) # calcola il valore corrispondente del polinomio

plot(y~x, col="black", pch=1, main="Regressione polinomiale di terzo grado", xlab="Valori prefissati (x)", ylab="Risultati delle misure effettuate (y)") # predispone il grafico e riporta i dati

lines(newx, newy, col="black", lty=2,lwd=1) # traccia il polinomio

#

options(scipen=0) # ripristina la notazione scientifica

#

Le prime due righe di codice servono semplicemente ad assegnare (<-) i valori inclusi nella funzione c() al vettore x e al vettore y che in questo modo conterranno i nostri dati.

Poi mediante la funzione options() con l'argomento scipen=999 facciamo si che i risultati siano espressi in formato fisso, cosa che li renderà meglio leggibili (questo almeno a mio modo di vedere: se volete lasciare i risultati nella notazione scientifica eliminate questa riga di codice o, forse meglio, anteponete a questa riga di codice il simbolo #).

La quarta riga di codice assegna (<-) all'oggetto politer i risultati del calcolo della regressione polinomiale y~poly(x) di terzo grado (degree=3), l'argomento raw=TRUE è impiegato perché non ci interessano i polinomi ortogonali (l'opzione di default è raw=FALSE).

Se siete interessati ad approfondimenti, ricordo che potete visualizzare ad esempio la struttura dell'oggetto politer digitando str(politer) e potete avere informazioni sulle funzioni impiegate nello script digitando help(nomedellafunzione).

Alla quinta riga di codice la funzione summary(politer) consente infine di visualizzare le statistiche generate:

Call:

lm(formula = y ~ poly(x, degree = 3, raw = TRUE))

Residuals:

       Min         1Q     Median         3Q        Max 

-0.0032624 -0.0012957  0.0006167  0.0014060  0.0019740 

Coefficients:

                                         Estimate       Std. Error t value             Pr(>|t|)    

(Intercept)                       0.6966205128205  0.0025292032241 275.431 < 0.0000000000000002 ***

poly(x, degree = 3, raw = TRUE)1 -0.0003180913177  0.0000662417487  -4.802             0.000552 ***

poly(x, degree = 3, raw = TRUE)2  0.0000097653837  0.0000004730826  20.642       0.000000000381 ***

poly(x, degree = 3, raw = TRUE)3 -0.0000000200739  0.0000000009739 -20.612       0.000000000387 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001865 on 11 degrees of freedom

Multiple R-squared:  0.9997,    Adjusted R-squared:  0.9996 

F-statistic: 1.081e+04 on 3 and 11 DF,  p-value: < 0.00000000000000022

La voce Residuals riporta le statistiche delle differenze che, dopo l'approssimazione dei punti con il polinomio di terzo grado, residuano tra i punti stessi e la curva calcolata: differenze piccole come quelle qui osservate indicano un'ottima approssimazione della curva ai dati.

La voce Coefficients riporta dall'alto verso il basso il valore stimato (Estimate) dell'intercetta a (0.6966205128205), del coefficiente di primo grado b (-0.0003180913177), del coefficiente di secondo grado c (0.0000097653837) e del coefficiente di terzo grado d (-0.0000000200739) e ci consente di scrivere l'equazione del polinomio che sarà quindi 

y = 0.6966205128205 - 0.0003180913177· x + 0.0000097653837· x² - 0.0000000200739· x³     

   

Per ciascun coefficiente viene calcolato il valore del test t di Student come rapporto tra il valore stimato (Estimate) e il suo errore standard (Std. Error) e viene infine riportato il valore di probabilità (Pr(>|t|)) di osservare per caso il valore di t riportato: se tale probabilità è sufficientemente bassa ci assumiamo il rischio di considerare significativo il coefficiente in questione. Nel nostro caso essendo i valori di p sempre molto piccoli (di gran lunga inferiori al valore p = 0.05 generalmente considerato come valore soglia per la significatività) possiamo affermare che tutti e tre i coefficienti del polinomio di terzo grado sono significativi e lo rendono un'equazione in grado di descrivere adeguatamente i dati forniti. 

Il coefficiente di correlazione multipla R² (Multiple R-squared: 0.9997) molto elevato  con una statistica F significativa (p-value: < 0.00000000000000022) conferma ulteriormente la bontà dell'approssimazione. E niente meglio del grafico dei punti e del polinomio che li approssima può confermarla.

Il grafico è tracciato in modo da garantire una risoluzione grafica elevata:

→ l'intervallo compreso tra il valore minimo min(x) e il valore massimo max(x) della variabile in ascisse viene suddiviso per ottenere 1000 valori (ma potete cambiarne il numero) salvati nell'oggetto newx e in corrispondenza dei quali calcolare il valore del polinomio;

→ con la funzione predict() sono calcolati e riportati nell'oggetto newy i valori delle ordinate del polinomio corrispondenti ai valori delle ascisse contenuti in newx ;

→ la funzione plot() viene impiegata per riportare sul grafico di dati contenuti in x e in y;

→ la funzione lines() viene impiegata per sovrapporre ai dati il polinomio.

Potete riportare i dati modificando il simbolo impiegato (pch=...) e il suo colore (col="..."), come pure tracciare il polinomio modificando l'aspetto o stile della linea tracciata (lty=...), il suo colore (col="...") e il suo spessore (lwd=...) [2].

Da notare che il grafico è stato tracciato esclusivamente all'interno dei dati, quindi in un'ottica di interpolazione: e questo ci ricorda anche che, trattandosi di un modello empirico, è da escludere l'impiego del polinomio per la estrapolazione, cioè per trarre conclusioni al di fuori dei dati impiegati per calcolarlo. 

A chiudere lo script è la funzione options() che con l'argomento scipen=0 ripristina l'espressione dei numeri nella notazione scientifica [3].

Lo script può essere riutilizzato semplicemente riportando manualmente in x e in y i nuovi dati, o meglio ancora importandoli da un file [4].

----------

[1] Il metodo dei minimi quadrati prevede - ed è quello che abbiamo assunto essere il nostro caso - che la variabile in ascisse (x) abbia valori prefissati, che la variabile in ordinate (y) sia rappresentata da misure (affette da un errore casuale distribuito normalmente) effettuate corrispondenza delle x e minimizza la somma dei quadrati delle differenze y - y' ovvero delle differenze tra ciascun valore misurato y e il valore y' = f(x) corrispondente calcolato mediante la regressione (lineare o polinomiale). Esempi possono essere la risposta y di un sensore misurata in corrispondenza di una serie di concentrazioni date x di una sostanza chimica, la misurazione di un evento y in corrispondenza di una serie data di tempi x, e così via.

[2] Trovate le relative indicazioni:

→ nel post Cambiare i simboli dei punti di R, 

→ nel post Visualizzare i colori disponibili in R 

→ nel post Cambiare gli stili delle linee di R

[3] La notazione scientifica è una notazione esponenziale, cioè una notazione che consente di scrivere un numero N assegnato come prodotto di un opportuno numero a per la potenza di un altro numero b^k essendo b la base della notazione esponenziale. In particolare la notazione scientifica è una notazione esponenziale in base 10 (b = 10). Tuttavia "notazione esponenziale" non è sinonimo di "notazione scientifica" in quanto esistono notazioni esponenziali che impiegano una base diversa da quella impiegata dalla notazione scientifica: ne sono un esempio la notazione esponenziale in base 2 (b = 2) impiegata in informatica, e la notazione esponenziale in base e (b = e) - essendo il numero di Nepero e = 2.718 281 828 459... - impiegata in campo matematico e in campo scientifico.

[Precisazione aggiuntiva: in R il separatore delle cifre decimali è il punto (.) e come già riportato altrove questa convenzione per ragioni di omogeneità viene adottata non solo negli script ma anche nei file di dati e in tutto il testo di questo sito, quindi anche nell'espressione del numero e di Nepero. Per le norme che regolano la punteggiatura all'interno dei numeri rimando al post Come separare i decimali e come raggruppare le cifre].

[4]  Alla pagina Indice nel capitolo R - gestione dei dati trovate i collegamenti ai post che contengono gli script per importare i dati da file esterni.

Ricercare una possibile correlazione (tra due variabili)

La ricerca di una possibile correlazione, intesa come la tendenza di due grandezze a variare congiuntamente in modo statisticamente significativo, richiede – oltre alla necessità ineludibile di spiegare per quale ragione l'una dipenda dall'altra – un approccio integrato che include il calcolo del coefficiente di correlazione, la regressione lineare e l'esecuzione dei test di normalità. Lo vediamo ora con tre script da tenere a portata di mano per un impiego immediato al bisogno. 

Gli script 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 impiegare per trarre dai propri dati le conclusioni in merito alla significatività della correlazione 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) i dati hanno una distribuzione gaussiana [1] in quanto si tratta di un test parametrico; (ii) tra le due variabili a confronto esiste una relazione lineare;

2) i due coefficienti di correlazione non parametrici (ρ e τ) forniscono una misura della correlazione valida anche quando i dati non sono distribuiti in modo gaussiano e tra le due variabili a confronto non esiste una relazione lineare; 

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.

Per eseguire gli script è necessario disporre dei pacchetti mblm e moments [2], se non li avete dovete scaricarli e installarli dal CRAN. Gli script prevedono tutti la stessa sequenza logica che include:

a) il calcolo dei tre coefficienti di correlazione per una immediata comparazione dei loro risultati; 

b) il calcolo della significatività dei coefficienti di correlazione;

c) un grafico xy con la retta di regressione;

d) i test per la valutazione della normalità (gaussianità) dei dati.

Con il primo script viene calcolato il coefficiente di correlazione lineare r di Pearson (test parametrico) e viene rappresentata la regressione lineare parametrica. Copiate e incollate nella Console di R lo script e premete ↵ Invio:

# 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

> 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

 e il test di curtosi 

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

Di nuovo il test di asimmetria

> 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

 e il test di curtosi

> 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 gaussiana 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 che le due variabili sono distribuite in modo gaussiano.

Integriamo i test statistici con un semplice 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 e quindi viene impiegato il coefficiente di correlazione per ranghi ρ (rho) di Spearman (test non parametrico) e viene rappresentata la retta di regressione lineare non parametrica.

Copiate e incollate nella Console di R lo script, che ricalca il precedente, e premete ↵ Invio:

# 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 abbiamo:

→ due variabili che non sono adeguatamente distribuite in modo gaussiano;

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

Le indicazioni che dobbiamo seguire sono pertanto:

→ 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 impiegando gli istessi dati potete esprimere la correlazione mediante il coefficiente di correlazione τ (tau) di Kendall (test non parametrico) e rappresentando sempre la retta di regressione lineare non parametrica [5]. Per farlo copiate e incollate nella Console di R quest'ultimo script e premete ↵ Invio:

# 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

tau <- 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(tau, digits=3)))) # aggiunge al grafico la legenda

#

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

#

La differenza sostanziale rispetto al ρ (rho) di Spearman è un coefficiente di correlazione τ (tau) di Kendall che risulta meno significativo, come previsto per un test più conservativo. 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 lineare r e la sua significatività.

Infatti è necessario ricordare i princìpi generali soggiacenti alla correlazione:

→ uno studio di correlazione ha un senso solamente quando sia supportato 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 le variabili;

→ è indispensabile verificare che gli assunti alla base del modello di correlazione impiegato siano rispettati;

→ esistono alternative al coefficiente di correlazione lineare r che è necessario utilizzare quando detti assunti non sono rispettati;

→ correlazione statisticamente significativa non significa causazione [6], l'assunzione implicita di un rapporto di causa-effetto in seguito ad un coefficiente di correlazione significativo è una delle trappole subdola e ricorrente nella statistica [7].

Nonostante questi limiti, in R sono disponibili gli strumenti che, impiegati in modo appropriato, possono fornire un valido aiuto quando si vogliono analizzare i dati alla ricerca di una possibile correlazione, che in ogni caso va interpretata con attenzione, cautela e spirito critico.  

----------

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

[6] “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.

[7] Vedere quanto ho già sottolineato nel post Correlazione e causazione e quanto ci ricorda Tyler Vigen nel suo Spurious correlations.

Ricercare una possibile correlazione (tra più variabili)

La ricerca di una possibile correlazione, intesa come la tendenza di grandezze a variare congiuntamente in modo statisticamente significativo, richiede – oltre alla necessità ineludibile di spiegare per quale ragione l'una dipenda dall'altra – un approccio integrato che include il calcolo del coefficiente di correlazione, la regressione lineare e l'esecuzione dei test di normalità. Lo vediamo ora con tre script da tenere a portata di mano per un impiego immediato al bisogno. 

Gli script consentono per la parte statistica di calcolare, nel caso del confronto tra più di due variabili:

→ 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 impiegare per trarre dai propri dati le conclusioni in merito alla significatività della correlazione 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) i dati hanno una distribuzione gaussiana [1] in quanto si tratta di un test parametrico; (ii) tra le due variabili a confronto esiste una relazione lineare;

2) i due coefficienti di correlazione non parametrici (ρ e τ) forniscono una misura della correlazione valida anche quando i dati non sono distribuiti in modo gaussiano e tra le due variabili a confronto non esiste una relazione lineare; 

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.

Per eseguire gli script è necessario disporre dei pacchetti car, corrgram, gclus, Hmisc, 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.

Gli script prevedono tutti la stessa sequenza logica che include:

a) il calcolo dello specifico coefficiente di correlazione per tutte le coppie di variabili;

b) il calcolo della significatività dei coefficienti di correlazione;

c) i grafici xy con la retta di regressione per una valutazione grafica della linearità dei dati messi a confronto;

d) la rappresentazione dei dati sotto forma di correlogrammi, per meglio evidenziare le correlazioni;

e) i test per la valutazione della normalità (gaussianità) dei dati.

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. 

Per importare i vostri dati è sufficiente sostituire ais[c(1:11)] con il nome della vostra tabella e rcc, wcc, hc, hg e così via con i nomi delle variabili contenute nella tabella.

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.

P

       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 non dimenticate 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 che trovate nei post.

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

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

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

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

mydata.o <- order.single(abs(rho)) # 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. 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

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

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

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

mydata.o <- order.single(abs(tau)) # 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

#

Lo script ricalca il precedente. 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

Lo sola differenza sostanziale, a parte ovviamente il differente coeffiiente di correlazione, è la sostituzione della funzione rcorr(), che non prevede il calcolo del test τ di Kendall, con la funzione cor(). Il fatto che questa non riporti i valori di p di significatività del τ di Kendall si risolve ricorrendo di volta in volta al confronto diretto tra le due variabili che interessano [3].

La conclusione?

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

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 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 forniscono 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 di correlazione lineare 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.

Ma non è tutto. 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.

C'è infine 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 ancora ai princìpi generali soggiacenti alla correlazione ricordandoci che:

→ uno studio di correlazione ha un senso solamente quando sia supportato 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 le variabili;

→ è indispensabile verificare che gli assunti alla base del modello di correlazione impiegato siano rispettati;

→ esistono alternative al coefficiente di correlazione lineare r che è necessario utilizzare quando detti assunti non sono rispettati;

→ correlazione statisticamente significativa non significa causazione [4], l'assunzione implicita di un rapporto di causa-effetto in seguito ad un coefficiente di correlazione significativo è una trappola subdola e ricorrente nella statistica [5].

Nonostante questi limiti, in R sono disponibili gli strumenti che, impiegati in modo appropriato, possono fornire un valido aiuto quando si vogliono analizzare i dati alla ricerca di una possibile correlazione, che in ogni caso va interpretata con attenzione, cautela e spirito critico. 

----------

[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] Vedere il post Ricercare una possibile correlazione (tra due variabili).

[4] “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.

[5] Vedere quanto ho già sottolineato nel post Correlazione e causazione e quanto ci ricorda Tyler Vigen nel suo Spurious correlations.