63  Modello bayesiano di regressione lineare bivariata

In questo capitolo imparerai a:

• Comprendere il modello di regressione bayesiano e come si differenzia dall’approccio frequentista.
• Interpretare i parametri stimati in un contesto bayesiano e confrontarli con quelli frequentisti.
• Familiarizzare con l’uso di brms nella regressione.
• Interpretare le previsioni del modello bayesiano e le verifiche predittive a posteriori.

Prerequisiti

• Leggere l’Appendice L.
• Leggere il capitolo Simple Normal Regression di Johnson et al. (2022).
• Consultare Regression and Other Stories (Gelman et al., 2021).
  • Il pdf del libro è consultabile gratuitamente su questo sito.
  • Prestare particolare attenzione ai capitoli 1 “Overeview, 6,”Background on Regression Modeling,” 7, “Linear Regression with a Single Predictor” e 8, “Fitting regression models”, che offrono una guida dettagliata al modello di regressione bivariato da una prospettiva bayesiana.
• Per utilizzare il pacchetto R brms, è necessario installare preliminarmente Stan o CmdStan sul proprio computer. Si consiglia di optare per CmdStan. Il metodo più semplice per installare CmdStan consiste nell’installare il pacchetto R cmdstanr e seguire le istruzioni fornite nella documentazione.

Preparazione del Notebook

here::here("code", "_common.R") |> 
  source()

# Load packages
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(cmdstanr, posterior, see, brms)

63.1 Introduzione

In questa sezione della dispensa esploreremo il modello di regressione lineare bivariata (cioè con una sola variabile indipendente), ponendo particolare attenzione alla formulazione bayesiana. Confronteremo questo approccio con quello frequentista, così da mettere in luce i principali vantaggi dell’inferenza bayesiana, senza introdurre tecnicismi inutilmente complessi. Il nostro obiettivo comprendere come si costruisce un modello di regressione e come si interpretano i risultati nel contesto dell’incertezza.

63.2 Il modello di regressione lineare semplice

Supponiamo di voler prevedere una variabile quantitativa \(y\) (per esempio, il livello di ansia) a partire da una variabile esplicativa \(x\) (per esempio, il numero di ore di sonno). Il modello di regressione lineare assume che la relazione tra le due variabili sia descritta da:

\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i \]

dove:

• \(\beta_0\) è l’intercetta (valore medio di \(y\) quando \(x = 0\)),
• \(\beta_1\) è il coefficiente di regressione (quanto cambia \(y\) per ogni unità di aumento di \(x\)),
• \(\varepsilon_i\) è un termine di errore casuale, che tiene conto delle deviazioni impreviste rispetto alla linea di regressione.

Assumiamo che gli errori \(\varepsilon_i\) siano indipendenti e distribuiti normalmente con media zero e varianza costante \(\sigma^2\). Questo implica che anche i valori di \(y_i\), condizionati ai rispettivi \(x_i\), seguano una distribuzione normale:

\[ y_i \sim \mathcal{N}(\mu_i, \sigma^2), \quad \text{con} \quad \mu_i = \beta_0 + \beta_1 x_i . \]

63.3 Cos’è la verosimiglianza?

La verosimiglianza è una funzione che misura quanto bene un certo insieme di parametri del modello (cioè \(\beta_0\), \(\beta_1\), \(\sigma\)) riesce a spiegare i dati osservati \((x_i, y_i)\). Si tratta di un concetto centrale sia nell’inferenza frequentista sia in quella bayesiana.

Nel nostro caso, la probabilità di osservare un dato \(y_i\) condizionato a \(x_i\) è data dalla densità di una normale. Assumendo che tutte le osservazioni siano indipendenti, la funzione di verosimiglianza congiunta si ottiene moltiplicando le densità per ciascuna osservazione:

\[ \mathcal{L}(\beta_0, \beta_1, \sigma \mid \mathbf{y}, \mathbf{x}) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{(y_i - \beta_0 - \beta_1 x_i)^2}{2\sigma^2}\right) . \]

Poiché lavorare con prodotti può essere complicato, spesso si usa la log-verosimiglianza, che trasforma i prodotti in somme:

\[ \log \mathcal{L}(\beta_0, \beta_1, \sigma \mid \mathbf{y}, \mathbf{x}) = -\frac{n}{2} \log(2\pi) - n \log \sigma - \frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2 . \]

63.4 Verosimiglianza: confronto tra approccio frequentista e bayesiano

• Nell’approccio frequentista, l’obiettivo è trovare i valori dei parametri che massimizzano la verosimiglianza: è il metodo della massima verosimiglianza.
• Nell’approccio bayesiano, la verosimiglianza non è sufficiente: viene combinata con una distribuzione a priori sui parametri, e il risultato è una distribuzione a posteriori che riflette sia le evidenze empiriche, sia le nostre ipotesi iniziali.

63.5 Le distribuzioni a priori

Uno degli aspetti distintivi del metodo bayesiano è la specificazione delle distribuzioni a priori, che riflettono ciò che si crede possibile per i parametri prima di osservare i dati. Esistono diversi tipi di priori:

• non informativi, che riflettono un’assenza di conoscenza iniziale. Sono distribuzioni molto larghe o piatte, progettate per influenzare il meno possibile la distribuzione a posteriori. Un esempio è una distribuzione normale con varianza molto elevata, come \(\mathcal{N}(0, 1000)\);

• debolmente informativi, che introducono un minimo di informazione strutturale per evitare stime estreme o non plausibili. Sono particolarmente utili nei modelli complessi o quando il numero di dati è limitato. Un esempio comune è \(\mathcal{N}(0, 2.5)\) per i coefficienti di regressione: consente ampio margine di variazione, ma previene valori irrealistici;

• informativi, che riflettono una conoscenza specifica accumulata prima della raccolta dei dati, per esempio da ricerche precedenti, meta-analisi, o teorie ben consolidate. Questi priori sono più ristretti e centrati su valori considerati plausibili. Possono migliorare l’efficienza dell’inferenza, ma devono essere giustificati accuratamente per evitare di introdurre distorsioni arbitrarie.

63.6 Le distribuzioni a posteriori

Una volta specificata:

1. la verosimiglianza (basata sul modello e sui dati), e
2. le distribuzioni a priori (che esprimono le nostre credenze iniziali),

è possibile applicare il teorema di Bayes per ottenere la distribuzione a posteriori di ciascun parametro. Questa distribuzione rappresenta la nostra conoscenza aggiornata dopo aver osservato i dati.

A differenza dell’approccio frequentista, che fornisce stime puntuali, l’approccio bayesiano restituisce intere distribuzioni, permettendo di valutare l’incertezza nei risultati (ad esempio, attraverso gli intervalli di credibilità).

63.7 Implementazione con brms: la formula di Wilkinson

Nel pacchetto brms (che usa Stan come motore di calcolo), non è necessario scrivere la funzione di verosimiglianza a mano. Basta specificare il modello nella classica forma di Wilkinson:

brm(y ~ x, data = dati)

Questa notazione compatta definisce:

• il modello di regressione,
• la verosimiglianza implicita,
• e consente a brms di costruire automaticamente il modello bayesiano.

Se non si specificano i priori, brms utilizza prior debolmente informativi di default.

63.8 Come vengono stimate le distribuzioni a posteriori?

Poiché la distribuzione a posteriori è spesso troppo complessa per essere calcolata esattamente, brms utilizza tecniche di campionamento numerico chiamate MCMC (Markov Chain Monte Carlo).

In particolare, utilizza l’algoritmo NUTS (No-U-Turn Sampler), una variante evoluta dell’algoritmo di Metropolis-Hastings, che esplora lo spazio dei parametri in modo efficiente e adattivo. Grazie a questo, otteniamo campioni dalla distribuzione a posteriori, dai quali è possibile calcolare medie, intervalli di credibilità e fare previsioni.

In sintesi, il modello di regressione bayesiano consente di incorporare in modo trasparente incertezze, conoscenze pregresse e informazioni contenute nei dati. Rispetto all’approccio classico, non restituisce una singola stima puntuale ma un’intera distribuzione per ogni parametro. Questo permette inferenze più flessibili e più ricche di informazioni, particolarmente utili nelle scienze psicologiche, dove l’incertezza è la regola più che l’eccezione.

63.9 Un Esempio Concreto

Definiamo i parametri e simuliamo i dati.

set.seed(123)

# Definizione delle variabili
x <- 1:100
n <- length(x)
a <- 1.5
b <- 0.5
sigma <- 10

# Generazione di y
y <- a + b * x + rnorm(n, 0, sigma)

# Creazione del dataframe
fake <- tibble(x = x, y = y)
head(fake)
#> # A tibble: 6 × 2
#>       x      y
#>   <int>  <dbl>
#> 1     1 -3.60 
#> 2     2  0.198
#> 3     3 18.6  
#> 4     4  4.21 
#> 5     5  5.29 
#> 6     6 21.7

Iniziamo adattando ai dati un modello frequentista:

fm1 <- lm(y ~ x, data = fake)

summary(fm1)
#> 
#> Call:
#> lm(formula = y ~ x, data = fake)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -24.536  -5.524  -0.346   6.485  20.949 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)   1.1360     1.8429    0.62     0.54
#> x             0.5251     0.0317   16.57   <2e-16
#> 
#> Residual standard error: 9.15 on 98 degrees of freedom
#> Multiple R-squared:  0.737,  Adjusted R-squared:  0.734 
#> F-statistic:  275 on 1 and 98 DF,  p-value: <2e-16

Per ottenere l’intervallo di confidenza (nel senso frequentista) della stima dei parametri usiamo:

confint(fm1, level = 0.95)
#>               2.5 % 97.5 %
#> (Intercept) -2.5212  4.793
#> x            0.4622  0.588

Adattiamo ora ai dati un modello di regressione bayesiano utilizzando brms. Si noti che, anche in questo caso, usiamo la sintassi di Wilkinson y ~ x, come per lm(). Eseguiamo il campionamento:

fm2 <- brm(
  y ~ x, 
  data = fake,
  backend = "cmdstanr"
)

Come discusso nell’analisi dell’algoritmo di Metropolis, il primo passo è esaminare le tracce dei parametri per verificare la convergenza dell’algoritmo. La convergenza può essere considerata raggiunta se le catene (nel caso di brm, sono 4 per impostazione predefinita) risultano ben mescolate. Questo si manifesta in un trace plot che mostra una distribuzione uniforme e casuale dei campioni attorno a un valore centrale, senza pattern evidenti o tendenze sistematiche.

Le tracce dei parametri si ottengono nel modo seguente:

mcmc_trace(
  fm2, 
  pars = c("b_Intercept", "b_x", "sigma"),
  facet_args = list(nrow = 3)
)

Gli istogrammi delle distribuzioni a posteriori dei parametri si generano nel modo seguente:

mcmc_hist(
  fm2, 
  pars =c("b_Intercept", "b_x", "sigma"),
  facet_args = list(nrow = 3)
)

Per valutare l’autocorrelazione tra i campioni a posteriori del parametro beta, possiamo utilizzare il seguente comando:

mcmc_acf(fm2, "b_x")

L’autocorrelazione fornisce informazioni sulla dipendenza tra campioni successivi nella catena di Markov. È normale che i campioni successivi non siano completamente indipendenti, poiché le catene di Markov generano campioni correlati per costruzione. Tuttavia, se l’algoritmo ha raggiunto la convergenza, l’autocorrelazione dovrebbe diminuire rapidamente e diventare trascurabile dopo un numero relativamente piccolo di lag. Questo significa che, dopo un certo numero di passi, i campioni diventano progressivamente meno correlati tra loro, comportandosi in modo simile a campioni indipendenti estratti dalla distribuzione target.

Un’elevata autocorrelazione su lag più lunghi potrebbe invece indicare problemi di mescolamento delle catene o una mancata convergenza, richiedendo ulteriori verifiche o aggiustamenti, come l’aumento del numero di iterazioni o una diversa parametrizzazione del modello.

Nel caso presente, notiamo una rapida diminuzione dell’autocorrelazione in funzione del numero di passi. Ciò è indicativo del fatto che la convergenza è stata raggiunta.

Una sintesi numerica dei risultati si trova nel modo seguente:

summary(fm2)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: y ~ x 
#>    Data: fake (Number of observations: 100) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept     1.11      1.85    -2.55     4.77 1.00     3976     2968
#> x             0.53      0.03     0.46     0.59 1.00     4028     3041
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     9.26      0.69     7.99    10.75 1.00     3753     3097
#> 
#> Draws were sampled using sample(hmc). 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).

Confrontiamo le stime ottenute con i valori reali dei parametri simulati. L’intercetta è stata stimata attorno a 1.14, con un’incertezza al 95% che varia tra -2.4 e 4.8. Questo risultato rientra negli intervalli di credibilità previsti, confermando l’accuratezza del modello. Analogamente, per la pendenza \(b\), l’intervallo di credibilità al 95% include il valore reale simulato, dimostrando come le stime bayesiane riflettano accuratamente l’incertezza sui parametri.

Se si utilizza la funzione conditional_effects() viene prodotto un grafico che rappresenta la relazione stimata tra il predittore \(x\) e la variabile di risposta \(y\).

conditional_effects(fm2) |>
  plot(points = TRUE)

1. Linea stimata (effetto medio):
   • La linea centrale del grafico rappresenta il valore medio previsto di \(y\) per ogni valore di \(x\), dato dalla relazione \(y = \alpha + \beta x\).
   • Questa linea è calcolata usando i valori medi a posteriori stimati per \(\alpha\) e \(\beta\).
2. Bande di incertezza (intervalli di credibilità):
   • Le bande attorno alla linea rappresentano gli intervalli di credibilità (ad esempio, al 95%). Questi mostrano l’incertezza associata alle stime del modello per ogni valore di \(x\).
   • Più strette sono le bande, maggiore è la certezza del modello riguardo alla relazione stimata.
3. Dati osservati:
   • I punti rappresentano i valori effettivi di \(y\) osservati nei dati. Questo consente di confrontare visivamente come i dati reali si allineano con le previsioni del modello.

Il grafico consente

• una verifica visiva della relazione stimata tra \(y\) e \(x\);
• di identificazione di eventuali discrepanze tra i dati osservati e le previsioni del modello;
• una rappresentazione dell’incertezza nelle stime.

Ad esempio, il grafico può mostrare se \(x\) ha un effetto credibile su \(y\) e con quale livello di incertezza. Se l’effetto di \(x\) è debole o nullo, la linea stimata sarà piatta (vicina a zero) e le bande di incertezza saranno ampie.

63.10 Simulazione di Livelli di Copertura

Verifichiamo la copertura degli intervalli di credibilità al 95% attraverso simulazioni ripetute.

set.seed(42)
# Parametri veri
a_true <- 0.2
b_true <- 0.3
sigma_true <- 0.5
# Numero di simulazioni
num_simulations <- 100
# Conteggio delle coperture
coverage_a <- 0
coverage_b <- 0
for (i in 1:num_simulations) {
  # Generazione dei dati
  x <- 1:20
  y <- a_true + b_true * x + sigma_true * rnorm(length(x))
  # Adattamento del modello
  fit <- lm(y ~ x)
  ci <- confint(fit) # Intervalli di confidenza
  # Verifica delle coperture
  if (ci[1,1] <= a_true & ci[1, 2] >= a_true) {
    coverage_a <- coverage_a + 1
  }
  if (ci[2,1] <= b_true & ci[2, 2] >= b_true) {
    coverage_b <- coverage_b + 1
  }
}

# Risultati
cat("Coverage for a:", coverage_a / num_simulations, "\n")
#> Coverage for a: 0.93
cat("Coverage for b:", coverage_b / num_simulations, "\n")
#> Coverage for b: 0.96

I risultati indicano che i livelli di copertura empirici ottenuti con l’approccio frequentista corrispondono strettamente ai livelli teorici attesi.

Per proseguire, ripeteremo la simulazione adottando un approccio bayesiano. Useremo la funzione brm() del pacchetto brms al posto di lm().

#| message: false
#| warning: false
#| output: false
#| 
set.seed(23)
n_fake   <- 100
cover_68 <- logical(n_fake)
cover_95 <- logical(n_fake)

# Veri parametri
a     <- 0.2    # intercetta vera
b     <- 0.3    # pendenza vera
sigma <- 0.5    # deviazione standard vera
x     <- 1:20
n     <- length(x)

# Priors con set_prior 
priors <- c(
  set_prior("normal(0, 2.5)", class = "Intercept"),
  set_prior("normal(0, 2.5)", class = "b", coef = "x"),
  set_prior("cauchy(0, 2.5)", class = "sigma")
)

set.seed(23)
n_fake   <- 1000
cover_68 <- logical(n_fake)
cover_95 <- logical(n_fake)

a     <- 0.2
b     <- 0.3
sigma <- 0.5
x     <- 1:20
n     <- length(x)

for (s in seq_len(n_fake)) {
  y    <- a + b * x + rnorm(n, 0, sigma)
  fake <- data.frame(x = x, y = y)

  fit <- brm(
    y ~ 1 + x,
    data    = fake,
    family  = gaussian(),
    prior   = priors,
    iter    = 2000,
    chains  = 2,
    refresh = 0,
    backend = "cmdstanr"
  )

  post <- summary(fit)$fixed
  b_hat <- post["x", "Estimate"]
  b_se  <- post["x", "Est.Error"]

  cover_68[s] <- abs(b - b_hat) < b_se
  cover_95[s] <- abs(b - b_hat) < 2 * b_se
}

cat("Coverage 68%:", mean(cover_68), "\n")
cat("Coverage 95%:", mean(cover_95), "\n")

Con solo 100 iterazioni, i risultati sono i seguenti:

> cat("Coverage 68%:", mean(cover_68), "\n")
Coverage 68%: 0.73 
> cat("Coverage 95%:", mean(cover_95), "\n")
Coverage 95%: 0.953 

Questa seconda simulazione evidenzia che anche i livelli di copertura empirici ottenuti con l’approccio bayesiano si avvicinano ai valori teorici previsti.

I risultati ottenuti confermano l’efficacia degli intervalli di confidenza e di credibilità stimati attraverso i modelli frequentisti e bayesiani.

63.11 Confronti, non Effetti

Gelman et al. (2021) sottolineano che i coefficienti di regressione sono spesso denominati “effetti”, ma questa terminologia può trarre in inganno. Gli “effetti”, infatti, implicano una relazione causale. Tuttavia, ciò che un modello di regressione stima non è necessariamente un effetto causale, ma piuttosto un pattern osservazionale. In particolare, ciò che osserviamo è che la media della variabile dipendente nella sottopopolazione con \(X = x + 1\) è spesso maggiore o minore (a seconda del segno di \(\beta\)) rispetto alla media della sottopopolazione con \(X = x\).

La regressione è uno strumento matematico utilizzato principalmente per fare previsioni. I coefficienti di regressione devono quindi essere interpretati come confronti medi. Solo in circostanze specifiche, quando la regressione descrive un processo causale ben definito, è possibile interpretarli come effetti. Tuttavia, questa interpretazione causale deve essere giustificata dal disegno dello studio e non può essere dedotta unicamente dall’uso del modello statistico.

63.12 Riflessioni Conclusive

In questo capitolo abbiamo esplorato il modello di regressione lineare bivariata adottando una prospettiva bayesiana. Abbiamo visto che, quando si utilizzano priori debolmente informativi, le stime bayesiane tendono a coincidere con quelle dell’approccio frequentista, specialmente in presenza di un numero consistente di dati. Tuttavia, il valore distintivo dell’approccio bayesiano non si limita alla stima dei parametri: esso risiede soprattutto nella possibilità di integrare conoscenze a priori e di rappresentare l’incertezza in modo esplicito e probabilistico.

Indipendentemente dalla scelta tra approccio bayesiano o frequentista, è importante sottolineare che i modelli statistici non offrono verità definitive. Come ricordato da Alexander (2023), i modelli non sono specchi della realtà, ma strumenti concettuali per darle significato. Agiscono come lenti interpretative, attraverso cui possiamo mettere a fuoco specifici aspetti del fenomeno studiato.

In particolare, i modelli statistici possono essere utilizzati per due scopi principali, entrambi fondamentali nella ricerca psicologica:

• Previsione: si concentra sulla capacità del modello di anticipare nuovi dati. È un uso pragmatico ed empirico della modellizzazione, in cui l’attenzione è rivolta alla bontà predittiva del modello, spesso valutata con tecniche di validazione incrociata.

• Inferenza: mira a comprendere le relazioni causali tra variabili. Per rendere credibili le inferenze causali, la regressione deve essere accompagnata da una solida progettazione dello studio (es. esperimenti, disegni longitudinali, controllo di variabili confondenti) e da ipotesi teoriche ben motivate.

La regressione lineare, in ogni caso, rappresenta una forma di media ponderata tra le osservazioni, il che implica alcune limitazioni:

• i risultati possono essere distorti dalla presenza di variabili confondenti non incluse nel modello;
• la qualità dei dati e la verifica delle ipotesi del modello (normalità degli errori, indipendenza, omoschedasticità) giocano un ruolo cruciale nella validità delle stime;
• l’assunzione di linearità potrebbe non essere adatta a descrivere alcune relazioni psicologiche, che spesso sono complesse o non lineari.

Pertanto, è essenziale non assumere il modello come un fine, ma considerarlo uno strumento di esplorazione e interpretazione, da integrare in una riflessione teorica più ampia. I risultati di una regressione non parlano da soli: vanno interpretati alla luce del disegno dello studio, della letteratura scientifica, e delle ipotesi formulate dal ricercatore.

Per chi desidera approfondire il modello bayesiano di regressione lineare, oltre al testo di Johnson et al. (2022), si consiglia la lettura del testo Regression and Other Stories. Quest’ultimo rappresenta una guida pratica e ricca di esempi applicati, ideale per comprendere come integrare i metodi bayesiani nell’analisi di regressione e interpretarne i risultati in contesti reali.

Problemi

1. Verosimiglianza
   Definire la funzione di verosimiglianza per il modello bayesiano di regressione lineare bivariata, esplicitando i parametri e la loro interpretazione.

2. Scelta dei prior
   Proporre un set di prior debolmente informativi differenti da quelli riportati, motivando la scelta delle distribuzioni.

3. Simulazione dati
   Scrivere un blocco di codice in R/Quarto per simulare un dataset con \(n=50\), parametri \(lpha=2\), \(eta=0.5\), \(\sigma=1\) e visualizzare un grafico dispersione con la retta di regressione vera.

4. Stima frequentista vs bayesiana
   Utilizzando i dati simulati, adattare un modello con lm() e uno con brm() (specificando i prior). Confrontare i risultati prodotti dai due approcci, riportando i valori stimati e gli intervalli di confidenza/credibilità.

5. Diagnosi MCMC
   Elencare e spiegare almeno tre controlli diagnostici da effettuare sulle catene MCMC per garantire la convergenza e un buon mescolamento.

6. Interpretazione dei coefficienti
   In un contesto osservazionale, discutere perché non è appropriato interpretare i coefficienti della regressione come effetti causali. Fare riferimento ai concetti di confondimento e disegno sperimentale.

Soluzioni

1. Verosimiglianza

   \[ \mathcal{L}(\alpha, \beta, \sigma \mid y, x) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{(y_i - (\alpha + \beta x_i))^2}{2\sigma^2}\right). \]

2. Scelta dei prior

   Ad esempio:

   • \(\alpha \sim \mathcal{N}(0, 5)\): consente maggiore variabilità iniziale;
     
   • \(\beta \sim \text{Student-}t(3, 0, 2)\): robuste alle code pesanti;
     
   • \(\sigma \sim \text{Half-}Cauchy(0, 1)\): prior leggermente più stretto sulle deviazioni.

3. Simulazione dati

   set.seed(42)
   n <- 50
   alpha <- 2
   beta  <- 0.5
   sigma <- 1
   x <- rnorm(n, 0, 1)
   y <- alpha + beta * x + rnorm(n, 0, sigma)
   plot(x, y, main = "Dati simulati", xlab = "x", ylab = "y")
   abline(a = alpha, b = beta, col = "blue", lwd = 2)

4. Stima frequentista vs bayesiana

   • Frequentista (lm()):

     fm_f <- lm(y ~ x)
     summary(fm_f)
     confint(fm_f, level = 0.95)

   • Bayesiano (brm()):

     fm_b <- brm(y ~ x, data = data.frame(x, y),
                 prior = c(set_prior("normal(0,5)", class="Intercept"),
                           set_prior("normal(0,5)", class="b"),
                           set_prior("cauchy(0,1)", class="sigma")),
                 iter = 2000, chains = 2)
     summary(fm_b)

     Confronto: i valori medi a posteriori e gli intervalli di credibilità dovrebbero includere quelli di confidenza di lm().

5. Diagnosi MCMC

   • Trace plots: verificare mescolamento e stazionarietà delle catene;
     
   • R-hat: valore vicino a 1 indica convergenza;
     
   • Effective Sample Size (ESS): numero di campioni indipendenti effettivi.

6. Interpretazione dei coefficienti
   La regressione osservazionale non controlla automaticamente i potenziali confondenti; senza randomizzazione o disegno sperimentale rigoroso, non si può inferire causalità. Bisogna considerare variabili confondenti e criteri di validità interna.

Informazioni sull’Ambiente di Sviluppo

sessionInfo()
#> R version 4.5.0 (2025-04-11)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.4.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/UTF-8/C/C/C/C
#> 
#> time zone: Europe/Rome
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] brms_2.22.0      Rcpp_1.0.14      posterior_1.6.1  cmdstanr_0.9.0  
#>  [5] thematic_0.1.6   MetBrewer_0.2.0  ggokabeito_0.1.0 see_0.11.0      
#>  [9] gridExtra_2.3    patchwork_1.3.0  bayesplot_1.12.0 psych_2.5.3     
#> [13] scales_1.4.0     markdown_2.0     knitr_1.50       lubridate_1.9.4 
#> [17] forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4      purrr_1.0.4     
#> [21] readr_2.1.5      tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.2   
#> [25] tidyverse_2.0.0  rio_1.2.3        here_1.0.1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1     farver_2.1.2         loo_2.8.0           
#>  [4] fastmap_1.2.0        tensorA_0.36.2.1     pacman_0.5.1        
#>  [7] digest_0.6.37        timechange_0.3.0     lifecycle_1.0.4     
#> [10] StanHeaders_2.32.10  processx_3.8.6       magrittr_2.0.3      
#> [13] compiler_4.5.0       rlang_1.1.6          tools_4.5.0         
#> [16] yaml_2.3.10          data.table_1.17.0    labeling_0.4.3      
#> [19] bridgesampling_1.1-2 htmlwidgets_1.6.4    pkgbuild_1.4.7      
#> [22] mnormt_2.1.1         curl_6.2.2           plyr_1.8.9          
#> [25] RColorBrewer_1.1-3   abind_1.4-8          withr_3.0.2         
#> [28] grid_4.5.0           stats4_4.5.0         colorspace_2.1-1    
#> [31] inline_0.3.21        cli_3.6.5            mvtnorm_1.3-3       
#> [34] rmarkdown_2.29       generics_0.1.3       RcppParallel_5.1.10 
#> [37] rstudioapi_0.17.1    reshape2_1.4.4       tzdb_0.5.0          
#> [40] rstan_2.32.7         parallel_4.5.0       matrixStats_1.5.0   
#> [43] vctrs_0.6.5          V8_6.0.3             Matrix_1.7-3        
#> [46] jsonlite_2.0.0       hms_1.1.3            glue_1.8.0          
#> [49] codetools_0.2-20     ps_1.9.1             distributional_0.5.0
#> [52] stringi_1.8.7        gtable_0.3.6         QuickJSR_1.7.0      
#> [55] pillar_1.10.2        htmltools_0.5.8.1    Brobdingnag_1.2-9   
#> [58] R6_2.6.1             rprojroot_2.0.4      evaluate_1.0.3      
#> [61] lattice_0.22-7       backports_1.5.0      rstantools_2.4.0    
#> [64] coda_0.19-4.1        nlme_3.1-168         checkmate_2.3.2     
#> [67] xfun_0.52            pkgconfig_2.0.3

Bibliografia

Alexander, R. (2023). Telling Stories with Data: With Applications in R. Chapman; Hall/CRC.

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.

Johnson, A. A., Ott, M., & Dogucu, M. (2022). Bayes Rules! An Introduction to Bayesian Modeling with R. CRC Press.