Effettuare una interpolazione lineare

Data una serie di valori su un piano cartesiano, il valore della variabile dipendente y corrispondente ad uno specifico valore della variabile indipendente x può essere ricavato mediante:

→ interpolazione quando può essere calcolato a partire dalla sua posizione tra due dati adiacenti;

→ approssimazione quando può essere calcolato dalla intersezione con una funzione che approssima i dati e ne fornisce una appropriata descrizione;

→ estrapolazione quando, nell'uno e nell'altro caso, si assume che possa essere calcolato anche al di fuori dell'intervallo dei dati osservati.

L'approssimazione è ampiamente utilizzata in statistica, come ad esempio nel caso del metodo dei minimi quadrati, che consente di calcolare la funzione (retta o polinomio) che, sotto una serie di ipotesi [1], "meglio" approssima i dati. In ogni caso la x della quale si cerca la y corrispondente deve cadere nell'ambito (range) dei valori osservati perché il calcolo mediante la funzione approssimante possa essere ritenuto affidabile.

L'estrapolazione cioè l'estensione delle conclusioni al di la dei valori osservati è un'operazione da proscrivere, a meno che si disponga di un modello ben consolidato. Così ad esempio le leggi della meccanica celeste applicate a una sequenza di posizioni recenti di un asteroide possono consentire una estrapolazione della sua posizione in un tempo futuro, ma i casi di leggi deterministiche di questo genere sono limitati e collegati a campi specifici molto particolari.

Qui vediamo il caso più semplice, la interpolazione lineare, che assume che i due punti adiacenti alla x della quale si cerca la y corrispondente possano essere collegati da un segmento di retta. Nel caso di una curva, più i due punti adiacenti sono ravvicinati, più corto è il segmento di retta, più accurata è l'interpolazione.

Non sono richiesti pacchetti aggiuntivi, in quanto tutte le funzioni che impieghiamo sono  incluse nell'installazione base di R. Ora copiate e incollate nella Console di R questa prima parte dello script e premete ↵ Invio.

# INTERPOLAZIONE LINEARE

#

x <- c(1.21, 1.32, 1.42, 1.53, 1.64, 1.75, 1.85, 1.96, 2.07, 2.18, 2.28, 2.39, 2.50, 2.61, 2.72, 2.82, 2.93, 3.04, 3.15, 3.25, 3.36, 3.47, 3.58, 3.68, 3.79) # valori in ascisse

y <- c(46.86, 47.08, 47.26, 47.42, 47.48, 47.46, 47.37, 47.22, 46.92, 46.61, 46.29, 45.98, 45.64, 45.29, 44.95, 44.64, 44.38, 44.15, 43.98, 43.92, 43.90, 43.93, 43.97, 44.09, 44.26) # valori in ordinate

#

plot(x, y, type="l") # traccia il grafico della funzione

#

Le prime due righe riportano le coppie di coordinate x,y dei dati osservati. Questo è il grafico risultante, tracciato nella terza riga con la funzione plot() che impiega per tutti gli argomenti i valori di default, specificando solamente che i punti devono essere collegati con una linea ("l") [2].

Ora supponiamo di avere per la x i valori 1.3, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.7 e di volere calcolare mediante interpolazione lineare i valori della y corrispondenti.

Per farlo copiate e incollate nella Console di R queste due altre righe di codice e premete ↵ Invio.

#

interp <- approx(x, y, xout=c(1.3, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.7), ties="ordered") # interpola una serie di valori specificati
#

round(cbind(interp$x, interp$y), digits=2) # mostra i risultati dell'interpolazione

#

L'interpolazione viene effettuata con la funzione approx() impiegando come argomenti:

→ le ascisse x dei punti che descrivono la funzione;

→ le ordinate y dei punti che descrivono la funzione;

→ le ascisse xout dei punti da interpolare; 

→ l'indicazione "ordered" che indica che i dati di partenza sono già ordinati.

I risultati dell'approssimazione sono salvati nell'oggetto interp dal quale sono tratti i valori delle ascisse inseriti (interp$x) e i valori delle ordinate calcolati per ciascuno di essi (interp$y). Questi valori sono combinati in una tabella con la funzione cbind(), sono arrotondati con la funzione round() a due decimali (digits=2) e mostrati come segue:

> round(cbind(interp$x, interp$y), digits=2) # mostra i risultati dell'interpolazione

      [,1]  [,2]

 [1,]  1.3 47.04

 [2,]  1.4 47.22

 [3,]  1.6 47.46

 [4,]  1.8 47.42

 [5,]  2.0 47.11

 [6,]  2.2 46.55

 [7,]  2.4 45.95

 [8,]  2.6 45.32

 [9,]  2.8 44.70

[10,]  3.0 44.23

[11,]  3.2 43.95

[12,]  3.4 43.91

[13,]  3.6 43.99

[14,]  3.7 44.12

Con questa riga di codice potete, per controllo, sovrapporre alla curva originaria i punti interpolati.

#

points(interp$x, interp$y, pch=3, col="red") # riporta sul grafico i punti interpolati
#

Come si vede i punti interpolati cadono tutti sulla curva e indicano una adeguata approssimazione. 

Per interpolare un singolo valore e riportarlo sul grafico il codice da impiegare risulta ovviamente più conciso come in questo esempio.

#

val <- approx(x, y, xout=1.7) # interpola un singolo valore

#

points(val, pch=4, col="blue") # riporta sul grafico il punto interpolato

#

Per visualizzare il risultato numerico di questa ulteriore interpolazione è sufficiente nella Console di R digitare

val 

che mostra il contenuto dell'oggetto, rappresentato dal valore della x inserita e da quello della y risultante:

> val

$x

[1] 1.7

$y

[1] 47.46909

----------

[1] Vedere il post La regressione lineare: assunti e modelli.

[2] Digitare help(plot) nella Console di R per la documentazione della funzione.

Barre dell'errore e limiti di confidenza

Questa volta vediamo come con ggplot2 sia possibile, impiegando una sola riga di codice, ottenere rappresentazioni multiple di una statistica (come ad esempio la media) con le corrispondenti barre dell'errore, che tecnicamente definiscono i limiti di confidenza della statistica.

Nel primo dei due esempi che seguono sottraiamo e sommiamo alla media di una variabile 1.96 volte la sua deviazione standard ottenendo i limiti di confidenza al 95% della distribuzione campionaria.

Copiate e incollate nella Console di R questo script e premete ↵ Invio.

# BARRE DELL'ERRORE E LIMITI DI CONFIDENZA della distribuzione campionaria

#

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

library(dplyr) # carica il pacchetto per i calcoli

library(ggplot2) # carica il pacchetto per la grafica

# 

# calcola media e deviazione standard separatamente per sport e per sesso

# 

ais_stat <- ais %>% group_by(sport, sex) %>% summarise(mean_hg=mean(hg), sd_hg=sd(hg)) %>% ungroup()

#

# crea il grafico ma alcuni dati risultano sovrapposti

#

ggplot(ais_stat, aes(x=sport, y=mean_hg, color=sex)) + geom_point(size=3) + geom_errorbar(aes(ymin=mean_hg-1.96*sd_hg, ymax=mean_hg+1.96*sd_hg), width=0.2) + labs(x="Sport praticato", y="Emoglobina in g/dL", color="Sesso dell'atleta") + theme_minimal() + theme(legend.position="top") 

#

Le prime tre righe servono a caricare i pacchetti necessari, che devono essere preventivamente scaricati dal CRAN, il primo dei quali contiene il set di dati ais impiegato come esempio [1]. 

La quarta riga di codice serve per preparare i dati impiegando l'operatore pipe (%>%) che trasferisce i dati da una funzione alla successiva in questo modo:

→ i dati (ais) sono trasferiti (%>%) alla funzione group_by();

→ la funzione group_by() raggruppa i dati in sottoinsiemi per sport praticato (sport) e per sesso dell'atleta (sex) e il risultato di questa suddivisione viene trasferito (%>%) alla funzione summarise();

→ la funzione summarise() calcola la media mean() e deviazione standard sd() della variabile hg che è la concentrazione nel sangue dell'emoglobina (espressa in g/dL, grammi per decilitro di sangue) dei 202 atleti australiani [1] suddivisi per sport e per sesso, e il risultato viene salvato (<-) nella tabella ais_stat;

→ i risultati sono trasferiti (%>%) alla funzione ungroup() che rimuove dalla tabella ais_stat una serie di valori impiegati provvisoriamente per generarla.

L'ultima riga di codice realizza un grafico:

→ impiegando le statistiche salvate nella variabile ais_stat;

→ riportando in ascisse lo sport praticato (sport), in ordinate la media dell'emoglobina (mean_hg), assegnando colori diversi in base al sesso (sex);

→ riportando un grafico a punti (geom_point) dei dati;

→ aggiungendo con geom_errorbar() le barre che indicano l'errore che nel nostro caso è rappresentato da 1.96 volte la deviazione standard (sd_hg);

→ aggiungendo con labs() le opportune etichette; 

→ impiegando con theme_minimal() una rappresentazione semplice;

→ posizionando infine la legenda in alto (legend.position="top").

Questo è il grafico.

Poiché le barre dell'errore nei due sessi risultano sovrapposte, riducendo la chiarezza del grafico, riprendiamo lo stesso codice con una piccola modifica. Copiate e incollate nella Console di R questo script e premete ↵ Invio.

# BARRE DELL'ERRORE E LIMITI DI CONFIDENZA della distribuzione campionaria

#

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

library(dplyr) # carica il pacchetto per i calcoli

library(ggplot2) # carica il pacchetto per la grafica

# 

# calcola media e deviazione standard separatamente per sport e per sesso

# 

ais_stat <- ais %>% group_by(sport, sex) %>% summarise(mean_hg=mean(hg), sd_hg=sd(hg)) %>% ungroup()

#

# crea il grafico evitando la sovrapposizione dei dati

ggplot(ais_stat, aes(x=sport, y=mean_hg, color=sex)) + geom_point(size=3, position=position_dodge(width=0.5)) + geom_errorbar(aes(ymin=mean_hg-1.96*sd_hg, ymax=mean_hg+1.96*sd_hg), width=0.2, position=position_dodge(width=0.5)) + labs(x="Sport praticato", y="Emoglobina in g/dL", color="Sesso dell'atleta") + theme_minimal() + theme(legend.position="top") 

#

L'unica differenza è che aggiungiamo alla funzione geom_point() e alla funzione geom_errorbar() l'argomento position=position_dodge(width=0.5): questo consente di eliminare la sovrapposizione delle barre dell'errore. 

Si può ora notare meglio che le distribuzioni delle concentrazioni di emoglobina nei due sessi si sovrappongono, ma che le medie appaiono differire in modo abbastanza consistente.

Se digitate ais_stat nella Console di R potete visualizzare le medie (mean_hg) e le deviazioni standard (sd_hg), suddivise per ciascuno sport (sport) e per ciascun sesso (sex), calcolate e impiegate per realizzare i grafici.

> ais_stat

# A tibble: 17 × 4

   sport   sex   mean_hg sd_hg

   <fct>   <fct>   <dbl> <dbl>

 1 B_Ball  f        13.1 0.878

 2 B_Ball  m        15.1 0.922

 3 Field   f        14.6 0.682

 4 Field   m        16.0 0.805

 5 Gym     f        13.6 0.860

 6 Netball f        12.8 0.567

 7 Row     f        14.0 0.740

 8 Row     m        15.4 0.711

 9 Swim    f        13.6 0.583

10 Swim    m        15.5 0.655

11 T_400m  f        13.8 1.04 

12 T_400m  m        15.3 0.824

13 T_Sprnt f        14.2 0.556

14 T_Sprnt m        16.2 1.49 

15 Tennis  f        13.5 1.10 

16 Tennis  m        15.6 1.48 

17 W_Polo  m        15.5 0.718

Completiamo quindi l'analisi dei dati: sottraiamo e sommiamo alla media della variabile 1.96 volte il suo errore standard ottenendo i limiti di confidenza al 95% della media campionaria.

Copiate e incollate nella Console di R questo script e premete ↵ Invio.

# BARRE DELL'ERRORE E LIMITI DI CONFIDENZA della media campionaria

#

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

library(dplyr) # carica il pacchetto per i calcoli

library(ggplot2) # carica il pacchetto per la grafica

library(plotrix) # carica il pacchetto per il calcolo dell'errore standard

# 

# calcola media e deviazione standard separatamente per sport e per sesso

# 

ais_stat <- ais %>% group_by(sport, sex) %>% summarise(mean_hg=mean(hg), es_hg=std.error(hg)) %>% ungroup()

#

# crea il grafico evitando la sovrapposizione dei dati

#

ggplot(ais_stat, aes(x=sport, y=mean_hg, color=sex)) + geom_point(size=3, position=position_dodge(width=0.5)) + geom_errorbar(aes(ymin=mean_hg-1.96*es_hg, ymax=mean_hg+1.96*es_hg), width=0.2, position=position_dodge(width=0.5)) + labs(x="Sport praticato", y="Emoglobina in g/dL", color="Sesso dell'atleta") + theme_minimal() + theme(legend.position="top") 

#

Da notare che deve essere installato anche il pacchetto plotrix che viene impiegato per calcolare l'errore standard della media. Per il resto viene impiegato esattamente lo stesso codice dello script precedente, semplicemente sostituendo la deviazione standard dell'emoglobina con l'errore standard della media (es_hg) calcolato con la funzione std.error(). Questo è il risultato. 

Come si può rilevare dal grafico, a parte il caso del tennis (Tennis), all'interno di ciascuno sport i limiti di confidenza al 95% delle medie dell'emoglobina negli atleti dei due sessi non si sovrappongono.

Se digitate ais_stat nella Console di R potete ora visualizzare la tabella che contiene media (mean_hg) ed errore standard (es_hg) calcolati per ciascuno sport (sport) e per ciascun sesso (sex).

> ais_stat

# A tibble: 17 × 4

   sport   sex   mean_hg es_hg

   <fct>   <fct>   <dbl> <dbl>

 1 B_Ball  f        13.1 0.243

 2 B_Ball  m        15.1 0.266

 3 Field   f        14.6 0.258

 4 Field   m        16.0 0.232

 5 Gym     f        13.6 0.430

 6 Netball f        12.8 0.118

 7 Row     f        14.0 0.158

 8 Row     m        15.4 0.184

 9 Swim    f        13.6 0.194

10 Swim    m        15.5 0.182

11 T_400m  f        13.8 0.312

12 T_400m  m        15.3 0.194

13 T_Sprnt f        14.2 0.278

14 T_Sprnt m        16.2 0.450

15 Tennis  f        13.5 0.414

16 Tennis  m        15.6 0.741

17 W_Polo  m        15.5 0.174

Conclusione: 

→ il termine barra dell'errore o errorbar indica semplicemente l'espediente grafico impiegato per rappresentare attorno al valore di una statistica i limiti di confidenza della statistica al livello di probabilità desiderato. In genere si impiega il livello di probabilità del 95% in modo che la probabilità che il valore della statistica – stimato mediante i dati campionari – possa uscire dai limiti di confidenza così calcolati sia "solamente" del 5%;

→ i limiti di confidenza al 95% della distribuzione campionaria – calcolati come media±1.96 ∙ ds (deviazione standard) – sono i limiti all'interno dei quali cade il 95% dei valori osservati nel campione esaminato, quindi la probabilità che un nuovo dato risulti al di fuori di questi limiti è del 5%;

→ i limiti di confidenza al 95% della media – calcolati come media±1.96 ∙ es (errore standard) – sono i limiti all'interno dei quali abbiamo il 95% di probabilità che cada la media, stante l'informazione contenuta nel nostro campione. Questi limiti forniscono un test di significatività della differenza tra le medie: quando i limiti di confidenza al 95% di due medie non si sovrappongono, possiamo affermare che tra le due medie esiste una differenza statisticamente significativa. Quindi nel nostro caso possiamo concludere che all'interno di ciascuno degli sport praticati da entrambi i sessi, con una sola eccezione, la concentrazione media dell'emoglobina nelle atlete è significativamente inferiore a quella degli atleti e che la probabilità che questa conclusione sia errata è del 5% o meno;

→ il risultato di una statistica dovrebbe essere fornito riportando anche la sua incertezza, espressa in termini dei suoi limiti di confidenza (opportunamente specificati, e.g.: limiti di confidenza al 90%, limiti di confidenza al 95%, media±ds, media±es, eccetera).

----------

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

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.

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

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