Regressione lineare multipla

La regressione lineare semplice (regressione bivariata) costruita sul piano cartesiano

y = a + b·x

può essere estesa a uno spazio n-dimensionale: abbiamo così la regressione lineare multipla nella quale la variabile dipendente y viene fatta dipendere non più da una, bensì da k variabili indipendenti, indicate come x₁, x₂, x₃ … xₖ secondo l'equazione

y = a + b₁·x₁ + b₂·x₂ + b₃·x₃ + ... + bₖ·xₖ

Il risultato è rappresentato quindi da una intercetta a e da tanti coefficienti angolari b (indicati con b₁, b₂, b₃ ... bₖ) quante sono le variabili indipendenti (x₁, x₂, x₃ ... xₖ) che contribuiscono a determinare il valore della variabile dipendente y.

Come esempio impieghiamo il set di dati ais nel quale, tra i dati ematologici rilevati in un gruppo di atleti australiani, abbiamo tre grandezze:

→ la concentrazione degli eritrociti (globuli rossi) espressa in milioni per microlitro di sangue (10⁶/µL) e contenuta nella variabile rcc;

→ il valore ematòcrito cioè la frazione del volume del sangue occupata dagli eritrociti espressa in percentuale (%) e contenuta nella variabile hc;

→ la concentrazione dell'emoglobina espressa in grammi per decilitro di sangue (g/dL) e contenuta nella variabile hg.

Le tre grandezze sono strettamente legate tra loro dal punto di vista biologico, in quanto gli eritrociti hanno ciascuno un certo volume e contengono ciascuno una data quantità di emoglobina: pertanto all'aumentare del numero di eritrociti aumenta la frazione percentuale del volume del sangue occupata dagli eritrociti (il valore ematòcrito) e aumenta la quantità di emoglobina. Mentre la diminuzione dell'ematòcrito e la diminuzione dell'emoglobina sono indice della diminuzione del numero degli eritrociti.

L'esempio è stato scomposto in una serie di script indipendenti - onde facilitare il riutilizzo del codice per analizzare i propri dati - che impiegano i pacchetti aggiuntivi gvlma, car e rgl, che devono essere scaricati dal CRAN. Accertatevi di avere installato inoltre il pacchetto DAAG che contiene i dati impiegati nello script, o in alternativa  procedete come indicato nel post Il set di dati ais. 

Iniziamo con una analisi esplorativa preliminare dei dati, copiate questo primo script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - analisi esplorativa preliminare dei dati

#

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

library(car) # carica il pacchetto per la grafica

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

windows() # apre e inizializza una nuova finestra grafica

#

scatterplotMatrix(~rcc+hc+hg, col="black", pch=20, regLine=list(method=lm, lty=1, lwd=2, col="chartreuse3"), smooth=FALSE, diagonal=list(method="histogram", breaks="FD"), main="Istogrammi e grafici xy con rette di regressione", data=mydata) # istogrammi e grafici xy con rette di regressione

#

Dopo avere caricato i pacchetti necessari, copiato nell'oggetto mydata le tre variabili che ci interessano, aperto e inizializzato una nuova finestra grafica, l'analisi preliminare dei dati viene effettuata con la funzione scatterplotMatrix() nella quale:

→ ~rcc+hc+hg specifica le variabili da rappresentare, che sono la concentrazione degli eritrociti rcc, il valore ematòcrito hc e la concentrazone di emoglobia hg;

→ col="..." specifica il colore dei grafici; 

→ pch=... specifica il tipo di simbolo da impiegare (un piccolo cerchio pieno);

→ per la retta di regressione regLine sono impiegati il modello lineare (lm), una linea continua (lty=1) con uno spessore doppio (lwd=2) e un colore (col="chartreuse3"), mentre l'argomento smooth=FALSE fa si che non venga aggiunta ai grafici la rappresentazione della funzione non lineare prevista di default nel pacchetto car;

→ l'argomento diagonal=list() indica la rappresentazione della variabile riportata nella diagonale, che qui è l'istogramma (method="histogram") con il numero di classi (breaks="FD") calcolato mediante la regola di Freedman-Diaconis. 

Da notare che le espressioni che è possibile impiegare per l'argomento  diagonal=... sono:

→ list(method="boxplot");

→ list(method="density", bw="nrd0", adjust=1, kernel="gaussian", na.rm=TRUE);

→ list(method="histogram");

→ list(method="oned");

→ list(method="qqplot").

Vista la relazione lineare molto buona che intercorre tra le variabili sembra ragionevole continuare ed effettuare il calcolo della regressione lineare multipla, copiate questo secondo script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - calcolo della regressione lineare multipla

#  

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

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

reglin <- lm(rcc~hg+hc, data=mydata) # calcolo della regressione 

#

coefficients(reglin) # intercetta e coefficienti angolari 

confint(reglin, level=0.95) # intervalli di confidenza al 95% 

#  

La funzione lm(...) viene impiegata per calcolare la regressione lineare multipla specificando con l'espressione rcc~hg+hc la variabile dipendente (rcc), le due variabili indipendenti (hg e hc) e l'equazione 

rcc = a + b₁·hg + b₂·hc

quindi possiamo ricavare con la funzione coefficients() i valori dell'intercetta (a), del coefficiente angolare di hg (b₁) e del coefficiente angolare di hc (b₂):

(Intercept)          hg          hc 

-0.24269563  0.03283747  0.10403395 

e riportare l'equazione della regressione lineare multipla come

rcc = -0.24269563 + 0.03283747·hg + 0.10403395·hc

Successivamente con la funzione confint() sono riportati gli intervalli di confidenza al 95% dell'intercetta e dei due coefficienti angolari:

                  2.5 %     97.5 %

(Intercept) -0.53163317 0.04624192

hg          -0.02459817 0.09027311

hc           0.08267073 0.12539718

Da notare che gli intervalli di confidenza forniscono un test di significatività per intercetta e coefficienti angolari:

→ dato un intervallo di confidenza al 95% che va da -0.53163317 a 0.04624192 (che quindi include il valore 0) l'intercetta non è significativamente diversa da 0, cosa molto logica in quanto al diminuire fino ad azzerarsi dei valori di emoglobina ed ematòcrito deve necessariamente diminuire e azzerarsi il numero degli eritrociti;

→ con un intervallo di confidenza al 95% che va da -0.02459817 a 0.09027311 (che quindi include il valore 0) il coefficiente angolare della concentrazione di emoglobina hg non è significativamente diversa da 0, cosa meno logica in quanto non solo dal punto di vista biologico deve esistere una proporzionalità tra numero degli eritrociti e concentrazione dell'emoglobina, ma anche nell'analisi esplorativa preliminare dei dati riportata con lo script precedente risulta evidente la proporzionalità tra numero degli eritrociti e concentrazione dell'emoglobina che è una dato di fatto dal punto di vista biologico;

→ con un intervallo di confidenza al 95% che va da 0.08267073 a 0.12539718 (che quindi non include il valore 0) il coefficiente angolare del valore ematòcrito hc risulta significativamente diverso da 0.

Come detto all'inizio, oltre ad effettuare una analisi esplorativa preliminare dei dati è opportuno sottoporre i risultati ad una valutazione critica a posteriori, che effettuiamo ora sotto forma di una diagnostica numerica della regressione lineare multipla, copiate il terzo script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - diagnostica numerica della regressione

#

library(gvlma) # carica il pacchetto per la diagnostica numerica 

library(car) # carica il pacchetto per la grafica 

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

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

reglin <- lm(rcc~hg+hc, data=mydata) # calcolo della retta di regressione lineare multipla

#

summary(gvlma(reglin)) #  test globale per l'assunto di linearità 

outlierTest(reglin) # valore p di Bonferroni per la presenza di dati aberranti (outliers) 

#

Ecco la prima parte del riepilogo dei risultati fornito dalla funzione summary():

Call:

lm(formula = rcc ~ hg + hc, data = mydata)

Residuals:

     Min       1Q   Median       3Q      Max 

-0.55791 -0.11284 -0.00699  0.09950  0.53595 

Coefficients:

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

(Intercept) -0.24270    0.14652  -1.656   0.0992 .  

hg           0.03284    0.02913   1.127   0.2609    

hc           0.10403    0.01083   9.603   <2e-16 ***

---

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

Residual standard error: 0.1744 on 199 degrees of freedom

Multiple R-squared:  0.8565,    Adjusted R-squared:  0.855 

F-statistic: 593.8 on 2 and 199 DF,  p-value: < 2.2e-16

Il valore di probabilità Pr(>|t|) ottenuto con il test t di Student conferma (ovviamente) i risultati della significatività ricavati dagli intervalli di confidenza riportati sopra, per i quali valgono le considerazioni già fatte. E se la significatività statistica non è coerente con una realtà biologica assodata, dovrà ovviamente essere quest'ultima a prevalere: per cui assumiamo che l'equazione calcolata sia, fino a prova contraria, una descrizione sufficientemente adeguata della relazione di funzione che intercorre tra le tre variabili esaminate.

Che questo sia vero ce lo conferma la funzione gvlma(), acronimo di global validation of linear models assumptions, che fornisce questi risultati:

→ i dati sono distribuiti in modo lineare (Global Stat), assunto accettabile;

→ i dati non presentano asimmetria (Skewness), assunto accettabile;

→ i dati non presentano curtosi (Skewness), assunto accettabile;

→ nell'ambito dei valori osservati la varianza è costante (Heteroscedasticity), assunto accettabile.

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        6.5801 0.15981    Assumptions acceptable.

Skewness           0.5149 0.47304    Assumptions acceptable.

Kurtosis           0.7559 0.38461    Assumptions acceptable.

Link Function      3.9844 0.04592 Assumptions NOT satisfied!

Heteroscedasticity 1.3249 0.24972    Assumptions acceptable.

L'assenza nei dati di asimmetria e di curtosi, e la loro omoschedasticità, sono condizioni sufficienti per garantire la loro adeguata descrizione mediante una regressione lineare. Il primo test (Global Stat) di valutazione globale degli assunti di linearità, che sostanzialmente comprende tutti e tre i test precedenti in forma compatta, li conferma. In questo contesto il test rimanente e l'unico discordante (Link Function) può essere ignorato.

Con la funzione outlierTest() viene  effettuato il test di Bonferroni, che non evidenzia dati anomali (outliers) ma segnala il dato 78 come quello che si discosta maggiormente dagli altri:

No Studentized residuals with Bonferroni p < 0.05

Largest |rstudent|:

   rstudent unadjusted p-value Bonferroni p

78 -3.28855          0.0011917      0.24073

In tema di valutazione critica a posteriori, aggiungiamo ora alla diagnostica numerica anche una diagnostica grafica della regressione lineare multipla, per effettuarla copiate questo quarto script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - diagnostica grafica della regressione

#

library(car) # carica il pacchetto per la grafica

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

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

reglin <- lm(rcc~hg+hc, data=mydata) # calcolo della retta di regressione lineare multipla

#

windows() # apre e inizializza una nuova finestra grafica

qqPlot(reglin, main="Grafico dei quantili per i residui", xlab="t-quantili", ylab="Residui studentizzati") # grafico globale per linearità e dati aberranti

#

windows() # apre e inizializza una nuova finestra grafica

ceresPlots(reglin, ask=FALSE, main="Grafico di Ceres", ylab="Residui di Ceres (rcc)") # grafico di Ceres per la valutazione separata della linearità di hg e hc vs. rcc

#

windows() # apre e inizializza una nuova finestra grafica

plot(mydata$rcc, reglin$fitted.values, xlim=c(4,7), ylim=c(4,7)) # grafico valori originari vs valori stimati

abline(0, 1, lty=2) # linea di indentità y = x

#

La funzione qqPlot() riporta una distribuzione dei residui lineare, evidenziando nel contempo i dati 78 e 161 come quelli che dei quali sarebbe utile controllare la validità, ovvero che sono situati in intervalli di valori per i quali potrebbe essere opportuno acquisire più dati. 

La funzione ceresPlot() genera due grafici di Ceres separati per emoglobina ed ematòcrito, e consente di  visualizzare una relazione lineare tra ematòcrito ed eritrociti migliore di quella che intercorre tra emoglobina ed eritrociti, confermando in un certo senso i valori di significatività riportati sopra.

La diagnostica grafica della regressione lineare multipla si conclude con questo grafico minimalista, ma che entra nell'essenza dei risultati.

Il grafico, generato con la funzione plot(), riporta in ascisse i valori originari della concentrazione degli eritrociti (mydata$rcc) e in ordinate i corrispondenti valori reglin$fitted.values ricalcolati da hg e hc impiegando l'equazione della regressione lineare multipla ove

reglin$fitted.values = -0.24269563 + 0.03283747·hg + 0.10403395·hc

La linea di identità y = x che viene riportata evidenzia, senza la necessità di ulteriori statistiche, che l'informazione acquisita misurando la concentrazione dell'emoglobina e l'ematòcrito può essere compressa in un unico valore, che fornisce una stima mediamente adeguata della concentrazione degli eritrociti.

Nota: con "mediamente adeguata" si intende che, data la misura diretta della concentrazione degli eritrociti riportata in ascisse, la sua misura indiretta - ottenuta misurando emoglobina ed ematòcrito ed applicando la regressione lineare multipla - comporta una perdita di informazione [aumento dell'incertezza della misura] che si riflette nella "dispersione media" dei valori sull'asse delle ordinate. Se non ci fosse perdita di informazione i punti sarebbero perfettamente allineati sulla retta y = x.

Ora, visto che con tre variabili la cosa è ancora fattibile, non resta che aggiungere un grafico tridimensionale (3D).

Copiate quest'ultimo script, incollatelo nella Console di R e premete ↵ Invio:

# GRAFICO 3D ANIMATO con il pacchetto rgl

#  

library(rgl) # carica il pacchetto per la grafica 3D

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

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

windows() # apre e inizializza una nuova finestra grafica

#  

# grafico 3D animato

#  

plot3d(mydata$rcc, y=mydata$hg, z=mydata$hc, type="p", col="red", size=3) # realizza il grafico

play3d(spin3d(axis=c(0,0,1), rpm=10), duration=10) # anima il grafico

#  

Dopo avere caricato il pacchetto necessario con il comando library(rgl) il grafico viene realizzato con la funzione plot3d().

Quindi con la funzione play3d() viene realizzata una animazione della durata di 10 secondi (duration=10) con i parametri specificati dalla funzione spin3d() e cioè l'asse sul quale fare ruotare l'animazione, nel nostro caso  l'asse delle z (0,0,1), e il numero di rotazioni al minuto (rpm=10).

Quando l'animazione che compare nella finestra "RGL device" aperta nella Console di R dallo script si ferma, se “afferrate” il grafico 3D facendo click con il tasto sinistro del mouse e tenendolo premuto senza rilasciarlo potete ruotarlo a vostro piacimento muovendo il mouse. 

Anche se non aggiunge nulla di nuovo, il grafico [1] ci conforta documentando la linearità della relazione tra le tre variabili.

Conclusione: mettere in relazione tra di loro troppe variabili contemporaneamente può portare a risultati di difficile comprensione. Tuttavia una analisi esplorativa preliminare dei dati e una attenta valutazione critica a posteriori dei risultati possono aiutare a identificare i casi nei quali l'impiego della regressione lineare multipla può essere appropriato.

----------

[1] Trovate altre modalità di rappresentazione sotto forma di grafici tridimensionali nel post Grafici 3D.

Regressione lineare semplice non parametrica

La caratteristica saliente dei metodi di regressione lineare non parametrica è costituita dalla loro robustezza. Questo significa che si tratta di metodi che risentono poco della eventuale presenza di dati apparentemente anomali, ma che non possono essere esclusi a priori dal calcolo della regressione.

Essendo n il numero delle coppie di dati/punti esistono n · (n - 1) / 2 modi di connettere due punti qualsiasi con una retta, cioè esistono n · (n - 1) / 2 possibili coefficienti angolari [1]. Nel caso più semplice il coefficiente angolare b della retta di regressione è la mediana di questi valori e l'intercetta a è la mediana degli n valori possibili calcolati mediante il coefficiente angolare trovato. Questi sono concettualmente i principi alla base della regressione lineare non parametrica, con possibili varianti da metodo a metodo.

Qui vediamo, confrontati con la regressione lineare semplice parametrica (regressione lineare ordinaria), tre modelli di regressione lineare non parametrica:

→ la regressione lineare di Thiel-Sen;

→ la regressione lineare di Siegel;

→ la regressione quantilica.

Da notare che questi tre modelli di regressione lineare non parametrica assumono la x come variabile indipendente analogamente alla regressione lineare ordinaria.

Le funzioni necessarie per il calcolo delle tre regressioni lineari non parametriche sono incluse in due pacchetti aggiuntivi, mblm e quantreg, che devono essere scaricati dal CRAN. Trovate come al solito la documentazione dei pacchetti e in particolare il loro manuale di riferimento sul repository della documentazione [2].

Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE ORDINARIA E NON PARAMETRICA DI THEIL-SEN, DI SIEGEL, QUANTILICA

#

library(mblm) # carica il pacchetto per le regressioni di Theil-Sen e Siegel

library(quantreg) # carica il pacchetto per la regressione quantilica

#

var_x <- c(18, 46, 104, 72, 24, 63, 23, 58, 36, 82, 94) # variabile in ascisse

var_y <- c(21, 42, 109, 64, 102, 78, 19, 66, 42, 89, 91) # variabile in ordinate

#

regpar_yx <- lm(var_y ~ var_x) # regressione lineare ordinaria x variabile indipendente

a_yx <- regpar_yx$coefficients[1] # intercetta a

b_yx <- regpar_yx$coefficients[2] # coefficiente angolare b

#

reg_TS <- mblm(var_y ~ var_x, repeated=FALSE) # regressione lineare di Theil-Sen

a_TS <- reg_TS$coefficients[1] # intercetta a

b_TS <- reg_TS$coefficients[2] # coefficiente angolare b

#

reg_S = mblm(var_y ~ var_x, repeated=TRUE) # regressione lineare di Siegel

a_S <- reg_S$coefficients[1] # intercetta a

b_S <- reg_S$coefficients[2] # coefficiente angolare b

#

reg_q <- rq(var_y ~ var_x, tau=0.5) # regressione lineare quantilica

a_q <- reg_q$coefficients[1] # intercetta a

b_q <- reg_q$coefficients[2] # coefficiente angolare b

#

# traccia il grafico dei dati con le rette di regressione e aggiunge una legenda

#

windows() # apre una nuova finestra

plot(var_x, var_y, xlim = c(0,120), ylim = c(0,120), pch=1, xlab="Variabile x", ylab="Variabile y", main="Regressioni ordinaria, Theil-Sen, Siegel, quantilica", cex.main = 0.9) # grafico dei dati

abline(a_yx, b_yx, col="blue", lty=4) # retta regressione ordinaria

abline(a_TS, b_TS, col="red", lty=4) # retta regressione di Theil-Sen

abline(a_S, b_S, col="green", lty=1) # retta regressione di Siegel

abline(a_q, b_q, col="orange", lty=1) # retta regressione quantilica

legend(60, 30, legend=c("Regressione ordinaria", "Regressione di Theil-Sen", "Regressione di Siegel", "Regressione quantilica"), col=c("blue", "red", "green", "orange"), lty=c(4,4,1,1), cex=0.8) # aggiunge al grafico la legenda

#

# crea una tabella con intercetta e coefficiente angolare delle quattro rette di regressione

#

cells <- c(a_yx,b_yx,a_TS,b_TS,a_S,b_S,a_q,b_q) # genera l'array cells con i valori numerici di a e di b

rnames <- c("Regressione ordinaria", "Regressione di Theil-Sen", "Regressione di Siegel", "Regressione quantilica") # genera l'array rnames con i nomi delle righe

cnames <- c("Intercetta a", "Coefficiente angolare b") # genera l'array cnames con i nomi delle colonne

tabris <- matrix(cells, nrow=4, ncol=2, byrow=TRUE, dimnames=list(rnames, cnames)) # costruisce la tabella con i risultati a partire dagli array cells, rnames, cnames

tabris # mostra la tabella con i risultati

#

Questo è il grafico di dispersione (xy) dei dati, con le rette calcolate:

I risultati delle quattro equazioni

y = a + b · x

sono così riportati nella Console di R:

> tabris # mostra la tabella con i risultati

                         Intercetta a Coefficiente angolare b

Regressione ordinaria       24.078677               0.7389267

Regressione di Theil-Sen     6.529412               0.9852941

Regressione di Siegel        5.873402               1.0042750

Regressione quantilica       2.581395               1.0232558

Le cose da notare sono:

→ è evidente la presenza di un dato che si discosta in modo importante da tutti gli altri;

→ i risultati dei tre metodi non parametrici sono simili tra loro;

→ i risultati delle regressione ordinaria (parametrica) si discostano da quelli delle regressioni non parametriche.

Per verificare quanto il dato, chiaramente anomalo, influenzi le regressioni, eseguiamo nuovamente lo script, sostituendo la terza e la quarta riga con queste due, nelle quali è stato eliminato il dato, il quinto, quello con le coordinate (x,y) uguali a (24, 102):

var_x <- c(18, 46, 104, 72, 63, 23, 58, 36, 82, 94) # variabile in ascisse

var_y <- c(21, 42, 109, 64, 78, 19, 66, 42, 89, 91) # variabile in ordinate

Ora i grafici delle rette sono praticamente sovrapposti (e difficilmente distinguibili tra loro)

mentre i risultati sono diventati i seguenti:

> tabris # mostra la tabella con i risultati

                         Intercetta a Coefficiente angolare b

Regressione ordinaria        1.337096                1.019512

Regressione di Theil-Sen     2.729167                1.020833

Regressione di Siegel        2.654334                1.022497

Regressione quantilica       2.581395                1.023256

Come si vede dal confronto con i risultati precedenti, con l'eliminazione del dato anomalo i risultati della regressione ordinaria sono cambiati: e questa fornisce ora risultati sovrapponibili a quelli dei tre metodi non parametrici. Poco sono cambiati, rispetti ai precedenti, i risultati della regressione di Thiel-Sen e di Siegel, mentre la regressione quantilica ha addirittura fornito identici valori nei due casi, un fatto molto interessante. I metodi non parametrici confermano la loro robustezza, mentre viene documentato quanto un singolo dato anomalo possa influire sulla regressione lineare ordinaria.

----------

[1] Per 100 coppie di dati il numero dei coefficienti angolari (sui quali calcolare la mediana) è uguale quindi a 4 950. Con 500 dati passiamo a 124 740. Con 1 000 dati passiamo a 499 500. 

[2] Available CRAN Packages By Name.

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