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).
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 di una serie di dati (le evidenze) che sono stati raccolti.
Vediamo a titolo di esempio come calcolare la regressione lineare bayesiana impiegando due differenti algoritmi impiegati in due pacchetti aggiuntivi.
Regressione lineare bayesiana con il pacchetto brms
Come ricordato anche nello script, oltre a installare il pacchetto [1] in questo caso è necessario installare in via preliminare dal sito di R il software Rtools che contiene il compilatore C++ necessario per generare il codice eseguibile impiegato nel calcolo.
# 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à a sempre con una media di 0 e una deviazione standard di 10.
Una media di 0 e una incertezza così ampia lasciano, per il calcolo dei risultati, completamente spazio all'informazione contenuta nei dati forniti, ovvero, se si preferisce, alle evidenze fornite dai dati.
Ccon la funzione brms() viene calcola la regressione bayesiana sui valori di emoglobina in funzione degli eritrociti (hg~rcc) contenuti nel set di dati ais impiegando gli a priori (apriori) definiti alla riga di codice precedente. Per i parametri impiegati nel metodo Monte Carlo utilizzato si rimanda al manuale del pacchetto [1].
Attenzione: il metodo impiegato richiede tempo, orientativamente con una macchina con processore Intel® Core™ i7 un paio di minuti, che potrebbero aumentare nel caso di processori meno potenti e RAM limitate.
Con la funzione summary() viene riportata una sintesi dei risultati.
> # sommario dei risultati
> summary(bayes_model)
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.08 0.46 1.16 2.98 1.00 8484 5871
rcc 2.65 0.10 2.45 2.84 1.00 8531 5659
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 8528 6201
Per giudicare il risultato ottenuto con la regressione lineare bayesiana che è quindi y = 2.08 + 2.65 · x il metodo più semplice e immediato è calcolare 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) e con la retta di regressione calcolata, inclusi 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 dipende dalla gestione dell'estetica del pacchetto impiegato per la grafica, ma non ha alcuna rilevanza.
Da ultimo è possibile generare anche un grafico (qui non riportato) che illustra le distribuzioni calcolate a posteriori di a, b e sigma
#
windows()
plot(bayes_reg) # grafico delle distribuzioni 'a posteriori'
#
per le quali si rimanda nuovamente a quanto riportato nella documentazione del pacchetto [1].
Regressione lineare bayesiana con il pacchetto BAS
Con il pacchetto BAS il codice che è necessario impiegare è addirittura 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:
→ con la funzione predict() viene calcolata l'intercetta a come valore della y corrispondente al valore 0 della x (data.frame(rcc=0));
→ il valore del coefficiente angolare b viene riportato mediante la funzione coef();
→ 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 è sostanzialmente identica a quelle riportate sopra per la regressione lineare bayesiana e per la regressione lineare ordinaria.
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, come nel nostro caso, il metodo impiegato nella regressione lineare ordinaria, un metodo analitico – il metodo dei minimi quadrati, così denominato nel 1805 da Legendre che lo descrisse formalmente per la prima volta [2] anche se Gauss avanzò la pretesa di averlo ideato qualche anno prima (ma senza pubblicarlo) – fornisca risultati identici a quelli di un metodo probabilistico.
Ma c'è una ragione sottile e insieme profonda. Quando l'informazione in ingresso è scarsa, perché i dati sono pochi, e ciascun dato è affetto da elevata incertezza di misura, il risultato è determinato in larga parte dagli assunti alla base del modello matematico-statistico impiegato; tanto maggiore è l'informazione contenuta nei dati, perché i dati sono molti, e sono affetti da incertezza di misura ridotta, tanto meno influiscono gli assunti sui risultati, che 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 le componenti fondamentali dell'informazione contenuta nei dati e che la qualità dei risultati in uscita non può essere migliore della qualità dei dati in ingresso.
----------
[1] 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
[2] Nouvelles méthodes pour la détermination des orbites des comètes, par A.-M.
Legendre ...
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..."
Nessun commento:
Posta un commento