Regressione lineare media geometrica

Per illustrare il tema, è opportuno ricordare gli assunti alla base della regressione lineare ordinaria cioè della regressione lineare x variabile indipendente con il metodo dei minimi quadrati [1] [2] che esprime la y in funzione della x: 

→ in corrispondenza di una serie di valori dati della x abbiamo altrettanti valori misurati della y;

→ in quanto misurati i valori della y sono affetti da un errore casuale che ha una distribuzione gaussiana;

→ per tutti i valori della x la distribuzione dell'errore casuale della y ha la medesima varianza;

→ la media delle distribuzioni gaussiane degli errori della y viene allineata su una retta avente equazione y = a₁ + b₁·x 

Esprimendo la x in funzione della y e avendo per una serie di valori dati della y altrettanti valori misurati della x si ha la regressione lineare y variabile indipendente, mentre l'equazione della retta assume la forma x = a₂ + b₂·y 

In condizioni usuali questo secondo modo di analizzare i dati non viene impiegato: per convenzione la variabile indipendente viene sempre riportata in ascisse e la variabile dipendente viene riportata in ordinate, quindi si procede con il calcolo della regressione lineare x variabile indipendente, o regressione lineare ordinaria, quella cui di fatto ci si riferisce quando si parla semplicemente di regressione lineare.

Tuttavia i due modi di vedere riportati sopra diventano fondamentali quando entrambe le variabili sono il risultato di un processo di misura, e quindi non esiste una variabile indipendente. In questo caso la soluzione più semplice ma anche la più robusta, che assegna lo stesso peso agli errori di misura della x e a quelli della y, è rappresentata dalla regressione lineare che impiega come coefficiente angolare b la media geometrica dei due coefficienti angolari, x variabile indipendente e y variabile indipendente [3],[4].

Nel set di dati ais contenuto nel pacchetto DAAG [5] abbiamo 202 valori di concentrazione degli eritrociti, espressa in migliaia di miliardi per litro di sangue (10¹²/L) – o milioni per microlitro di sangue (10⁶/µL), le due espressioni sono numericamente equivalenti – e i corrispondenti valori di emoglobina, la cui concentrazione è espressa in grammi per decilitro di sangue (g/dL). Entrambe le grandezze sono misurate quindi nessuna delle due può essere assunta come variabile indipendente, e procediamo con il calcolo della regressione lineare media geometrica.

Innanzitutto calcoliamo la regressione lineare x variabile indipendente y = a₁ + b₁·x 

 

# REGRESSIONE LINEARE media geometrica 

#

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

x <- ais$rcc # eritrociti in ascisse

y <- ais$hg # emoglobina in ordinate

#

# regressione y = a + b·x (rcc ascisse, hg in ordinate)

regpar <- lm(y ~ x) # regressione lineare ordinaria 

#

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

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

#

Dopo avere caricato il pacchetto che include il set di dati ais per chiarezza semantica gli eritrociti sono ridenominati x e l'emoglobina è ridenominata y. 

Con la funzione lm() viene effettuato il calcolo della regressione lineare ordinaria con la y espressa in funzione della x (y ~ x) e sono memorizzati i valori così calcolati del coefficiente angolare (a_yx) e dell'intercetta (b_yx) della retta di regressione y = a₁ + b₁·x

Calcoliamo ora la regressione lineare y variabile indipendente x = a₂ + b₂·y in due passaggi: 

→ prima calcoliamo la regressione lineare ordinaria invertendo tra loro le due variabili; 

→ poi ripristiniamo la x in ascisse e la y in ordinate in modo che il grafico possa essere sovrapposto a quello della regressione lineare x variabile indipendente.

  

#

# regressione x = a + b·y (hg in ascisse, rcc in ordinate)

regpar_1 <- lm(x ~ y) # regressione lineare ordinaria (con variabili invertite)

#

a_xy <- - regpar_1$coefficients[1] / regpar_1$coefficients[2] # intercetta a ripristinando rcc in ascisse e hg in ordinate

b_xy <- 1 / regpar_1$coefficients[2] # coefficiente angolare b ripristinando rcc in ascisse e hg in ordinate)

#

Con la funzione lm() viene effettuato il calcolo della regressione lineare ordinaria con la x espressa in funzione della y (x ~ y), quindi sono memorizzati i valori del coefficiente angolare a_xy = - a₂/b₂ e dell'intercetta b_xy = 1/b₂ ricavati dalla retta di regressione x = a₂ + b₂·y esplicitando quest'ultima rispetto alla y essendo quindi

x = a₂ + b₂·y ovvero b₂·y = - a₂ + x ovvero  y = - a₂/b₂ + (1/b₂)·x 

A questo punto per la regressione lineare media geometrica con equazione y = a₃ + b₃·x si calcola il coefficiente angolare b₃ come media geometrica dei due coefficienti angolari b_yx (b₁) e b_xy (1/b₂) ovvero come radice quadrata del loro prodotto (prima riga) 

 

#

# regressione media geometrica 

b_mg <- sqrt(b_yx*b_xy) # coefficiente angolare b media geometrica

a_mg <- mean(y) - (b_mg * (mean(x))) # intercetta a = ȳ - b·x̄

#

quindi dal coefficiente angolare b₃ così calcolato si ricava (seconda riga) l'intercetta a₃ come ȳ - b₃·x̄ (media delle y - b₃ volte la media delle x).

Per rappresentare il tutto riportiamo infine un grafico con le tre rette di regressione e aggiungiamo una tabella con i tre coefficienti angolari e le tre intercette calcolate. 

# grafico

plot(x, y, xlim = c(0,8), ylim = c(0,20), xaxs="i", yaxs="i", xlab="rcc", ylab="hg", cex.main = 0.9) # grafico dei dati

abline(a_yx, b_yx, col="blue", lty=4) # retta di regressione x variabile indipendente

abline(a_xy, b_xy, col="red", lty=4) # retta di regressione y variabile indipendente

abline(a_mg, b_mg, col="green", lty=1) # retta di regressione media geometrica

# legenda

legend(3.5, 7.5, legend=c("Regressione x variabile indipendente", "Regressione y variabile indipendente", "Regressione media geometrica"), col=c("blue", "red", "green"), lty=c(4,4,1), cex=0.8) # aggiunge al grafico la legenda

#

# genera la tabella con i risultati delle regressioni

interc_a <- c(a_yx, a_xy, a_mg) # array con le intercette

cofang_b <- c(b_yx, b_xy, b_mg) # array con i coefficienti angolari

tabris <- data.frame(interc_a, cofang_b) # combina gli array in una tabella

rownames(tabris) <- c("Regressione x variabile indipendente", "Regressione y variabile indipendente", "Regressione media geometrica") # aggiunge i nomi delle righe

colnames(tabris) <- c("Intercetta a", "Coefficiente angolare b") # aggiunge i nomi delle colonne

tabris # mostra la tabella con i risultati

#

Con la funzione plot() viene realizzato il grafico impiegando, accanto agli usuali argomenti, gli argomenti xaxs="i" e yaxs="i" che allineano inizio e fine delle scale al riquadro grafico (gli argomenti di default sono xaxs="r" e yaxs="r",  e lasciano il 6% di spazio aggiuntivo tra inizio e fine delle scale e riquadro grafico, si può fare una prova con questi valori per visualizzare la differenza).

Con la funzione abline() sono poi aggiunte al grafico le rette di regressione. e con la funzione legend() una legenda per identificarle.

Dai valori delle intercette a_yx, a_xy e a_mg e dai valori dei coefficienti angolari b_yx, b_xy e b_mg viene infine costruita la tabella con la sintesi dei risultati.

> tabris # mostra la tabella con i risultati

                                     Intercetta a Coefficiente angolare b

Regressione x variabile indipendente    2.0897326                2.644125

Regressione y variabile indipendente   -1.2275298                3.347141

Regressione media geometrica            0.5287507                2.974938

Si può fare una verifica dei risultati impiegando il pacchetto lmodel2 e, dopo averlo installato, eseguendo 

 

#

# regressione lineare con più modelli

library(lmodel2)

lmodel2(y ~ x, nperm=1000)

#

che dopo alcune righe di commenti ci fornisce

Regression results

  Method  Intercept    Slope Angle (degrees) P-perm (1-tailed)

1    OLS  2.0897326 2.644125        69.28353       0.000999001

2     MA -0.8840492 3.274348        73.01710       0.000999001

3    SMA  0.5287507 2.974938        71.42037                NA

Confidence intervals

  Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope

1    OLS      1.1885079       2.9909574   2.454020    2.834229

2     MA     -2.0733368       0.1584401   3.053417    3.526390

3    SMA     -0.3969115       1.3971493   2.790902    3.171111

ovvero intercetta e coefficiente angolare della regressione lineare ordinaria (OLS) e della regressione lineare media geometrica con le variabili standardizzate (SMA) – identici a quelli ottenuti con il nostro codice – e l'intervallo di confidenza al 95% dell'intercetta della regressione lineare media geometrica -0.3969115  1.3971493 che ci conferma che l'intercetta 0.5287507 non è significativamente diversa da 0 (zero) in quanto entrambi i valori sono inclusi nell'intervallo.

Domanda: ma perché non impiegare addirittura subito il pacchetto lmodel2? Risposta: perché si sono voluti illustrare i concetti e i passaggi alla base della regressione lineare media geometrica. Una volta confermato che il pacchetto ci fornisce i risultati attesi, nulla osta che si proceda per vie brevi e si impieghino le sue funzioni, ma a questo punto con una maggior consapevolezza del tema.

Conclusione:

→ l'emoglobina, la proteina che trasporta l'ossigeno dai polmoni ai tessuti, è contenuta esclusivamente nei globuli rossi (o eritrociti) del sangue e pertanto ci si attende una relazione di proporzionalità lineare tra la concentrazione di questi e la concentrazione dell'emoglobina;

→ dati empirici inoppugnabili ci dicono che in assenza di eritrociti deve necessariamente essere assente l'emoglobina, quindi la retta di regressione dovrebbe passare per l'origine degli assi (x,y = 0,0).

→ eritrociti ed emoglobina sono entrambi il risultato di misure quindi nessuna delle due grandezze può essere assunta come variabile indipendente;

→ la regressione media geometrica, che media tra regressione x variabile indipendente e regressione y variabile indipendente assegnando uguale peso all'errore di misura loro proprio, fornisce un valore dell'intercetta che non è significativamente diverso da 0; il fatto che il valore non sia esattamente uguale a 0 non dipende dal modello di regressione, bensì risiede nell'informazione contenuta nei dati, che per il fatto di essere molto concentrati in un ambito di valori ristretto e molto lontani dallo 0, forniscono una informazione insufficiente per rendere accurata una estrapolazione a così elevata distanza dai dati sperimentali;

→ se, come bisognerebbe fare, si evita l'estrapolazione e si traggono le conclusioni all'interno dei valori effettivamente osservati, in assenza di una variabile indipendente la regressione media geometrica è quella che fornisce i risultati concettualmente più validi;

→ a fronte della granitica certezza che due variabili debbano raggiungere congiuntamente il valore 0 è possibile impiegare un modello di regressione alternativo avente equazione della retta y = b·x che non prevede una intercetta in quanto impone il passaggio della retta per l'origine degli assi [6].

----------

[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] Ling Leng et al. Ordinary least square regression, orthogonal regression, geometric mean regression and their applications in aerosol science. J. Phys. 2007;Conf. Ser. 78 012084. DOI: 10.1088/1742-6596/78/1/012084

https://iopscience.iop.org/article/10.1088/1742-6596/78/1/012084

[4] Shaoji Xu. A Property of Geometric Mean Regression. The American Statistician 2014;68(4);277-281. DOI: 10.1080/00031305.2014.962763

https://www.tandfonline.com/doi/abs/10.1080/00031305.2014.962763 

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

[6] Vedere il post Regressione lineare passante per l'origine.