Correlazione in R: cor(), cor.test() e matrici di correlazione

Correlazione tra due variabili: cor()

La correlazione si chiede: quando una variabile sale, anche l'altra tende a salire (positiva), a scendere (negativa) o fa per conto suo (vicina a zero)? In R basta una chiamata:

print(cor(mtcars$wt, mtcars$mpg))

La risposta è circa −0.87: le auto più pesanti percorrono meno chilometri con un litro, e la relazione è forte. Quel singolo numero è il coefficiente di correlazione di Pearson, che si scrive ovunque r.

Leggere r

Il coefficiente cade sempre tra −1 e +1. Il segno indica la direzione; il valore assoluto indica la forza:

| |r| | Lettura tipica | | --- | --- | | da 0.0 a 0.2 | trascurabile | | da 0.2 a 0.4 | debole | | da 0.4 a 0.6 | moderata | | da 0.6 a 0.8 | forte | | da 0.8 a 1.0 | molto forte |

Considera queste fasce come spunti di discussione, non come legge: in fisica un r di 0.6 è deludente, in psicologia è il momento clou di una carriera. Due proprietà da interiorizzare: r non ha unità di misura (correlare il peso in tonnellate con mpg dà lo stesso r del peso in chilogrammi con mpg, perché r si calcola su valori standardizzati), e r misura solo l'associazione lineare: una relazione perfetta a forma di U può avere r ≈ 0.

E la frase che va detta: correlazione non significa causalità. Le vendite di gelati e le morti per annegamento sono fortemente correlate lungo i mesi dell'anno, non perché il gelato faccia annegare le persone, ma perché l'estate fa crescere entrambe. Una correlazione ti dice che due variabili si muovono insieme; non dice nulla sul perché. Forse x influenza y, forse y influenza x, forse una terza cosa (la stagione, in questo caso) influenza entrambe. Per scegliere tra queste ipotesi servono esperimenti o un ragionamento causale attento, non un r più grande.

Spearman e Kendall: quando Pearson è lo strumento sbagliato

Pearson lavora sui valori grezzi, e questo lo rende sensibile agli outlier e cieco alle relazioni curve. L'argomento method passa ad alternative basate sui ranghi:

x <- c(1, 2, 3, 4, 5, 6, 7, 8)
y <- c(1, 4, 9, 16, 25, 36, 49, 300)  # increasing, with a wild last value

print(cor(x, y))                        # Pearson - dragged by the outlier
print(cor(x, y, method = "spearman"))   # ranks: perfectly monotonic
print(cor(x, y, method = "kendall"))

Spearman sostituisce ogni valore con il suo rango e poi calcola Pearson sui ranghi. Siccome y cresce sempre quando cresce x, tutti i ranghi si allineano e Spearman riporta esattamente 1: la grandezza dell'outlier smette di contare, conta solo la sua posizione. Scegli Spearman quando i dati sono ordinali (scale di sondaggio), molto asimmetrici, o quando la relazione è monotona ma non rettilinea. Kendall risponde a una domanda simile partendo dalle coppie concordanti e discordanti; è più robusto nei piccoli campioni ma più lento, e Spearman è la scelta predefinita più comune.

La matrice di correlazione

Per analizzare le relazioni tra molte variabili in una volta, passa a cor() diverse colonne numeriche:

m <- cor(mtcars[, c("mpg", "wt", "hp", "disp")])
print(round(m, 2))

Ogni variabile contro ogni altra, con degli 1 sulla diagonale (ogni cosa è perfettamente correlata con sé stessa) e un'immagine speculare ai due lati. Arrotondare a due decimali conta più di quanto sembri: la matrice non arrotondata è un muro di cifre, e lo scopo di una matrice è poterla scorrere con lo sguardo. Qui l'occhiata mostra che mpg è correlata negativamente con tutte e tre (auto più pesanti, più potenti e con motori più grandi consumano più carburante), mentre wt, hp e disp sono tutte fortemente positive tra loro: un gruppo di variabili da "auto grande" che conterà quando arriverai alla regressione lineare e ai suoi grattacapi di multicollinearità.

Valori mancanti: l'argomento use

Con dati mancanti, il comportamento predefinito di cor() è restituire NA invece di tirare a indovinare:

d <- mtcars[, c("mpg", "wt", "hp")]
d$wt[c(3, 7)] <- NA  # knock two values out

print(cor(d)[1, 2])
print(cor(d, use = "complete.obs")[1, 2])
print(cor(d, use = "pairwise.complete.obs")[1, 2])

• use = "complete.obs" elimina ogni riga che contiene qualsiasi NA, poi calcola l'intera matrice con le righe rimaste. Coerente, ma spreca dati: un wt mancante toglie quella riga anche dalla coppia tra mpg e hp.
• use = "pairwise.complete.obs" calcola ogni cella da tutte le righe in cui quella coppia è presente. Conserva più dati, ma celle diverse si basano su sottoinsiemi diversi, il che molto di rado può produrre una matrice non coerente al suo interno.

Per un paio di NA sparsi vanno bene entrambe; basta dire quale hai usato.

È significativa? cor.test()

cor() dà un numero ma nessuna idea se possa essere rumore. cor.test() esegue il test d'ipotesi:

print(cor.test(mtcars$wt, mtcars$mpg))

Leggi l'output blocco per blocco:

• t = −9.56, df = 30: la statistica del test. L'ipotesi nulla è che la correlazione vera sia zero; l'r osservato viene convertito in una statistica t con n − 2 gradi di libertà (32 auto meno 2).
• p-value = 1.29e-10: se la correlazione vera fosse zero, la probabilità di vedere un r così lontano da zero in un campione di 32 sarebbe circa 0.0000000001. È un'evidenza schiacciante che l'associazione sia reale, ma ricorda che il p-value dice se r è diverso da zero, non se la relazione sia grande o causale.
• 95 percent confidence interval: −0.93 to −0.74: l'intervallo plausibile per la correlazione vera. Spesso è più utile del p-value: anche l'estremo più ottimista di questo intervallo è una correlazione negativa forte.
• sample estimates: cor = −0.87: lo stesso numero che ti ha dato cor().

I piccoli campioni meritano qui un rispetto in più: con n = 10, correlazioni di ±0.5 compaiono per caso con una frequenza allarmante, e l'ampio intervallo di confidenza te lo dirà. Riporta l'intervallo, non solo il p-value.

Vederla: fai sempre un grafico

Un coefficiente di correlazione comprime un'intera relazione in un solo numero, e la compressione può nascondere curvature, gruppi, o un unico punto che fa tutto il lavoro. Prima di fidarti di un r, guarda il grafico a dispersione:

plot(mtcars$wt, mtcars$mpg)              # one pair
pairs(mtcars[, c("mpg", "wt", "hp", "disp")])  # every pair in the matrix

pairs() disegna una griglia di grafici a dispersione che corrisponde alla tua matrice di correlazione: il modo più rapido per controllare che i numeri significhino ciò che pensi. Per grafici di matrici curati in stile heatmap, il pacchetto corrplot (install.packages("corrplot"), poi corrplot(cor(m))) è lo strumento standard.

Cosa ti porti a casa

• cor(x, y) dà l'r di Pearson: il segno è la direzione, il valore assoluto la forza, sempre in [−1, 1], senza unità.
• La correlazione misura il movimento congiunto lineare e non dice nulla sulla causalità: una terza variabile nascosta è sempre una candidata.
• method = "spearman" per i ranghi: dati ordinali, outlier, relazioni monotone ma curve.
• cor(df) su colonne numeriche dà la matrice; arrotondala con round(, 2) e fai attenzione all'argomento use = quando mancano dati.
• cor.test(x, y) aggiunge il p-value e un intervallo di confidenza: riporta l'intervallo.
• Guarda sempre il grafico a dispersione prima di credere al numero.

Prossimo passo: quando la domanda si fa più precisa, da "si muovono insieme?" a "la media di questo gruppo è diversa da quella dell'altro?", arriva il t-test.