Regressione lineare bayesiana

Il teorema di Bayes è in sostanza un motore inferenziale in grado di aggiornare la probabilità a priori di un evento con le evidenze empiriche per fornire la probabilità a posteriori dell'evento stesso (date le evidenze acquisite).

In questo senso l'approccio bayesiano può essere impiegato anche per calcolare una regressione lineare. Ed è possibile farlo a partire da un'ipotesi minima, assumendo ad esempio che per la retta di equazione y = a + b·x sia a=0 e sia b=0 e quindi aggiornando questa "credenza a priori" sulla base dei dati (le evidenze) che sono stati raccolti. 

Vediamo a titolo di esempio come calcolare la regressione lineare bayesiana impiegando le funzioni predisposte in due pacchetti aggiuntivi.

Regressione lineare bayesiana con il pacchetto brms

Come ricordato anche nello script, oltre a installare il pacchetto brms in questo caso è necessario installare dal sito di R il software Rtools che contiene il compilatore C++ che genera il codice eseguibile impiegato nei calcoli.

 

# REGRESSIONE LINEARE BAYESIANA y = a + b · x (con il pacchetto brms)

#

# Richiede download e installazione della versione di Rtools appropriata da

# https://cran.r-project.org/bin/windows/Rtools/

#

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

library(brms) # carica il pacchetto per la regressione lineare bayesiana

#

apriori <- c(prior(normal(0, 10), class="Intercept"), prior(normal(0, 10), class="b"), prior(student_t(3, 0, 10), class="sigma")) # definisce gli 'a priori'

#

bayes_reg <- brm(formula=hg~rcc, data=ais, prior=apriori, family=gaussian(), chains=4, iter=4000, warmup=2000, seed=123) # calcola la regressione lineare bayesiana

#

summary(bayes_reg) # sommario dei risultati

#

Come dati di esempio impieghiamo i soliti 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), valori che sono contenuti nel set di dati ais del pacchetto DAAG [1] che al bisogno deve essere scaricato e installato insieme al pacchetto brms.

Per il calcolo della regressione bayesiana è necessario innanzitutto impiegare la funzione priors() per fissare la probabilità a priori delle grandezze da calcolare, in questo modo:

→ per l'intercetta a (class="Intercept") della retta  y = a + b · x assumiamo a priori una distribuzione gaussiana (normal) con media 0 e deviazione standard 10;

→ lo stesso facciamo per il coefficiente angolare b (class="b");

→ per la deviazione standard residua (sigma) assumiamo a priori una distribuzione t di Student (student_t) con 3 gradi di libertà, con una media di 0 e una deviazione standard di 10.

Una media di 0 e una incertezza così ampia lasciano, nel calcolo dei risultati, completamente spazio all'informazione contenuta nei dati forniti, ovvero, se si preferisce, alle evidenze fornite dai dati.

Quindi con la funzione brm() viene calcolata la regressione bayesiana sui valori di emoglobina in funzione degli eritrociti (hg~rcc) contenuti nel set di dati ais impiegando i valori a priori (apriori) definiti alla riga precedente, mentre per gli altri parametri impiegati e per il metodo di calcolo (metodo Monte Carlo basato su catena di Markov) si rimanda al manuale del pacchetto [2].

Attenzione: il metodo impiegato richiede tempo, orientativamente con una macchina con processore Intel® Core™ i7 quasi due minuti, che potrebbero aumentare nel caso di processori meno potenti e RAM limitate. 

Con la funzione summary() viene riportata una sintesi dei risultati.

> summary(bayes_reg) # sommario dei risultati

 Family: gaussian 

  Links: mu = identity 

Formula: hg ~ rcc 

   Data: ais (Number of observations: 202) 

  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;

         total post-warmup draws = 8000

Regression Coefficients:

          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS

Intercept     2.10      0.46     1.20     2.99 1.00     8802     5734

rcc           2.64      0.10     2.45     2.83 1.00     8804     5739

Further Distributional Parameters:

      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS

sigma     0.63      0.03     0.57     0.70 1.00     8965     6033

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS

and Tail_ESS are effective sample size measures, and Rhat is the potential

scale reduction factor on split chains (at convergence, Rhat = 1).

Per giudicare il risultato ottenuto che è quindi y = 2.10 + 2.64 · x il metodo più semplice e immediato è confrontarlo con la regressione lineare ordinaria x variabile indipendente 

 

#

lm(ais$hg ~ ais$rcc) # regressione lineare ordinaria

#

che ci riporta il risultato praticamente identico y = 2.09 + 2.64 · x 

> lm(ais$hg ~ ais$rcc)

Call:

lm(formula = ais$hg ~ ais$rcc)

Coefficients:

(Intercept)      ais$rcc  

      2.090        2.644 

Il pacchetto consente anche di riportare un grafico con i dati (points=TRUE), con la retta di regressione calcolata e con i suoi limiti di confidenza

 

#

plot(conditional_effects(bayes_reg), points=TRUE) # grafico della regressione con i limiti di confidenza

#

Nota: il messaggio che compare dopo la generazione del grafico è determinato dalla particolare gestione dell'estetica da parte delle funzioni grafiche impiegate, e non ha rilevanza pratica.

Da ultimo con il codice che segue è possibile generare anche un grafico (qui non riportato) che illustra le distribuzioni calcolate a posteriori di a, b e sigma e l'andamento dei calcoli

 

#

windows()

plot(bayes_reg) # grafico delle distribuzioni 'a posteriori'

# 

per i quali si rimanda nuovamente a quanto illustrato nella documentazione del pacchetto [2]. 

Regressione lineare bayesiana con il pacchetto BAS

Con il pacchetto BAS il codice che è necessario impiegare è ancora più semplice.

 

# REGRESSIONE LINEARE BAYESIANA y = a + b · x (con il pacchetto BAS)

#

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

library(BAS) # carica il pacchetto per la regressione lineare bayesiana

#

bayes_model <- bas.lm(hg~rcc, data=ais) # calcola la regressione lineare bayesiana

#

a <- round(predict(bayes_model, newdata=data.frame(rcc=0))$fit[1], digits=2) # intercetta a

b <- round(coef(bayes_model)$postmean[2], digits=2) # coefficiente angolare b

paste("y =", a, "+", b, "· x") # equazione della retta

#

plot(ais$rcc, ais$hg, pch = 19, col = "gray") # grafico dei dati

abline(a, b, col="blue", lwd=2) # retta di regressione

#

Dopo avere caricato i due pacchetti necessari il calcolo viene effettuato con la funzione bas.lm() che richiede solamente di specificare le variabili (hg~rcc) e il set di dati (ais) che le contiene.

Quindi:

→ l'intercetta a viene calcolata con la funzione predict() come valore della y corrispondente al valore 0 della x (data.frame(rcc=0)) e questo artifizio si rende necessario in quanto la funzione coef(bayes_model) fornisce come intercetta la media delle y – per chiarimenti vedere la documentazione del pacchetto [3];

→ come valore del coefficiente angolare b viene riportato quello fornito dalla funzione coef(bayes_model);

→ a e b sono arrotondati a due cifre decimali con la funzione round();

→ con la funzione paste() intercetta e coefficiente angolare sono combinati nell'equazione della retta di regressione 

> paste("y =", a, "+", b, "· x") # equazione della retta

[1] "y = 2.12 + 2.64 · x"

che nella sostanza non differisce da quelle riportate sopra per la regressione lineare bayesiana e per la regressione lineare ordinaria (per capirci e senza la necessità di impiegare test di significatività, tutti e tre i coefficienti angolari sono uguali a 2.64 e se arrotondiamo i risultati a una cifra decimale abbiamo tutte e tre le intercette uguali a 2.1).

A questo punto non resta che impiegare la funzione plot() per rappresentare il gafico e tracciare la retta di regressione con la funzione abline().

 

Conclusione: ci si potrebbe stupire del fatto che la regressione lineare ordinaria, che impiega un approccio analitico – il metodo dei minimi quadrati, così denominato da Legendre che lo descrisse formalmente per la prima volta nel 1805 [4], mentre Gauss, che polemizzò con Legendre dicendo di essere stato lui a idearlo qualche anno prima, lo pubblicò solamente nel 1809 [5] – fornisca risultati identici ai due metodi bayesiani, che impiegano un approccio probabilistico. 

Ma c'è una ragione sottile e insieme profonda:

→ quando l'informazione in ingresso è scarsa – perché i dati sono pochi e sono affetti da un consistente errore di misura – i risultati sono influenzati dagli assunti alla base del modello matematico-statistico impiegato; 

→ tanto maggiore è l'informazione contenuta nei dati – cioè quando i dati sono molti e sono affetti da incertezza di misura ridotta (come nel nostro caso) – tanto meno gli assunti influiscono sui risultati, che anche se ottenuti con metodi differenti finiscono per convergere (stiamo ovviamente parlando di metodi rigorosi, scientificamente fondati e matematicamente dimostrati). 

Il che ci ricorda ancora una volta, alla luce del noto adagio "garbage in, garbage out", che la numerosià campionaria, così come l'errore di misura, sono componenti fondamentali dell'informazione fornita dai dati (le evidenze) e che la qualità dei risultati in uscita non può essere migliore della qualità dei dati in ingresso. 

----------

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

[2] Vedere la documentazione del pacchetto brms:

https://cran.r-project.org/web/packages/brms/brms.pdf

https://cran.r-project.org/web/packages/brms/vignettes/brms_overview.pdf

[3] BAS: Bayesian Variable Selection and Model Averaging using Bayesian Adaptive Sampling.

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

[4] Nouvelles méthodes pour la détermination des orbites des comètes, par A.-M.

Legendre, 1805 ...

"Nuovi metodi per la determinazione delle orbite delle comete"

https://gallica.bnf.fr/ark:/12148/bpt6k15209372/

"... la méthode qui me paroit la plus simple et la plus générale, consiste à rendre minimum la somme des quarrés des erreurs. On obtient ainsi autant d'équations qu'il y a de coefficiens inconnus; ce qui achève de déterminer tous les élémens de !'orbite. Comme la méthode dont je viens de parler, et que j'appelle méthode des moindres quarrés, peut ètre d'une grande utilité dans toutes !es questions de physique et d'astronomie où il s'agil de tirer de l'observation les résultats les plus exacts qu'elle peut offrir; j'ai ajouté, dans une appendice, des détails particuliers sur cette méthode ..."

"... il metodo che mi sembra più semplice e generale consiste nel minimizzare la somma dei quadrati degli errori. Questo fornisce tante equazioni quanti sono i coefficienti incogniti; ccosa che consente la determinazione di tutti gli elementi dell'orbita. Poiché il metodo che ho appena descritto, che chiamo metodo dei minimi quadrati, può essere di grande utilità in tutte le questioni di fisica e astronomia nelle quali si tratta di ricavare dalle osservazioni i risultati più accurati che questa possa offrire, ho aggiunto, in appendice, dettagli specifici su questo metodo..."

[5] Theoria motus corporum coelestium in sectionibus conicis Solem ambientium. Auctore Carolo Friderico Gauss, 1809, pp. 220-221.

"Teoria del moto dei corpi celesti in sezioni coniche attorno al Sole" 

https://archive.org/details/bub_gb_ORUOAAAAQAAJ

Grafici di dispersione (scatterplot) per fattore con ggplot

La statistica è l'arte sottile che consente di ricavare informazione dai dati. Vediamo come ciò sia possibile impiegando anche semplici grafici di dispersione ma in questo caso separati per fattore.

Il punto di partenza è un grafico xy della distribuzione dei valori di concentrazione nel sangue dei globuli rossi (in ascisse) e dell'emoglobina (in ordinate) misurati in 202 atleti australiani di entrambi i sessi che praticano dieci differenti sport e contenuti nel set di dati ais.

 

# SCATTERPLOT semplice con ggplot

#

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

library(ggplot2) # carica il pacchetto per la grafica

#

# scatterplot semplice

ggplot(ais, aes(x=rcc, y=hg, color=sport)) + geom_point()

#

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

Per realizzare il grafico viene impiegata la funzione ggplot() che prevede come argomenti:

→ il nome della tabella che contiene i dati (ais);

→ la funzione aes() che specifica la variabile da riportare in ascisse (rcc) e la variabile da riportare in ordinate (hg), che sono rispettivamente la concentrazione degli eritrociti (o globuli rossi) nel sangue espressa in milioni per microlitro di sangue (10^6/µL) e la concentrazione di emoglobina espressa in grammi per decilitro di sangue (g/dL);

→ il colore dei punti (color) separato per ciascuno degli sport praticati e codificati nella variabile sport.

Alla funzione ggplot() viene quindi collegata con il segno più [+] la funzione geom_point() che realizza il grafico impiegando i valori previsti di default e individuando gli sport praticati con un colore.

Visto il grado di sovrapposizione dei dati, vogliamo provare a separarli generando dieci scatterplot, uno per ciascuno degli sport. Ecco lo script necessario e i grafici risultanti.

 

# SCATTERPLOT separati per fattore

#

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

library(ggplot2) # carica il pacchetto per la grafica

#

# scatterplot separati per il fattore sport

ggplot(ais, aes(x=rcc, y=hg, color=sport)) + geom_point() + facet_wrap(~sport) 

#

Quindi per scomporre i dati nei dieci sottoinsiemi rappresentativi di ciascuno sport è stato sufficiente, al codice necessario per realizzare il grafico di dispersione che riporta tutti i dati, aggiungere (+) la funzione facet_wrap() specificando che i dati devono essere suddivisi per sport (~sport).

Possiamo anche migliorare il grafico come segue.

 

# SCATTERPLOT riorganizzati per righe

#

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

library(ggplot2) # carica il pacchetto per la grafica

#

# scatterplot riorganizzati per righe

ggplot(ais, aes(x=rcc, y=hg, color=sport)) + geom_point(show.legend=FALSE) + facet_wrap(~sport, nrow=2) 

#

In questo caso:

→ visto che la rappresentazione fornita contiene già le etichette che individuano i dati, è stata eliminata la legenda (show.legend=FALSE);

→ i grafici sono stati riorganizzati su due righe con nrow=2, con questo risultato.

Dei nostri dati abbiamo quindi realizzato una grafica sintetica (primo script) e poi abbiamo visto come realizzare una grafica analitica: ma che dire se potessimo visualizzarle contemporaneamente entrambe? Ecco come è possibile farlo.

 

# SCATTERPLOT sovrapposti ai dati globali

#

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

library(ggplot2) # carica il pacchetto per la grafica

#

# scatterplot separati sovrapposti alla distribuzione globale dei dati

ggplot(ais, aes(x=rcc, y=hg, color=sport)) + geom_point(data=transform(ais, sport=NULL), colour="grey82") + geom_point(show.legend = FALSE) + facet_wrap(~sport, nrow=2, strip.position="bottom")

#

Nota: i grafici mettono in rilievo l'informazione che la distribuzione dei valori non è omogenea in tutti gli sport – ad esempio gli atleti che praticano lo sport Netball hanno concentrazioni di eritrociti ed emoglobina relativamente basse e all'estremo opposto di quelle relativamente alte che si riscontrano negli atleti che praticano pallanuoto (W_Polo).

In questo caso viene aggiunta una nuova funzione geom_point() che elimina dai dati la variabile sport (sport=NULL) mediante la funzione transform() per riportarli tutti nel grafico – in colore grigio (colour="grey82") – prima che ad essi con le successive funzioni  geom_point() e facet_wrap() siano sovrapposti separatamente i dati dei singoli sport. 

Da notare che la funzione facet_wrap() tra le altre cose:

→  consente con l'argomento strip.position di personalizzare la posizione delle etichette, con il valore predefinito "top", che può essere cambiato in "bottom" (come in questo caso), "left" o "right";

→  consente con l'argomento scales di specificare il tipo di scale impiegate, con il valore predefinito "fixed" che assume come in questo caso scale degli assi x e y uguali in tutti i grafici, e con i valori "free", "free_x" o "free_y" da impiegare qualora siano previste scale differenti degli assi.

Per completare l'analisi grafica dei dati la eseguiamo anche per il sesso degli atleti sostituendo sport con sex, il nome della variabile fattore che contiene l'indicazione del sesso dell'atleta [1].

 

# SCATTERPLOT sovrapposti ai dati globali

#

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

library(ggplot2) # carica il pacchetto per la grafica

#

# scatterplot separati sovrapposti alla distribuzione globale dei dati

ggplot(ais, aes(x=rcc, y=hg, color=sex)) + geom_point(data=transform(ais, sex=NULL), colour="grey82") + geom_point(show.legend=FALSE) + facet_wrap(~sex, nrow=2, strip.position="bottom")

#

In questo caso il grafico mette in evidenza che la distribuzione del valori non è omogenea nei due sessi, un fatto che rispecchia i differenti "valori normali" dei globuli rossi e dell'emoglobina nelle donne e negli uomini (valori che sono fisiologicamente mediamente inferiori nelle prime rispetto ai secondi).

Conclusione: opportunamente impiegati i grafici xy consentono di ricavare dai dati informazioni che vanno al di la della semplice relazione tra due variabili.

---------

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