La
regressione lineare tout court, quella illustrata in tutti i
testi di statistica e qui sviluppata con R, è una regressione
lineare semplice (in contrapposizione alla
regressione lineare multipla), parametrica (in
contrapposizione alla regressione lineare non parametrica), x variabile indipendente (in contrapposizione ad alternative che non prevedono questo assunti). Dato che la denominazione razionale, quella completa, diverrebbe chilometrica, e dato che l'espressione regressione lineare è però troppo generica, laddove è opportuno semplificare impiegherò la denominazione regressione lineare ordinaria, con riferimento al fatto che è quella non solo dovunque illustrata ma anche più ampiamente impiegata.
Gli
assunti alla base del modello matematico-statistico implicano una
serie di requisiti che devono essere soddisfatti dai dati, che sono
ben chiariti ad esempio in Marubini [1] e in Snedecor [2], ma anche
online [3]. Per il caso speciale, ma non infrequente, nel quale
nessuna delle due variabili confrontate abbia i requisiti richiesti
per essere considerata come variabile indipendente, vedere anche il
post La regressione lineare: assunti e modelli.
Essendo
x la variabile indipendente posta sull'asse delle ascisse, y
la variabile dipendente posta sull'asse delle ordinate, ed essendo
soddisfatti i requisiti previsti per i dati, il metodo dei minimi
quadrati consente di calcolare l'intercetta a e il
coefficiente angolare b della retta di regressione di
equazione
y
= a + b∙x
che
meglio approssima la distribuzione dei dati sperimentali.
In R
il calcolo dell'equazione della retta di regressione viene effettuato
mediante la funzione lm(), che
può essere applicata anche al caso di più variabili indipendenti,
consentendo quindi il calcolo della regressione lineare multipla [4].
Ma la cosa più interessante è che R, oltre al calcolo
dell'equazione della retta di regressione, un calcolo di per sé semplice, fornisce una serie molto interessante di strumenti che
consentono di valutare quanto i dati soddisfano i requisiti
richiesti, o in altre parole di valutare se la regressione lineare
descrive in modo adeguato la relazione tra le due variabili.
Qui
impieghiamo come esempi due diversi set di dati, per uno solo dei
quali, come vedremo, la regressione lineare fornisce risultati
adeguati.
Il
primo è il set di dati ais, nel quale prendiamo in considerazione
la concentrazione degli eritrociti (espressa in 1012/L) e
il valore ematòcrito (espresso in %), due variabili che, analizzate
mediante i coefficienti di correlazione parametrici e non parametrici, hanno
mostrato di essere correlate, e che, all'ispezione visiva di una serie di grafici di dispersione (xy) realizzati mediante i correlogrammi, hanno mostrato valori ben allineati su una possibile retta.
I
pacchetti aggiuntivi DAAG, gvlma e car, se non
l'avete già fatto, devono essere preventivamente scaricati e
installati dal CRAN. Copiate e incollate nella
Console di R questo script e premete
↵ Invio:
#
REGRESSIONE LINEARE SEMPLICE PARAMETRICA y = a + b · x
#
library(DAAG)
# carica il pacchetto che include il set di dati ais
library(gvlma)
# carica il pacchetto per la funzione gvlma()
library(car)
# carica il pacchetto per analisi ouliers e analisi grafica
str(ais)
# mostra la struttura dei dati
#
var_x
<- ais$rcc # variabile eritrociti x in ascisse
var_y
<- ais$hc # variabile ematòcrito y in ordinate
reglin
<- lm(var_y ~ var_x) # calcola intercetta (a) e
coefficiente angolare (b)
coefficients(reglin)
# mostra i coefficienti dell'equazione hc = a + b · rcc
confint(reglin,
level=0.95) # calcola gli intervalli di confidenza
dell'intercetta e del coefficiente angolare
#
#
analisi statistica della adeguatezza della regressione lineare
#
summary(reglin)
# mostra un riepilogo dei risultati
t.test(residuals(reglin))
# verifica che la media degli errori non sia significativamente
diversa da zero
shapiro.test(residuals(reglin))
# verifica la normalità della distribuzione degli errori
summary(gvlma(reglin))
# test globale per l'assunto di linearità
outlierTest(reglin)
# valore p di Bonferonni per la presenza di dati anomali (outliers)
#
#
analisi grafica della adeguatezza della regressione lineare
#
windows()
# apre una nuova finestra
par(mfrow=c(2,2))
# predispone la suddivisione della finestra in quattro quadranti, uno
per grafico
#
newx
= seq(min(var_x), max(var_x), by = 0.01) # valori della x
per i quali calcolare l'intervallo di confidenza
conf_interval
<- predict(reglin, newdata=data.frame(var_x=newx),
interval="confidence", level = 0.95) # calcola
gli intervalli di confidenza
plot(var_x,
var_y, xlab="Eritrociti (10^12/L)", ylab="Ematòcrito
(%)", main="Regressione lineare y = a + b·x")
# grafico dei dati
abline(reglin,
col="lightblue") # retta di regressione
lines(newx,
conf_interval[,2], col="blue", lty=2) # limite
di confidenza inferiore
lines(newx,
conf_interval[,3], col="blue", lty=2) # limite
di confidenza superiore
#
plot(var_y,
residuals(reglin), xlab="Ematòcrito (%)", ylab="Ematòcrito
osservato - calcolato (%)", main="Analisi delle differenza
residue") # grafico delle differenza tra ematòcrito
osservato e ematòcrito calcolato con l'equazione della retta
#
influencePlot(reglin,
xlab="t-quantili (il diametro dei cerchi", sub="è
proporzionale alla distanza di Cook)", ylab="Residui
studentizzati", main="Influenza dei dati")
# grafico dell'influenza dei dati sulle conclusioni
#
qqPlot(reglin,
xlab="t-quantili", ylab="Residui studentizzati",
main="Quantili vs. residui") # mostra il grafico
dei quantili per i residui studentizzati
#
Dopo
avere caricato i pacchetti aggiuntivi necessari, e dopo avere
mostrato con str(ais) la
struttura dei dati, con var_x <- ais$rcc e con var_y <- ais$hc sono
memorizzate negli oggetti var_x e
var_y rispettivamente la
variabile indipendente x, posta in ascisse, e la variabile
dipendente y, posta in ordinate. In questo modo sarà
possibile riutilizzare per intero lo script semplicemente sostituendo
ad ais$rcct e ad ais$hc
i vettori contenenti i propri dati.
Con
reglin <- lm(var_y ~ var_x)
l'equazione della retta di regressione che esprime la y
in funzione della x
viene calcolata e memorizzata nell'oggetto reglin,
che a questo punto diventa l'argomento chiave, l'argomento contenente
i dati in ingresso impiegati nelle funzioni successive [5].
Con
le funzioni coefficients() e
confint() sono infine mostrati i
coefficienti della retta di regressione e i loro intervalli di
confidenza al 95%
> coefficients(reglin) # mostra i coefficienti dell'equazione hc = a + b · rcc
(Intercept) var_x
8.183033 7.398052
> confint(reglin, level=0.95) # calcola gli intervalli di confidenza dell'intercetta e del coefficiente angolare
2.5 % 97.5 %
(Intercept) 6.173717 10.192350
var_x 6.974206 7.821898
e
pertanto questa risulta essere l'equazione della retta di regressione:
hc =
8.183033
+ 7.398052
· rcc
A
questo punto seguono due blocchi di codice, il primo per effettuare
una analisi statistica, e il secondo per effettuare una analisi
grafica della regressione, entrambe allo scopo di valutare, come già
detto, se la regressione lineare descrive in modo adeguato la
relazione tra le due variabili.
Per
quanto concerne l'analisi statistica, le principali conclusioni sono
quelle tratte con il test globale per l'assunto di linearità
effettuato mediante la funzione gvlma(),
che conferma il fatto che gli assunti che stanno alla base del
modello di regressione lineare sono tutti soddisfatti:
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = reglin)
Value p-value Decision
Global Stat 3.479904 0.4809 Assumptions acceptable.
Skewness 2.233166 0.1351 Assumptions acceptable.
Kurtosis 0.958154 0.3277 Assumptions acceptable.
Link Function 0.005239 0.9423 Assumptions acceptable.
Heteroscedasticity 0.283345 0.5945 Assumptions acceptable.
Il
test di Bonferroni evidenzia la presenza di un possibile dato anomalo
(dato aberrante o outlier), il dato numero
68:
> outlierTest(reglin) # valore p di Bonferonni per la presenza di dati anomali (outliers)
No Studentized residuals with Bonferonni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferonni p
68 -3.569525 0.00044828 0.090552
Le
conclusioni dell'analisi grafica confermano, con i valori delle
differenze residue dispersi in modo casuale, che la regressione
lineare descrive in modo adeguato la relazione tra le due variabili
fare click per ingrandire
e
confermano anche la possibile presenza di dati anomali con il
grafico delle distanze di Cook e con il grafico dei quantili, dati
che sono evidenziati sia direttamente nei grafici, sia nei rispettivi
riepiloghi riportati nella
Console di R:
> influencePlot(reglin, xlab="t-quantili (il diametro dei cerchi", sub="è proporzionale alla distanza di Cook)", ylab="Residui studentizzati", main="Influenza dei dati") # grafico dell'influenza dei dati sulle conclusioni
StudRes Hat CookD
68 -3.569525 0.006301040 0.03815682
78 3.186317 0.009724284 0.04766681
161 -2.179789 0.039758811 0.09655642
166 1.363972 0.099962743 0.10287127
> #
> qqPlot(reglin, xlab="t-quantili", ylab="Residui studentizzati", main="Quantili vs. residui") # mostra il grafico dei quantili per i residui studentizzati
[1] 68 78
I
dati per i quali i metodi di identificazione meglio concordano sono
il numero
68 e il numero
78, ma anche i dati numero
161 e numero
166 meritano di essere valutati. Questi dati
corrispondono in effetti a quelli che più attirano l'attenzione in
un grafico di dispersione impiegato per illustrare la funzione di R
che consente di identificare i dati in un grafico con un click. E
stanno a indicare i casi che
influenzano in modo importante la regressione. Si tratta di casi dei
quali sarebbe importante controllare la validità, o che potrebbero
essere situati in intervalli di valori per i quali potrebbe essere
opportuno acquisire più dati.
Vediamo
ora la stessa identica analisi applicata al set di dati galton. I
pacchetti aggiuntivi psych, gvlma e car, se non
l'avete già fatto, devono essere preventivamente scaricati e
installati dal CRAN. Copiate e incollate nella
Console di R questo script e premete
↵ Invio:
#
REGRESSIONE LINEARE SEMPLICE PARAMETRICA y = a + b · x
#
library(psych)
# carica il pacchetto che include il set di dati galton
library(gvlma)
# carica il pacchetto per la funzione gvlma()
library(car)
# carica il pacchetto per analisi outliers e analisi grafica
str(galton)
# mostra la struttura dei dati
#
var_x
<- galton$parent # variabile altezza dei padri x in
ascisse
var_y
<- galton$child # variabile altezza dei figli y in
ordinate
reglin
<- lm(var_y ~ var_x) # calcola intercetta (a) e
coefficiente angolare (b)
coefficients(reglin)
# mostra i coefficienti dell'equazione hc = a + b · rcc
confint(reglin,
level=0.95) # calcola gli intervalli di confidenza
dell'intercetta e del coefficiente angolare
#
#
analisi statistica della adeguatezza della regressione lineare
#
summary(reglin)
# mostra un riepilogo dei risultati
t.test(residuals(reglin))
# verifica che la media degli errori non sia significativamente
diversa da zero
shapiro.test(residuals(reglin))
# verifica la normalità della distribuzione degli errori
summary(gvlma(reglin))
# test globale per l'assunto di linearità
outlierTest(reglin)
# valore p di Bonferonni per la presenza di dati anomali (outliers)
#
#
analisi grafica della adeguatezza della regressione lineare
#
windows()
# apre una nuova finestra
par(mfrow=c(2,2))
# predispone la suddivisione della finestra in quattro quadranti, uno
per grafico
#
newx
= seq(min(var_x), max(var_x), by = 0.01) # valori della x
per i quali calcolare l'intervallo di confidenza
conf_interval
<- predict(reglin, newdata=data.frame(var_x=newx),
interval="confidence", level = 0.95) # calcola
gli intervalli di confidenza
plot(var_x,
var_y, xlab="Altezza dei padri (pollici)", ylab="Altezza
dei figli (pollici)", main="Regressione lineare y = a +
b·x") # grafico dei dati
abline(reglin,
col="lightblue") # retta di regressione
lines(newx,
conf_interval[,2], col="blue", lty=2) # limite
di confidenza inferiore
lines(newx,
conf_interval[,3], col="blue", lty=2) # limite
di confidenza superiore
#
plot(var_y,
residuals(reglin), xlab="Altezza dei figli (pollici)",
ylab="Altezza osservata - calcolata (pollici)",
main="Analisi delle differenza residue") #
grafico delle differenza tra altezza osservata e altezza calcolata
con l'equazione della retta
#
influencePlot(reglin,
xlab="t-quantili (il diametro dei cerchi", sub="è
proporzionale alla distanza di Cook)", ylab="Residui
studentizzati", main="Influenza dei dati")
# grafico dell'influenza dei dati sulle conclusioni
#
qqPlot(reglin,
xlab="t-quantili", ylab="Residui studentizzati",
main="Quantili vs. residui") # mostra il grafico
dei quantili per i residui studentizzati
#
Anche
in questo caso, per quanto concerne l'analisi statistica, le
principali conclusioni sono quelle tratte con il test globale per
l'assunto di linearità effettuato mediante la funzione gvlma(),
che però questa volta ci dice che gli assunti che stanno alla base
del modello di regressione lineare, a parte la curtosi, non sono per
niente soddisfatti:
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = reglin)
Value p-value Decision
Global Stat 25.489 4.011e-05 Assumptions NOT satisfied!
Skewness 8.979 2.731e-03 Assumptions NOT satisfied!
Kurtosis 1.965 1.610e-01 Assumptions acceptable.
Link Function 4.632 3.138e-02 Assumptions NOT satisfied!
Heteroscedasticity 9.913 1.641e-03 Assumptions NOT satisfied!
Il
test di Bonferroni evidenzia la presenza di un possibile dato anomalo
(dato aberrante o outlier), il dato numero
1:
> outlierTest(reglin) # valore p di Bonferonni per la presenza di dati anomali (outliers)
No Studentized residuals with Bonferonni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferonni p
1 -3.512671 0.00046509 0.43161
Le
conclusioni dell'analisi grafica confermano, con i valori delle
differenze residue che mostrano una proporzionalità diretta con l'altezza, che i dati in questione violano uno degli assunti cruciali del modello di regressione lineare
fare click per ingrandire
e
confermano anche la possibile presenza di dati anomali con il grafico
delle distanze di Cook e con il grafico dei quantili, dati che sono evidenziati sia direttamente nei grafici, sia nei rispettivi
riepiloghi riportati nella
Console di R:
> influencePlot(lm(galton$child ~ galton$parent), xlab="t-quantili (il diametro dei cerchi", sub="è proporzionale alla distanza di Cook)", ylab="Residui studentizzati", main="Influenza dei dati") # grafico dell'influenza dei dati sulle conclusioni
StudRes Hat CookD
1 -3.5126706 0.002699826 0.016499436
2 -2.9226403 0.001090010 0.004622768
857 0.4839886 0.008511029 0.001006222
897 2.6611087 0.003740530 0.013207269
898 0.9327550 0.008511029 0.003734739
> #
> qqPlot(lm(galton$child ~ galton$parent), xlab="t-quantili", ylab="Residui studentizzati", main="Quantili vs. residui") # mostra il grafico dei quantili per i residui studentizzati
[1] 1 2
I
dati numero
1,
2, 857, 897, 898
forniti
dalla funzione influencePlot()
stanno di nuovo a indicare i casi che influenzano in modo importante
la regressione, casi dei
quali sarebbe importante controllare la validità, o che potrebbero
essere situati in intervalli di valori per i quali potrebbe essere
opportuno acquisire più dati.
Anche la funzione qqPlot()
fornisce, oltre al grafico, l'indicazione di due punti da
controllare, che sono punti
1,
2 già
indicati dalla funzione precedente.
Il
set di dati galton, oltre a non soddisfare i requisiti fondamentali
che si richiedono ai dati per l'applicazione della regressione, è
anche un esempio di dati nei quali non è chiaro, a dispetto delle
conclusioni tratte da Francis Galton di una "regressione verso la media" delle altezze dei figli rispetto a quelle dei padri, quale delle due variabili debba essere
considerata la variabile indipendente. Le implicazioni di
questo fatto, le conseguenze che esso determina nelle conclusioni
tratte dalla regressione lineare, e il calcolo della regressione
lineare con modelli alternativi, sono discussi a parte nel post La regressione lineare: assunti e modelli.
Si rammentano infine, con queste cinque righe di codice,
fitted(reglin)
# ricalcola i valori mediante l'equazione della retta di regressione
residuals(reglin)
# calcola le differenze residue tra valore osservato e valore
calcolato
anova(reglin)
# analisi della varianza per le differenze spiegate dalle x
vcov(reglin)
# matrice di covarianza dell'intercetta e del coefficiente angolare
influence(reglin)
# influenza dei singoli dati sulla regressione
altrettante
funzioni di R per l'analisi dei dati della regressione. Sono
molto analitiche, con un estremo dettaglio numerico, e sebbene siano
state qui aggiunte principalmente per completezza, potrebbero essere
utili in casi particolari, per una analisi ancor più dettagliata
della regressione, a chi volesse approfondire ulteriormente il tema.
Possono essere impiegate, copiandole e incollandole nella
Console di R, in un
qualsiasi punto degli script, ma necessariamente dopo che è stato
generato l'oggetto reglin.
Come potete notare entrambi gli script sono stato realizzati in modo da rendere immediato il loro riutilizzo: è sufficiente assegnare alla variabile x (var_x <-) e alla variabile y (var_y <-) i vostri nuovi dati (e personalizzare opportunamente titoli e legende).
[1]
Bossi A, Cortinovis I, Duca PG, Marubini E. Introduzione
alla statistica medica.
La Nuova Italia Scientifica, Roma, 1994, ISBN 88-430-0284-8. Il
modello statistico nella regressione lineare,
pp-305-308.
[2]
Snedecor GW, Cochran WG. Statistical
Methods.
The Iowa State University Press, 1980, ISBN 0-8138-1560-6. The
matematical model in linear regression,
pp.153-157.
[3]
Regressione lineare. Da Wikipedia, l'enciclopedia libera. URL
consultato il 06/02/2019: https://goo.gl/YDtkZW
[4]
Digitate help(lm) nella
Console di R per la documentazione della funzione lm().
[5]
Digitate help(nomedellafunzione)
nella Console di R per la
documentazione delle funzioni impiegate. Informazioni più esaurienti
sono contenute nei manuali di riferimento dei pacchetti aggiuntivi
qui impiegati, che trovate alla pagina Available CRAN Packages By
Name. URL consultato il 01/02/20198: https://goo.gl/hLC9BB