63  Zucchero sintattico

In questo capitolo imparerai a
  • utilizzare brm per costruire e adattare modelli;
  • analizzare i risultati con brm.
Prerequisiti
Preparazione del Notebook 
here::here("code", "_common.R") |> 
  source()

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

63.1 Introduzione

I modelli lineari sono così ampiamente utilizzati che sono stati sviluppati appositamente una sintassi, dei metodi e delle librerie per la regressione. Una di queste librerie è brms (Bayesian Regression Models using Stan), già introdotta nel Capitolo 62. brms è un pacchetto R progettato per adattare modelli gerarchici generalizzati lineari (di cui il modello lineare bivariato è un caso particolare), utilizzando una sintassi simile a quella presente nei pacchetti R, come lme4, nlme, rstanarm. brms si basa su Stan, ma offre un’API di livello superiore.

In questo capitolo esploreremo in maniera dettagliata come condurre un’analisi di regressione utilizzando brms invece di Stan.

63.2 Interfaccia brms

Per fare un esempio, applicheremo il modello di regressione bivariato alla relazione tra altezza e peso. I dati contenuti nel file Howell_18.csv sono parte di un censimento parziale della popolazione !Kung San dell’area di Dobe, raccolti tramite interviste condotte da Nancy Howell alla fine degli anni ’60 (McElreath, 2020). I !Kung San sono una delle popolazioni di raccoglitori-cacciatori più conosciute del ventesimo secolo e sono stati oggetto di numerosi studi antropologici. In questa analisi, consideriamo un sottocampione di dati relativi alla popolazione adulta (di età superiore ai 18 anni).

Importiamo i dati contenuti nel file Howell_18.csv.

df <- rio::import(here::here("data", "Howell_18.csv"))
df |> 
  head()
#>   height weight age male
#> 1  151.8  47.83  63    1
#> 2  139.7  36.49  63    0
#> 3  136.5  31.86  65    0
#> 4  156.8  53.04  41    1
#> 5  145.4  41.28  51    0
#> 6  163.8  62.99  35    1

Generiamo un diagramma a dispersione tra le variabili height (altezza) e weight (peso):

ggplot(df, aes(x = weight, y = height)) +
  geom_point() +  
  labs(x = "Weight", y = "Height")

brms si concentra sui modelli di regressione, e questa specializzazione permette di adottare una sintassi più semplice, conosciuta come sintassi di Wilkinson (Wilkinson & Rogers, 1973).

Ad esempio, il modello \(y = \alpha + \beta x + \varepsilon\) si implementa come segue:

a_model = brm(y ∼ 1 + x, data = df)

Nella sintassi di Wilkinson, il simbolo tilde (∼) separa la variabile dipendente (a sinistra) dalle variabili indipendenti (a destra). In questo caso, stiamo specificando solo la media (\(\mu\)) della \(y\).

brms assume di default che la distribuzione di verosimiglianza sia gaussiana, ma è possibile modificarla tramite l’argomento family.

La notazione 1 si riferisce all’intercetta. L’intercetta viene inclusa di default. Per cui il modello precedente si può anche scrivere, in maniera equivalente, come

a_model = brm(y ∼ x, data = df)

Se desideriaamo escludere l’intercetta dal modello, possiamo farlo in questo modo

no_intercept_model = brm(y ∼ 0 + x, data = df)

oppure in questo modo

no_intercept_model = brm(y ∼ -1 + x, data = df)

Per includere ulteriori variabili nel modello, possiamo procedere così:

model_2 = brm("y ∼ x + z", data)

brms consente anche di includere effetti a livello di gruppo (gerarchici). Ad esempio, se desideriamo un modello ad effetti misti nel quale abbiamo un effetto diverso di \(x\) in ciascun gruppo g, possiamo usare la seguente sintassi:

model_h = brm(y ∼ x + z + (x | g), data = df)

La sintassi di Wilkinson non specifica le distribuzioni a priori, ma solo come le variabili dipendenti e indipendenti sono collegate. brms definirà automaticamente delle distribuzioni a priori debolmente informative per noi, rendendo superflua la loro definizione esplicita. Tuttavia, se preferiamo avere un maggiore controllo, possiamo specificarle manualmente, come vedremo in seguito.

63.2.1 Centrare le Variabili

Per interpretare più facilmente l’intercetta, centriamo la variabile weight rispetto alla media del campione:

df$weight_c <- df$weight - mean(df$weight)

Ora, l’intercetta (\(\alpha\)) rappresenterà l’altezza media quando il peso corrisponde alla media del campione.

Adattiamo un modello lineare con la variabile weight centrata e esaminiamo i risultati:

fit_1 = brm(
  bf(height ~ 1 + weight_c, center = FALSE), 
  data = df, 
  backend = "cmdstanr", 
  silent = 0
)

Utilizzando center = FALSE nel modello Bayesiano garantiamo che il centraggio venga mantenuto e non applicato nuovamente da brms.

Le tracce dei parametri si ottengono nel modo seguente:

mcmc_trace(
  fit_1, 
  pars = c("b_Intercept", "b_weight_c", "sigma"),
  facet_args = list(nrow = 3)
)

summary(fit_1)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: height ~ 1 + weight_c 
#>    Data: df (Number of observations: 352) 
#>   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   154.59      0.27   154.06   155.11 1.00     3947     2871
#> weight_c      0.91      0.04     0.82     0.99 1.00     3945     2996
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     5.10      0.20     4.74     5.51 1.00     4902     3275
#> 
#> 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).

La stima dell’intercetta \(\alpha\) = 154.60 suggerisce che, per le persone con un peso corrispondente al valore medio del campione analizzato, l’altezza prevista è di 154.60 cm.

Possiamo confrontare i risultati ottenuti da brm() con quelli prodotti dall’appriccio frequentista:

fit_2 <- lm(height ~ 1 + weight_c, data = df)
summary(fit_2)
#> 
#> Call:
#> lm(formula = height ~ 1 + weight_c, data = df)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -19.746  -2.884   0.022   3.142  14.774 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)  154.597      0.271   570.2   <2e-16
#> weight_c       0.905      0.042    21.5   <2e-16
#> 
#> Residual standard error: 5.09 on 350 degrees of freedom
#> Multiple R-squared:  0.57,   Adjusted R-squared:  0.568 
#> F-statistic:  463 on 1 and 350 DF,  p-value: <2e-16

Anche in questo caso, l’uso di prior debolmente informativi fa in modo che i risultati dei due approcci siano praticamente equivalenti.

63.3 Visualizzazione dei Risultati

Per comprendere visivamente la relazione stimata tra peso e altezza nel nostro modello bayesiano, utilizziamo la funzione conditional_effects:

conditional_effects(fit_1, effects = "weight_c")

Il grafico generato fornisce una rappresentazione completa della relazione stimata:

  • Linea centrale (media posteriore): Rappresenta la stima più probabile dell’altezza per ciascun valore del peso centrato.
  • Area colorata (intervallo di credibilità): Mostra l’intervallo di densità più alta (HDI) al 95%, indicando l’incertezza attorno alla stima centrale.
  • Punti sul grafico: Rappresentano i dati reali osservati, permettendo di confrontare il modello con i dati originali.

È possibile adattare il livello di incertezza mostrato modificando l’argomento prob:

# Visualizzazione con intervallo di credibilità all'89%
conditional_effects(fit_1, effects = "weight_c", prob = 0.89)

Ridurre il valore di prob (es. a 0.80 o 0.50) produce un intervallo più stretto, mentre aumentarlo (es. a 0.99) amplia l’area di incertezza visualizzata.

63.3.1 Interpretazione Pratica del Grafico

Nel grafico:

  • Il punto dove la linea attraversa weight_c = 0 corrisponde all’altezza prevista per un individuo con peso medio (poiché abbiamo centrato la variabile).
  • La pendenza della linea indica quanto ci aspettiamo che l’altezza aumenti per ogni kg di peso aggiuntivo.
  • La larghezza dell’intervallo di credibilità indica la nostra certezza sulle stime: più stretto l’intervallo, maggiore la certezza.

63.4 Due Tipi di Incertezza nei Modelli Bayesiani

Quando visualizziamo i risultati di un modello bayesiano, esistono due tipi fondamentali di incertezza che possiamo rappresentare:

  1. Incertezza del parametro (metodo predefinito): Mostra quanto siamo incerti riguardo alla “vera” relazione tra variabili.
  2. Incertezza predittiva (method = "predict"): Mostra quanto siamo incerti riguardo a singole osservazioni future.

La visualizzazione che incorpora l’incertezza sui parametri si ottiene con:

conditional_effects(fit_1, effects = "weight_c")

Questa visualizzazione standard mostra:

  • La linea di regressione stimata (media posteriore)
  • L’intervallo di credibilità attorno alla linea (solitamente al 95%)

Ciò che vediamo è l’incertezza della media condizionale - ovvero, quanto siamo incerti sul valore medio dell’altezza per ogni valore del peso.

La visualizzazione predittiva si ottiene con:

conditional_effects(fit_1, effects = "weight_c", method = "predict")

Questa visualizzazione alternativa mostra:

  • La stessa linea di regressione stimata.
  • Un intervallo di predizione molto più ampio.

L’intervallo mostrato rappresenta:

  • L’incertezza dei parametri (come nella visualizzazione standard).
  • PIÙ la variabilità residua del modello (la dispersione naturale dei dati).

63.4.1 Perché sono Diverse e Quale Scegliere?

  • Visualizzazione standard: Risponde alla domanda “Quanto siamo sicuri della relazione media tra peso e altezza?”
  • Visualizzazione predittiva: Risponde alla domanda “Dove ci aspettiamo che cadrà la prossima osservazione individuale?”

63.4.1.1 Differenze Visive Notevoli:

  1. Ampiezza dell’intervallo: L’intervallo predittivo è sempre considerevolmente più ampio dell’intervallo di credibilità per la media.
  2. Interpretazione pratica: L’intervallo predittivo include la variabilità naturale delle osservazioni individuali.

63.4.2 Quando Usare Ciascun Metodo:

  • Usa il metodo predefinito quando sei interessato alla relazione tra variabili.
  • Usa method = "predict" quando vuoi fare previsioni su nuove osservazioni individuali.

63.4.2.1 Esempio Concreto

Immagina di voler predire l’altezza di una persona specifica con un dato peso:

  • La visualizzazione standard ti dirà quanto sei sicuro del valore medio per persone con quel peso.
  • La visualizzazione predittiva ti dirà in quale intervallo probabilmente cadrà un’osservazione individuale, incorporando sia l’incertezza del modello che la variabilità naturale.

In sintesi, la differenza principale sta nell’inclusione o meno della variabilità residua. La visualizzazione predittiva è più onesta riguardo alla nostra capacità di fare previsioni per singole osservazioni, mentre quella standard è più adatta per comprendere la relazione generale tra le variabili.

63.5 Distribuzione a Posteriori dei Parametri

Per esaminare la distribuzione a posteriori dei parametri usiamo la funzione mcmc_plot():

mcmc_plot(fit_1, type = "dens")

Esaminiamo in maggiori dettagli il sommario numerico dei parametri stimati:

draws <- posterior::as_draws(fit_1, variable = "^b_", regex = TRUE)
posterior::summarise_draws(draws, "mean", "sd", "mcse_mean", "mcse_sd")
#> # A tibble: 2 × 5
#>   variable       mean     sd mcse_mean  mcse_sd
#>   <chr>         <dbl>  <dbl>     <dbl>    <dbl>
#> 1 b_Intercept 155.    0.269   0.00428  0.00391 
#> 2 b_weight_c    0.906 0.0423  0.000675 0.000681

Utilizziamo la funzione as_draws() che trasforma un oggetto R in un formato compatibile con posterior. Gli argomenti variable = "^b_" e regex = TRUEconsentono di selezionare solo i parametri il cui nome inizia con b_: nel nostro caso saranno l’intercetta e la pendenza del modello di regressione lineare.

Successivamente usiamo la funzione summarise_draws() con gli argomenti specificati per un sommario della distribuzione a posteriori dei parametri prescelti.

63.5.1 Spiegazione di mcse_mean e mcse_sd

I valori mcse_mean e mcse_sd sono le Monte Carlo Standard Errors (errori standard Monte Carlo) per la stima della media (mean) e della deviazione standard (sd), rispettivamente. Questi valori quantificano l’incertezza associata al processo di campionamento effettuato durante l’analisi bayesiana, in particolare quando si utilizzano algoritmi Monte Carlo come MCMC (Markov Chain Monte Carlo).

63.5.1.1 mcse_mean

  • Rappresenta l’errore standard Monte Carlo per la stima della media.
  • Indica quanto la media stimata (\(\text{mean}\)) potrebbe variare a causa della finitezza dei campioni generati dall’algoritmo MCMC.
  • Un valore di mcse_mean basso rispetto alla deviazione standard (\(\text{sd}\)) suggerisce che il numero di campioni generati è sufficiente per ottenere una stima accurata della media.

63.5.1.2 mcse_sd

  • Rappresenta l’errore standard Monte Carlo per la stima della deviazione standard.
  • Indica quanto potrebbe variare la stima della deviazione standard (\(\text{sd}\)) a causa del numero finito di campioni generati.
  • Anche qui, un valore basso di mcse_sd rispetto alla sd suggerisce che l’incertezza introdotta dal campionamento è trascurabile.

63.5.2 Come interpretarli?

  • Proporzione rispetto alla sd:
    • mcse_mean e mcse_sd dovrebbero essere molto più piccoli rispetto ai rispettivi parametri (mean e sd), idealmente almeno un ordine di grandezza inferiore.
    • Ad esempio, per b_Intercept, mcse_mean = 0.0044 è molto più piccolo rispetto a sd = 0.2695, indicando che la stima della media è robusta.
  • Indicazione della qualità del campionamento:
    • Valori alti di mcse_mean o mcse_sd rispetto alla sd potrebbero indicare che il numero di iterazioni MCMC non è sufficiente, che le catene non sono ben mescolate o che ci sono problemi di convergenza.

In sintesi, mcse_mean e mcse_sd sono utili per valutare l’affidabilità delle stime derivate dal campionamento Monte Carlo. Se questi valori sono bassi, possiamo essere confidenti che il numero di campioni è sufficiente per rappresentare accuratamente la distribuzione a posteriori.

63.6 Specificare i Priors

Se vogliamo personalizzare i priors, possiamo utilizzare la funzione get_prior per esplorare quelli predefiniti:

get_prior(height ~ 1 + weight_c, data = df)
#>                     prior     class     coef group resp dpar nlpar lb ub
#>  student_t(3, 154.3, 8.5) Intercept                                     
#>                    (flat)         b                                     
#>                    (flat)         b weight_c                            
#>      student_t(3, 0, 8.5)     sigma                                 0   
#>        source
#>       default
#>       default
#>  (vectorized)
#>       default
  1. prior: descrive il prior predefinito assegnato a ciascun parametro del modello. Ad esempio:
    • student_t(3, 154.3, 8.5): prior t di Student per l’intercetta con 3 gradi di libertà, una media di 154.3, e una scala di 8.5.
    • (flat): prior piatto (non informativo) per i coefficienti delle variabili predittive, come \(b\) e \(b_{weight_c}\).
    • student_t(3, 0, 8.5): prior t di Student per il parametro \(\sigma\) (deviazione standard residua), centrato su 0 con una scala di 8.5.
  2. class: identifica la classe di parametro a cui il prior si applica:
    • Intercept: prior per l’intercetta (\(\alpha\)).
    • b: prior per i coefficienti delle variabili predittive (\(\beta\)).
    • sigma: prior per il parametro della deviazione standard residua (\(\sigma\)).
  3. coef: specifica a quale predittore si riferisce il prior, se applicabile. Ad esempio:
    • Vuoto per l’intercetta (poiché non dipende da un predittore specifico).
    • weight_c per il coefficiente relativo al predittore weight_c.
  4. lb e ub: rappresentano rispettivamente i limiti inferiori (lower bound) e superiori (upper bound) per il prior, se specificati. Ad esempio:
    • Per sigma, il limite inferiore è \(0\), dato che la deviazione standard non può essere negativa.
  5. source: indica l’origine del prior. Se il prior è predefinito (default), il valore sarà default. Se un prior è specificato manualmente dall’utente, sarà indicato come tale.

Ora impostiamo priors espliciti e adattiamo un nuovo modello:

prior_guassian <-
  prior(normal(160, 10), class = "b", coef = "Intercept") +
  prior(normal(0, 5), class = "b", coef = "weight_c") +
  prior(cauchy(0, 5), class = "sigma")
fit_2 = brm(
  bf(height ~ 1 + weight_c, center = FALSE), 
  prior = prior_guassian,
  data = df, 
  backend = "cmdstanr", 
  silent = 0
)

Otteniamo un sommario numerico dei parametri stimati:

draws <- posterior::as_draws(fit_2, variable = "^b_", regex = TRUE)
posterior::summarise_draws(draws, "mean", "sd", "mcse_mean", "mcse_sd")
#> # A tibble: 2 × 5
#>   variable       mean     sd mcse_mean  mcse_sd
#>   <chr>         <dbl>  <dbl>     <dbl>    <dbl>
#> 1 b_Intercept 155.    0.273   0.00425  0.00409 
#> 2 b_weight_c    0.905 0.0428  0.000703 0.000678

I prior che abbiamo specificato non cambiano in maniera rilevante la soluzione a posteriori.

63.7 Predizioni Predittive a Posteriori

Un aspetto fondamentale nella valutazione di un modello statistico, sia frequentista che bayesiano, è verificare quanto bene i dati osservati siano rappresentati dalle predizioni del modello. Tuttavia, l’approccio e l’interpretazione differiscono tra i due paradigmi.

63.7.1 Confronto Frequentista

Nel caso frequentista, si confrontano i valori predetti dal modello, \(\hat{y} = \hat{\alpha} + \hat{\beta}x\), con i dati osservati. Questo confronto si basa sull’analisi di: - La vicinanza della retta di regressione stimata ai dati osservati. - L’eventuale presenza di pattern nei dati che si discostano da un andamento lineare. - La variazione della dispersione dei valori di \(y\) rispetto a \(x\) (ad esempio, per verificare l’ipotesi di omoschedasticità).

63.7.2 Approccio Bayesiano

Nell’approccio bayesiano, si eseguono le stesse verifiche di base, ma l’analisi si arricchisce attraverso l’uso delle Predizioni Predittive a Posteriori (Posterior Predictive Checks, PPCs). Questo metodo consente di confrontare i dati osservati con dati simulati dal modello, utilizzando l’incertezza stimata nelle distribuzioni a posteriori dei parametri.

63.7.3 Costruzione delle Predizioni Predittive a Posteriori

Nel caso di un modello bivariato, il processo per generare un grafico delle Predizioni Predittive a Posteriori è il seguente:

  1. Dati osservati: Si parte dall’istogramma lisciato dei dati osservati, che rappresenta la distribuzione empirica di \(y\).

  2. Simulazione di dati predetti:

    • Si estrae un campione casuale di valori \(\alpha'\), \(\beta'\), e \(\sigma'\) dalle distribuzioni a posteriori dei parametri (\(\alpha\), \(\beta\), e \(\sigma\)).
    • Usando questi valori, si calcolano dati simulati da una distribuzione normale: \[ y_{\text{sim}} \sim \mathcal{N}(\alpha' + \beta'x, \sigma') \] dove \(x\) sono i valori predittori osservati.

  3. Creazione di istogrammi: Per ogni campione simulato, si costruisce un istogramma lisciato che rappresenta la distribuzione predetta dal modello.

  4. Ripetizione: Il processo viene ripetuto più volte, generando molti istogrammi lisciati.

  5. Confronto: Tutti gli istogrammi predetti vengono sovrapposti all’istogramma dei dati osservati. Questo consente di confrontare visivamente la capacità del modello di rappresentare la distribuzione dei dati.

63.7.4 Interpretazione

  • Buona corrispondenza: Se gli istogrammi lisciati dei dati simulati si sovrappongono bene all’istogramma dei dati osservati, significa che il modello è in grado di rappresentare adeguatamente il campione corrente.
  • Discrepanze: Se vi sono discrepanze sistematiche (ad esempio, picchi o code mancanti nei dati predetti rispetto agli osservati), ciò indica che il modello potrebbe non essere adeguato o che vi sono aspetti dei dati non catturati dal modello.

L’approccio delle Predizioni Predittive a Posteriori è particolarmente potente perché:

  • Integra l’incertezza nei parametri del modello.
  • Permette di verificare non solo la bontà di adattamento complessiva, ma anche specifici aspetti delle distribuzioni predette.
  • È visivo e intuitivo, facilitando l’identificazione di discrepanze tra modello e dati.

In conclusione, le Predizioni Predittive a Posteriori forniscono un modo robusto per valutare l’adeguatezza di un modello bayesiano rispetto ai dati osservati. Se il modello riproduce bene la distribuzione dei dati osservati, si può concludere che è adatto almeno per il campione corrente. In caso contrario, potrebbe essere necessario rivedere le specifiche del modello, come i priors o la struttura delle variabili.

Verifichiamo dunque le predizioni del modello confrontandole con i dati osservati del campione corrente:

pp_check(fit_2)

Nel caso presente, vi è una buona corrispondenza tra i dati simulati dal modello e i dati osservati.

Il grafico seguente analizza gli errori del modello rispetto alla retta di regressione stimata.

pp_check(fit_1, type = "error_scatter_avg")

Questo comando utilizza la funzione pp_check() per produrre un grafico che mostra i residui bayesiani, ovvero le differenze tra i dati osservati e quelli predetti dal modello. Nel tipo specifico di grafico scelto ("error_scatter_avg"), i residui sono rappresentati rispetto ai valori predetti, consentendo di valutare visivamente se sono distribuiti in modo uniforme.

Dal grafico, si osserva che i residui bayesiani appaiono distribuiti in modo omogeneo rispetto alla retta di regressione (che non è direttamente mostrata nel grafico). Questo suggerisce che:

  • Il modello cattura correttamente la relazione tra la variabile predittiva e la variabile di risposta.
  • Non ci sono pattern sistematici nei residui, come deviazioni non lineari o variazioni della dispersione (eteroschedasticità).

Se fossero presenti pattern evidenti nei residui (ad esempio, una struttura curva o una variazione sistematica della dispersione), ciò indicherebbe che il modello potrebbe non essere adeguato, richiedendo una rivalutazione della sua struttura (ad esempio, aggiungendo termini non lineari o trasformando le variabili).

63.8 Regressione Robusta

In questa sezione introduciamo la regressione robusta. Lo scopo è quello di mostrare quanto sia facile modificare il modello definito da brm per specificare una diversa distribuzione degli errori. Questo non è possibile nel caso dell’approccio frequentista.

I modelli robusti sono utili in presenza di outlier. Ad esempio, introduciamo un outlier nei dati:

df_outlier <- df
df_outlier$height[1] <- 200
df_outlier$weight_c[1] <- -15
df_outlier |> 
  ggplot(aes(x = weight_c, y = height)) +
    geom_point() +  
    labs(x = "Weight", y = "Height")

Notiamo come la presenza di un solo outlier introduce una distorsione nei risultati:

fit_3 = brm(
  bf(height ~ 1 + weight_c, center = FALSE), 
  prior = prior_guassian,
  data = df_outlier, 
  backend = "cmdstanr", 
  silent = 0
)
draws <- posterior::as_draws(fit_3, variable = "^b_", regex = TRUE)
posterior::summarise_draws(draws, "mean", "sd", "mcse_mean", "mcse_sd")
#> # A tibble: 2 × 5
#>   variable       mean     sd mcse_mean  mcse_sd
#>   <chr>         <dbl>  <dbl>     <dbl>    <dbl>
#> 1 b_Intercept 155.    0.319   0.00489  0.00481 
#> 2 b_weight_c    0.845 0.0494  0.000792 0.000709

Adattiamo ora un modello robusto utilizzando una distribuzione \(t\) di Student:

fit_4 = brm(
  bf(height ~ 1 + weight_c, center = FALSE), 
  prior = prior_guassian,
  family = student(),
  data = df_outlier, 
  backend = "cmdstanr", 
  silent = 0
)

I risultati mostrano che il modello \(t\) è meno influenzato dagli outlier rispetto al modello gaussiano.

draws <- posterior::as_draws(fit_4, variable = "^b_", regex = TRUE)
posterior::summarise_draws(draws, "mean", "sd", "mcse_mean", "mcse_sd")
#> # A tibble: 2 × 5
#>   variable       mean     sd mcse_mean  mcse_sd
#>   <chr>         <dbl>  <dbl>     <dbl>    <dbl>
#> 1 b_Intercept 155.    0.267   0.00426  0.00429 
#> 2 b_weight_c    0.921 0.0395  0.000643 0.000593

Il parametro \(\nu\) della \(t\) di Student viene stimato dal modello. Nel caso presente

draws <- posterior::as_draws(fit_4, variable = "nu")
posterior::summarise_draws(draws, "mean", "sd", "mcse_mean", "mcse_sd")
#> # A tibble: 1 × 5
#>   variable  mean    sd mcse_mean mcse_sd
#>   <chr>    <dbl> <dbl>     <dbl>   <dbl>
#> 1 nu        6.09  1.69    0.0303  0.0341

Con un parametro \(\nu\) = 6, la \(t\) di Student ha delle “code” molto maggiori di una gaussiana, e questo le consene di “assorbire” gli outliers in maniera maggiore che la gaussiana.

63.9 Indice di Determinazione Bayesiano

Con il pacchetto brms, possiamo calcolare il Bayes \(R^2\), che rappresenta l’equivalente bayesiano del classico indice di determinazione \(R^2\). Questo indice quantifica la proporzione di varianza spiegata dal modello, tenendo conto dell’incertezza intrinseca delle stime bayesiane.

Il comando per calcolarlo è:

bayes_R2(fit_4)
#>    Estimate Est.Error   Q2.5  Q97.5
#> R2   0.5006   0.01975 0.4599 0.5362

Il comando restituisce un tibble (una tabella ordinata) con le seguenti informazioni:

  • Estimate: La stima media del Bayes \(R^2\), cioè la proporzione di varianza spiegata dal modello, basata sulle distribuzioni a posteriori dei parametri.
  • Est.Error: L’errore standard associato alla stima del \(R^2\).
  • Q2.5 e Q97.5: I limiti inferiore e superiore dell’intervallo di credibilità al 95% per il Bayes \(R^2\). Questi valori indicano l’incertezza sul \(R^2\), riflettendo la distribuzione a posteriori.

Nel caso presente

  • Stima del \(R^2\): Il modello spiega in media circa il 50% della varianza osservata nella variabile dipendente.
  • Errore Standard: L’incertezza sulla stima è relativamente bassa (±0.02).
  • Intervallo di Credibilità: C’è un 95% di probabilità che il vero valore del \(R^2\) si trovi tra 0.457 e 0.537.

63.9.1 Differenze rispetto al Frequentista \(R^2\)

  1. Incertezza: Il Bayes \(R^2\) include un’intera distribuzione a posteriori, permettendo di rappresentare l’incertezza attraverso l’intervallo di credibilità. Questo non è possibile con il \(R^2\) frequentista, che fornisce una stima puntuale.
  2. Priors: Il Bayes \(R^2\) è influenzato dai priors scelti per i parametri del modello, il che consente una maggiore flessibilità e incorpora conoscenze preesistenti.

In conclusione, il Bayes \(R^2\) è uno strumento potente per valutare l’adattamento di un modello bayesiano, permettendo di quantificare non solo la proporzione di varianza spiegata, ma anche l’incertezza associata alla stima.

63.10 Approfondimento sulla manipolazione della distribuzione a posteriori con brms

Di seguito illustriamo come accedere e manipolare i campioni generati dal modello bayesiano in brms. Supponiamo di aver costruito un modello lineare semplice, dove vogliamo predire la variabile height in funzione di weight_c:

fit_1 <- brm(
  bf(height ~ 1 + weight_c, center = FALSE), 
  data = df, 
  backend = "cmdstanr", 
  silent = 0
)

Una volta che il modello è stato fit, possiamo estrarre le draws (ossia i campioni) della distribuzione a posteriori tramite la funzione as_draws():

posterior_1 <- as_draws(fit_1)

63.10.1 Struttura dei campioni a posteriori

L’oggetto posterior_1 ottenuto è di tipo draws (una struttura definita dal pacchetto posterior), che internamente può essere rappresentato come una lista o un array in cui sono memorizzati i campioni MCMC:

str(posterior_1)
#> List of 4
#>  $ 1:List of 5
#>   ..$ b_Intercept: num [1:1000] 155 155 154 155 154 ...
#>   ..$ b_weight_c : num [1:1000] 0.886 0.84 0.91 0.888 0.98 ...
#>   ..$ sigma      : num [1:1000] 5.36 5.21 5.07 4.97 5.36 ...
#>   ..$ lprior     : num [1:1000] -2.7 -2.68 -2.67 -2.66 -2.7 ...
#>   ..$ lp__       : num [1:1000] -1077 -1074 -1073 -1072 -1075 ...
#>  $ 2:List of 5
#>   ..$ b_Intercept: num [1:1000] 154 155 155 154 155 ...
#>   ..$ b_weight_c : num [1:1000] 0.967 0.924 0.865 0.929 0.891 ...
#>   ..$ sigma      : num [1:1000] 5.25 5.25 4.88 4.94 4.99 ...
#>   ..$ lprior     : num [1:1000] -2.69 -2.69 -2.66 -2.66 -2.67 ...
#>   ..$ lp__       : num [1:1000] -1074 -1073 -1073 -1073 -1072 ...
#>  $ 3:List of 5
#>   ..$ b_Intercept: num [1:1000] 155 155 154 155 155 ...
#>   ..$ b_weight_c : num [1:1000] 0.953 0.955 0.809 0.933 0.939 ...
#>   ..$ sigma      : num [1:1000] 4.97 4.94 5.12 5.3 5.7 ...
#>   ..$ lprior     : num [1:1000] -2.66 -2.66 -2.68 -2.69 -2.73 ...
#>   ..$ lp__       : num [1:1000] -1073 -1073 -1075 -1075 -1077 ...
#>  $ 4:List of 5
#>   ..$ b_Intercept: num [1:1000] 155 155 154 154 154 ...
#>   ..$ b_weight_c : num [1:1000] 0.905 0.899 0.878 0.913 0.938 ...
#>   ..$ sigma      : num [1:1000] 5.24 4.98 4.79 5.11 5.09 ...
#>   ..$ lprior     : num [1:1000] -2.69 -2.66 -2.65 -2.68 -2.67 ...
#>   ..$ lp__       : num [1:1000] -1072 -1072 -1074 -1072 -1073 ...
#>  - attr(*, "class")= chr [1:3] "draws_list" "draws" "list"

Questa ispezione ti permette di vedere che l’oggetto contiene le catene e i parametri campionati.

63.10.2 Estrazione di un parametro specifico

Possiamo estrarre i nomi dei parametri del modello dall’oggetto creato da brm() nel modo seguente:

variables(fit_1)
#> [1] "b_Intercept" "b_weight_c"  "sigma"       "lprior"      "lp__"

Nel nostro esempio, siamo interessati al coefficiente di regressione associato a weight_c, che in brms è etichettato come b_weight_c. Per semplificare la manipolazione e l’analisi, possiamo servirci di tidybayes, un pacchetto che fornisce funzioni utili per trasformare i campioni in formati “tidy”.

library(tidybayes)

b_slope_draws <- posterior_1 |> 
  spread_draws(b_weight_c)

La funzione spread_draws() estrae e “srotola” le draw dei parametri in un tibble, che risulta più comodo da esplorare:

head(b_slope_draws)
#> # A tibble: 6 × 4
#>   .chain .iteration .draw b_weight_c
#>    <int>      <int> <int>      <dbl>
#> 1      1          1     1      0.886
#> 2      1          2     2      0.840
#> 3      1          3     3      0.910
#> 4      1          4     4      0.888
#> 5      1          5     5      0.980
#> 6      1          6     6      0.884

In questo modo, ogni riga corrisponde a un singolo campione della catena MCMC, per quel parametro.

63.10.3 Calcolo di statistiche riassuntive

Una volta estratti i campioni del parametro di interesse (qui b_weight_c), possiamo calcolare facilmente statistiche come quantili e medie:

quantile(b_slope_draws$b_weight_c, probs = c(0.03, 0.50, 0.97))
#>     3%    50%    97% 
#> 0.8263 0.9060 0.9845
mean(b_slope_draws$b_weight_c)
#> [1] 0.9056
  • I quantili a 0.03, 0.50 e 0.97 forniscono, rispettivamente, un limite inferiore al 94% (0.03–0.97), la mediana a posteriori (0.50) e un limite superiore.
  • La funzione mean() restituisce la media a posteriori del coefficiente, un’altra statistica utile per la stima puntuale.

63.10.4 Visualizzazione della distribuzione a posteriori

Per comprendere meglio la forma della distribuzione a posteriori, è buona prassi tracciarne la densità. Con tidyverse e tidybayes, possiamo creare un grafico elegante in poche righe:

tibble(beta = b_slope_draws$b_weight_c) %>%
  ggplot(aes(x = beta)) +
  stat_halfeye(fill = "skyblue", alpha = 0.6) +
  labs(
    title = "Distribuzione a posteriori di β (b_weight_c)",
    x = "Valore di β",
    y = "Densità a posteriori"
  )

  • stat_halfeye() mostra la densità stimata dei campioni, evidenziando in modo intuitivo i valori più probabili.
  • Il colore e l’alpha permettono di personalizzare l’aspetto del grafico.

In sintesi, utilizzando l’oggetto restituito da brm() e la funzione as_draws(), possiamo estrarre i campioni della distribuzione a posteriori e analizzare:

  • statistiche di sintesi, come media, mediana e quantili;
  • distribuzioni di densità, per una visualizzazione più immediata della variabilità e della forma della distribuzione a posteriori di un parametro.

La combinazione di posterior, tidybayes e tidyverse rende l’intero flusso di lavoro — dall’estrazione dei campioni fino alla creazione di grafici e statistiche — semplice e molto flessibile.

63.11 Riflessioni Conclusive

Questo capitolo ha mostrato come utilizzare brms per costruire e interpretare modelli lineari, evidenziando le sue capacità di gestione dei priors, diagnostica e modellizzazione robusta. Grazie alla sua semplicità e flessibilità, brms rappresenta un potente strumento per l’inferenza bayesiana.

Informazioni sull’Ambiente di Sviluppo

sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.3.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
#> 
#> 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] tidybayes_3.0.7     rstan_2.32.6        StanHeaders_2.32.10
#>  [4] cmdstanr_0.8.1      posterior_1.6.1     brms_2.22.0        
#>  [7] Rcpp_1.0.14         thematic_0.1.6      MetBrewer_0.2.0    
#> [10] ggokabeito_0.1.0    see_0.10.0          gridExtra_2.3      
#> [13] patchwork_1.3.0     bayesplot_1.11.1    psych_2.4.12       
#> [16] scales_1.3.0        markdown_1.13       knitr_1.49         
#> [19] lubridate_1.9.4     forcats_1.0.0       stringr_1.5.1      
#> [22] dplyr_1.1.4         purrr_1.0.4         readr_2.1.5        
#> [25] tidyr_1.3.1         tibble_3.2.1        ggplot2_3.5.1      
#> [28] tidyverse_2.0.0     rio_1.2.3           here_1.0.1         
#> 
#> loaded via a namespace (and not attached):
#>  [1] mnormt_2.1.1         inline_0.3.21        sandwich_3.1-1      
#>  [4] rlang_1.1.5          magrittr_2.0.3       multcomp_1.4-28     
#>  [7] matrixStats_1.5.0    compiler_4.4.2       loo_2.8.0           
#> [10] vctrs_0.6.5          reshape2_1.4.4       arrayhelpers_1.1-0  
#> [13] pkgconfig_2.0.3      fastmap_1.2.0        backports_1.5.0     
#> [16] labeling_0.4.3       utf8_1.2.4           rmarkdown_2.29      
#> [19] tzdb_0.4.0           ps_1.9.0             xfun_0.51           
#> [22] jsonlite_1.9.1       parallel_4.4.2       R6_2.6.1            
#> [25] stringi_1.8.4        estimability_1.5.1   zoo_1.8-13          
#> [28] pacman_0.5.1         R.utils_2.13.0       Matrix_1.7-2        
#> [31] splines_4.4.2        timechange_0.3.0     tidyselect_1.2.1    
#> [34] rstudioapi_0.17.1    abind_1.4-8          yaml_2.3.10         
#> [37] codetools_0.2-20     curl_6.2.1           processx_3.8.6      
#> [40] pkgbuild_1.4.6       lattice_0.22-6       plyr_1.8.9          
#> [43] withr_3.0.2          bridgesampling_1.1-2 coda_0.19-4.1       
#> [46] evaluate_1.0.3       survival_3.8-3       RcppParallel_5.1.10 
#> [49] ggdist_3.3.2         pillar_1.10.1        tensorA_0.36.2.1    
#> [52] checkmate_2.3.2      stats4_4.4.2         distributional_0.5.0
#> [55] generics_0.1.3       rprojroot_2.0.4      hms_1.1.3           
#> [58] rstantools_2.4.0     munsell_0.5.1        xtable_1.8-4        
#> [61] glue_1.8.0           emmeans_1.10.7       tools_4.4.2         
#> [64] data.table_1.17.0    mvtnorm_1.3-3        grid_4.4.2          
#> [67] QuickJSR_1.6.0       colorspace_2.1-1     nlme_3.1-167        
#> [70] cli_3.6.4            svUnit_1.0.6         Brobdingnag_1.2-9   
#> [73] V8_6.0.1             gtable_0.3.6         R.methodsS3_1.8.2   
#> [76] digest_0.6.37        TH.data_1.1-3        htmlwidgets_1.6.4   
#> [79] farver_2.1.2         htmltools_0.5.8.1    R.oo_1.27.0         
#> [82] lifecycle_1.0.4      MASS_7.3-65

Bibliografia

Alter, U., Too, M. A., & Cribbie, R. A. (2025). Navigating the Bayes maze: The psychologist’s guide to Bayesian statistics, a hands-on tutorial with R code. International Journal of Psychology, 60(1), e13271.
McElreath, R. (2020). Statistical rethinking: A Bayesian course with examples in R and Stan (2nd Edition). CRC Press.
Wilkinson, G., & Rogers, C. (1973). Symbolic description of factorial models for analysis of variance. Journal of the Royal Statistical Society Series C: Applied Statistics, 22(3), 392–399.