67  Confronto tra le medie di due gruppi

In questo capitolo imparerai a
  • condurre un confronto bayesiano tra le medie di due gruppi utilizzando la funzione brm() del pacchetto brms.
Prerequisiti
  • Consultare l’articolo “Bayesian estimation supersedes the t test” (Kruschke, 2013).
Preparazione del Notebook 
here::here("code", "_common.R") |> 
  source()

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

67.1 Introduzione

Nella ricerca psicologica, capita spesso di voler confrontare due gruppi per verificare se, ad esempio, uno ha una media più alta dell’altro su una certa variabile. Tuttavia, le differenze osservate nei dati possono essere dovute non solo a reali differenze tra i gruppi, ma anche a fluttuazioni casuali o a rumore di misurazione. Per questo motivo, è fondamentale ricorrere a un modello statistico per valutare se la differenza osservata è effettivamente significativa o semplicemente dovuta al caso.

Tradizionalmente, il confronto tra due gruppi viene effettuato attraverso un test di ipotesi, che si basa sulla formulazione di un’ipotesi nulla (ad esempio: “non ci sono differenze tra i gruppi”) e sull’uso di una statistica test per valutare se i dati osservati sono compatibili con tale ipotesi. Se la statistica supera una certa soglia (detta livello di significatività), l’ipotesi nulla viene rifiutata.

Tuttavia, questo approccio presenta diverse criticità. Come osservato da Johnson (1999) e Goodman (1999), l’uso dei test di ipotesi è spesso guidato da convenzioni arbitrarie, come la scelta del livello di significatività o del tipo di test da applicare. Inoltre, l’output di questi test tende a fornire informazioni indirette, che possono essere facilmente fraintese o sovrainterpretate.

67.1.1 Un’alternativa più informativa: la stima bayesiana

Un approccio più utile e trasparente consiste nel passare dal test alla stima: invece di chiederci “Esiste una differenza?”, ci chiediamo “Quanto è grande la differenza?” e “Quanto siamo incerti su questa stima?”. Questo approccio è reso naturale dal paradigma bayesiano, che permette di esprimere esplicitamente l’incertezza (sia epistemica che aleatoria) e di incorporare conoscenze pregresse tramite distribuzioni a priori.

67.2 Regressione bayesiana per il confronto tra due gruppi

Il confronto tra due medie può essere affrontato con un modello di regressione, che ha il vantaggio di essere facilmente estendibile. Invece di calcolare direttamente la differenza tra medie, introduciamo una variabile indicatrice (dummy) \(D_i\), che assume valore 0 per il gruppo di riferimento e 1 per il gruppo di confronto:

\[ y_i = \alpha + \gamma D_i + \varepsilon_i , \]

dove:

  • \(y_i\) è la variabile osservata per l’unità \(i\),
  • \(\alpha\) rappresenta la media del gruppo con \(D = 0\),
  • \(\gamma\) rappresenta la differenza tra le medie dei due gruppi,
  • \(\varepsilon_i \sim \mathcal{N}(0, \sigma^2)\) rappresenta l’errore casuale.

Espresso in forma bayesiana, il modello diventa:

\[ \begin{aligned} y_i &\sim \mathcal{N}(\mu_i, \sigma) \\ \mu_i &= \alpha + \gamma D_i . \end{aligned} \]

Interpretazione:

  • per il gruppo di riferimento (\(D_i = 0\)): \(\mu_i = \alpha\);
  • per il gruppo di confronto (\(D_i = 1\)): \(\mu_i = \alpha + \gamma\).

Quindi, \(\gamma\) misura la differenza tra le medie. L’inferenza bayesiana si concentra sulla distribuzione a posteriori di \(\gamma\), che rappresenta in modo esplicito ciò che sappiamo (e quanto siamo incerti) sulla differenza tra i due gruppi.

67.3 Approccio Frequentista

Nel paradigma frequentista, l’inferenza sulla differenza tra due gruppi si basa sulla distribuzione campionaria della differenza tra le medie. Supponiamo che i dati provengano da due popolazioni normali:

\[ Y_1 \sim \mathcal{N}(\mu_1, \sigma_1^2) \quad \text{e} \quad Y_2 \sim \mathcal{N}(\mu_2, \sigma_2^2) \]

e che i campioni siano indipendenti. Se assumiamo inoltre che le varianze siano uguali \((\sigma_1^2 = \sigma_2^2 = \sigma^2)\), possiamo usare un modello semplificato.

67.3.1 Statistica di interesse

La quantità che vogliamo stimare è:

\[ \mu_1 - \mu_2 . \]

La corrispondente statistica campionaria è la differenza tra le medie osservate:

\[ \bar{Y}_1 - \bar{Y}_2 . \]

67.3.2 Proprietà della statistica campionaria

Valore atteso della differenza tra le medie campionarie:

\[ E(\bar{Y}_1 - \bar{Y}_2) = \mu_1 - \mu_2 . \]

Supponiamo di avere due campioni indipendenti:

  • \(Y_{11}, Y_{12}, \dots, Y_{1n_1} \sim \text{i.i.d.} \; \mathcal{D}_1\), con \(E(Y_{1i}) = \mu_1\)
  • \(Y_{21}, Y_{22}, \dots, Y_{2n_2} \sim \text{i.i.d.} \; \mathcal{D}_2\), con \(E(Y_{2j}) = \mu_2\)

Le medie campionarie sono:

  • \(\bar{Y}_1 = \dfrac{1}{n_1} \sum_{i=1}^{n_1} Y_{1i}\)
  • \(\bar{Y}_2 = \dfrac{1}{n_2} \sum_{j=1}^{n_2} Y_{2j}\)

Vogliamo calcolare:

\[ E(\bar{Y}_1 - \bar{Y}_2) \]

Per la linearità dell’operatore di valore atteso, abbiamo:

\[ E(\bar{Y}_1 - \bar{Y}_2) = E(\bar{Y}_1) - E(\bar{Y}_2) \]

Calcoliamo i due termini:

1. Valore atteso di \(\bar{Y}_1\)

\[ E(\bar{Y}_1) = E\left( \frac{1}{n_1} \sum_{i=1}^{n_1} Y_{1i} \right) = \frac{1}{n_1} \sum_{i=1}^{n_1} E(Y_{1i}) = \mu_1 \]

2. Valore atteso di \(\bar{Y}_2\)

\[ E(\bar{Y}_2) = E\left( \frac{1}{n_2} \sum_{j=1}^{n_2} Y_{2j} \right) = \frac{1}{n_2} \sum_{j=1}^{n_2} E(Y_{2j}) = \mu_2 \]

Conclusione:

\[ E(\bar{Y}_1 - \bar{Y}_2) = \mu_1 - \mu_2 \]

Varianza della differenza tra le medie campionarie (campioni indipendenti):

\[ \operatorname{Var}(\bar{Y}_1 - \bar{Y}_2) = \frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2} . \]

Abbiamo:

  • \(Y_{1i} \sim \text{i.i.d.}\), con \(E(Y_{1i}) = \mu_1\) e \(\operatorname{Var}(Y_{1i}) = \sigma_1^2\)
  • \(Y_{2j} \sim \text{i.i.d.}\), con \(E(Y_{2j}) = \mu_2\) e \(\operatorname{Var}(Y_{2j}) = \sigma_2^2\)

I due campioni sono indipendenti.

Medie campionarie:

\[ \bar{Y}_1 = \frac{1}{n_1} \sum_{i=1}^{n_1} Y_{1i}, \quad \bar{Y}_2 = \frac{1}{n_2} \sum_{j=1}^{n_2} Y_{2j} \]

Poiché \(\bar{Y}_1\) e \(\bar{Y}_2\) sono indipendenti:

\[ \operatorname{Var}(\bar{Y}_1 - \bar{Y}_2) = \operatorname{Var}(\bar{Y}_1) + \operatorname{Var}(\bar{Y}_2) \]

Per la proprietà della varianza della media campionaria:

\[ \operatorname{Var}(\bar{Y}_1) = \frac{\sigma_1^2}{n_1}, \quad \operatorname{Var}(\bar{Y}_2) = \frac{\sigma_2^2}{n_2} \]

Conclusione:

\[ \operatorname{Var}(\bar{Y}_1 - \bar{Y}_2) = \frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2} \]

Se assumiamo varianze uguali \((\sigma_1 = \sigma_2 = \sigma)\), la formula si semplifica:

\[ \operatorname{Var}(\bar{Y}_1 - \bar{Y}_2) = \sigma^2 \left( \frac{1}{n_1} + \frac{1}{n_2} \right) . \]

Questa quantità può essere stimata con la varianza pooled (varianza combinata dei due campioni):

\[ s_p^2 = \frac{(n_1 - 1)s_1^2 + (n_2 - 1)s_2^2}{n_1 + n_2 - 2} , \]

dove \(s_1^2\) e \(s_2^2\) sono le varianze campionarie calcolate con:

\[ s_j^2 = \frac{1}{n_j - 1} \sum_{i=1}^{n_j} (y_{j,i} - \bar{y}_j)^2, \quad j = 1,2 . \]

67.3.3 Distribuzione della statistica

Se entrambe le popolazioni sono normali e indipendenti, la statistica \(\bar{Y}_1 - \bar{Y}_2\) segue (approssimativamente) una distribuzione normale:

\[ \bar{Y}_1 - \bar{Y}_2 \sim \mathcal{N}\left( \mu_1 - \mu_2,\ \sigma \sqrt{\frac{1}{n_1} + \frac{1}{n_2}} \right). \]

67.4 Sintesi

Abbiamo quindi due approcci distinti:

  • l’approccio frequentista fornisce una stima puntuale e un intervallo di confidenza per la differenza tra le medie;
  • l’approccio bayesiano, invece, fornisce una distribuzione a posteriori della differenza, che quantifica direttamente l’incertezza e permette una interpretazione probabilistica dei risultati.

Nelle sezioni successive vedremo come implementare entrambi gli approcci in R, e come visualizzare le rispettive inferenze.

67.5 Un Esempio Illustrativo

Consideriamo il dataset relativo al quoziente di intelligenza (QI) di un campione di bambini, distinguendo tra quelli le cui madri hanno completato la scuola superiore e quelli le cui madri non l’hanno completata. Vediamo come implementare l’analisi descritta sopra passo per passo.

67.5.1 Esplorazione iniziale dei dati

Carichiamo i dati e osserviamo una sintesi delle prime righe per capire la struttura del dataset:

kidiq <- rio::import(here::here("data", "kidiq.dta"))
kidiq |> 
  head()
#>   kid_score mom_hs mom_iq mom_work mom_age
#> 1        65      1 121.12        4      27
#> 2        98      1  89.36        4      25
#> 3        85      1 115.44        4      27
#> 4        83      1  99.45        3      25
#> 5       115      1  92.75        4      27
#> 6        98      0 107.90        1      18

Successivamente, analizziamo la distribuzione dei bambini nei due gruppi, in base all’educazione delle madri:

kidiq |> 
  group_by(mom_hs) |> 
  summarize(
    avg = mean(kid_score),
    std = sd(kid_score),
    n = n()
)
#> # A tibble: 2 × 4
#>   mom_hs   avg   std     n
#>    <dbl> <dbl> <dbl> <int>
#> 1      0  77.5  22.6    93
#> 2      1  89.3  19.0   341

I risultati mostrano che:

  • 93 bambini hanno madri che non hanno completato la scuola superiore,
  • 341 bambini hanno madri diplomate,

con le medie e deviazioni standard del QI riportate sopra.

La differenza tra le medie del QI dei due gruppi può essere calcolata direttamente come:

mean(kidiq[kidiq$mom_hs == 1, ]$kid_score) - mean(kidiq[kidiq$mom_hs == 0, ]$kid_score)
#> [1] 11.77

Questa analisi preliminare evidenzia la differenza media tra i gruppi.

67.6 Confronto tra Approccio Frequentista e Bayesiano

Il confronto tra le medie di due gruppi indipendenti può essere affrontato sia con un approccio frequentista, basato sulla distribuzione campionaria delle medie, sia con un approccio bayesiano, che stima direttamente la probabilità della differenza tra gruppi. Vediamo i due approcci nel dettaglio, confrontandoli attraverso un esempio pratico.

67.6.1 Approccio Frequentista

Nel paradigma frequentista, l’inferenza si fonda sulla distribuzione campionaria della differenza tra le medie osservate:

\[ \bar{Y}_1 - \bar{Y}_2 \sim \mathcal{N}\left( \mu_1 - \mu_2,\ \sqrt{ \frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2} } \right). \]

Quando le varianze sono considerate uguali (\(\sigma_1 = \sigma_2 = \sigma\)), la varianza comune può essere stimata con la varianza poolata:

\[ s_p^2 = \frac{ (n_1 - 1)s_1^2 + (n_2 - 1)s_2^2 }{n_1 + n_2 - 2}. \]

67.6.1.1 Esempio pratico

Supponiamo di confrontare i punteggi di QI tra figli di madri diplomate e non diplomate:

  • diplomate: \(n_1 = 341\), \(\bar{Y}_1 = 89.3\), \(s_1 = 19.0\),
  • non diplomate: \(n_2 = 93\), \(\bar{Y}_2 = 77.5\), \(s_2 = 22.6\).

Calcoliamo la probabilità che la differenza osservata tra le medie sia superiore a 5 punti:

mean_1 <- 89.3
std_1 <- 19.0
n_1 <- 341

mean_2 <- 77.5
std_2 <- 22.6
n_2 <- 93

pooled_var <- (((n_1 - 1) * std_1^2) + ((n_2 - 1) * std_2^2)) / (n_1 + n_2 - 2)
std_diff <- sqrt(pooled_var / n_1 + pooled_var / n_2)

# Probabilità che la differenza superi 5 punti, sotto H0
pnorm(5, mean = 0, sd = std_diff, lower.tail = FALSE)
#> [1] 0.01553

Questa probabilità (p-value) quantifica quanto è raro osservare una differenza ≥ 5 se le medie fossero uguali.

67.6.2 Approccio Bayesiano

Nel paradigma bayesiano, il confronto è formulato come un modello di regressione con una variabile dummy \(\text{mom\_hs}\), che indica se la madre ha un diploma (1) o no (0):

\[ Y_i \sim \mathcal{N}(\mu_i, \sigma), \quad \mu_i = \beta_0 + \beta_1 \cdot \text{mom\_hs}_i, \]

  • \(\beta_0\): media del gruppo di riferimento (madri non diplomate),
  • \(\beta_1\): differenza tra le medie dei due gruppi,
  • \(\sigma\): deviazione standard dell’errore.

67.6.2.1 Specifica e stima con brms

fit_1 <- brm(
  kid_score ~ mom_hs, 
  data = kidiq, 
  backend = "cmdstanr",
  silent = 0
)

67.6.2.2 Stima della probabilità che la differenza sia maggiore di 5

posterior_samples <- as_draws_df(fit_1)
posterior_samples |> head()
#> # A draws_df: 6 iterations, 1 chains, and 6 variables
#>   b_Intercept b_mom_hs sigma Intercept lprior  lp__
#> 1          79     10.9    20        87   -7.9 -1917
#> 2          79     10.6    20        87   -7.9 -1917
#> 3          78     12.3    21        88   -7.9 -1918
#> 4          78     10.8    20        86   -7.8 -1917
#> 5          74     16.5    20        87   -7.9 -1919
#> 6          81      7.4    20        87   -7.9 -1918
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
mom_hs_coef <- posterior_samples$b_mom_hs
mean(mom_hs_coef > 5)
#> [1] 0.9975

Questa è la proporzione di campioni posteriori in cui \(\beta_1 > 5\): una stima diretta della probabilità che la differenza tra le medie superi 5 punti.

67.6.3 Visualizzazione della distribuzione dei punteggi

Per esplorare le differenze tra i gruppi anche visivamente:

ggplot(kidiq, aes(x = as.factor(mom_hs), y = kid_score)) +
  geom_violin(trim = FALSE, fill = "lightblue") +
  geom_boxplot(width = 0.1, outlier.shape = NA, fill = "white", color = "black") +
  labs(
    x = "Livello di istruzione della madre",
    y = "QI del bambino",
    title = "Distribuzione dei punteggi QI in base all'istruzione materna"
  ) +
  scale_x_discrete(labels = c("0" = "Non diplomata", "1" = "Diplomata"))

Il grafico mostra che i bambini di madri diplomate tendono ad avere punteggi QI più alti, e che le due distribuzioni sono sovrapposte ma asimmetriche.

67.6.4 Differenze chiave tra gli approcci

Aspetto Frequentista Bayesiano
Ipotesi di partenza \(\mu_1 = \mu_2\) (ipotesi nulla) Nessuna ipotesi nulla necessaria
Obiettivo Verificare se una differenza osservata è improbabile sotto \(H_0\) Stimare la distribuzione a posteriori della differenza
Interpretazione P-value: probabilità di osservare un dato risultato sotto \(H_0\) Probabilità che \(\beta_1 > 5\) dati i dati
Incertezza Intervallo di confidenza Intervallo di credibilità
Integrazione conoscenza pregressa Non prevista Ammessa tramite priors

67.6.5 Quale Approccio è Più Utile?

Il metodo frequentista è adatto se si vuole testare la compatibilità dei dati con l’ipotesi \(\mu_1 = \mu_2\). Tuttavia, questa ipotesi non è mai vera nel mondo reale: due gruppi non sono mai esattamente identici.

Il metodo bayesiano, al contrario:

  • non assume che la differenza sia nulla;
  • fornisce una stima diretta della probabilità che la differenza sia maggiore di una soglia;
  • permette di quantificare l’incertezza in modo esplicito;
  • consente di incorporare conoscenze pregresse (o di specificare priors debolmente informativi in loro assenza).

67.6.6 Sintesi

In questo confronto, entrambi gli approcci conducono a risultati coerenti nel senso che indicano una differenza tra i due gruppi. Ma solo il metodo bayesiano consente di rispondere in modo diretto a una domanda intuitiva: quanto è probabile che la differenza superi un certo valore rilevante?

In sintesi, l’approccio bayesiano è particolarmente vantaggioso quando:

  • vogliamo quantificare l’incertezza in modo probabilistico,
  • ci interessa la grandezza della differenza, non solo la sua esistenza,
  • desideriamo una modellazione flessibile e informativa.

67.6.7 Intervallo di credibilità

Sviluppiamo l’analisi bayesiana. Calcoliamo l’intervallo di credibilità all’89% per \(\beta_1\):

bayestestR::hdi(fit_1, parameters = "mom_hs", ci = 0.89)
#> Highest Density Interval
#> 
#> Parameter |       89% HDI
#> -------------------------
#> mom_hs    | [8.36, 15.66]

67.6.8 Distribuzione Predittiva a Posteriori

Effettuiamo un controllo grafico per confrontare i dati osservati con quelli predetti dal modello:

pp_check(fit_1)

67.6.9 R² bayesiano

Infine, calcoliamo il coefficiente di determinazione (Bayesiano \(R^2\)) per valutare la capacità del modello di spiegare la variabilità nei dati:

bayes_R2(fit_1)
#>    Estimate Est.Error    Q2.5  Q97.5
#> R2  0.05766   0.02055 0.02176 0.1007

Questo esempio dimostra come l’analisi di regressione bayesiana consenta di replicare e approfondire i risultati ottenuti in precedenza. Il modello non solo stima la differenza tra i gruppi, ma offre anche un quadro completo dell’incertezza associata, evidenziando la flessibilità e la potenza di questo approccio.

67.7 Una Parametrizzazione Alternativa

Poiché il posterior predictive check (pp-check) ha evidenziato una leggera discrepanza tra i valori osservati (\(y\)) e quelli predetti dal modello, consideriamo una parametrizzazione alternativa. Adottiamo un modello gaussiano esteso con un parametro aggiuntivo per modellare l’asimmetria nella distribuzione, utilizzando una famiglia di distribuzioni skew-normal.

Ecco come definiamo e stimiamo il nuovo modello:

fit_2 <- brm(
  kid_score ~ mom_hs, 
  family = skew_normal(),
  backend = "cmdstanr", 
  silent = 0,
  data = kidiq
)

67.7.1 Verifica del modello

Dopo aver stimato il modello, eseguiamo nuovamente il pp-check per confrontare i dati osservati con quelli predetti:

pp_check(fit_2)

I risultati mostrano un miglioramento nell’adattamento del modello, indicando che l’aggiunta del parametro per l’asimmetria ha contribuito a ridurre le discrepanze.

67.7.2 Valutazione delle stime

Analizziamo le stime posteriori dei parametri per verificare eventuali variazioni rispetto al modello precedente:

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 79.2   1.99    0.0325  0.0290
#> 2 b_mom_hs     9.65  2.24    0.0361  0.0321

67.7.3 Intervallo di credibilità

Calcoliamo l’intervallo di credibilità a densità massima (HDI) per il parametro associato a mom_hs:

bayestestR::hdi(fit_2, parameters = "mom_hs", ci = 0.89)
#> Highest Density Interval
#> 
#> Parameter |       89% HDI
#> -------------------------
#> mom_hs    | [6.10, 13.27]

67.7.4 Valutazione della variabilità spiegata

Infine, calcoliamo il coefficiente di determinazione Bayesiano (\(R^2\)) per quantificare la capacità del modello di spiegare la variabilità nei dati:

bayes_R2(fit_2)
#>    Estimate Est.Error    Q2.5   Q97.5
#> R2  0.04003   0.01711 0.01137 0.07688

In conclusioni, il modello con distribuzione skew-normal offre un adattamento migliore rispetto al modello gaussiano standard, come evidenziato dal pp-check. Tuttavia, le stime posteriori dei parametri differiscono solo marginalmente rispetto al modello precedente. Questo suggerisce che, sebbene l’aggiunta dell’asimmetria migliori l’adattamento, l’effetto sulla stima della relazione tra kid_score e mom_hs è limitato.

67.8 Prior Predictive Checks

Un passaggio fondamentale nella modellazione bayesiana consiste nel valutare le implicazioni predittive dei prior prima di osservare i dati. Questo si realizza attraverso i controlli predittivi a priori (prior predictive checks), ovvero simulazioni che mostrano cosa il modello, sulla base dei soli prior, “si aspetta” di vedere nei dati.

In brms, questi controlli si eseguono in modo simile ai posterior predictive checks, con un’unica differenza: si istruisce brms a campionare solo dai prior, ignorando temporaneamente i dati osservati. Questo richiede che tutti i parametri del modello abbiano prior espliciti — e idealmente, che tali prior siano ragionevoli (cioè non eccessivamente vaghi o implausibili).

67.8.1 Specifica dei prior

Applichiamo questo procedimento a un modello gaussiano sui dati kidiq. Iniziamo ispezionando i prior predefiniti scelti automaticamente da brms:

get_prior(kid_score ~ mom_hs, data = kidiq)
#>                   prior     class   coef group resp dpar nlpar lb ub
#>  student_t(3, 90, 19.3) Intercept                                   
#>                  (flat)         b                                   
#>                  (flat)         b mom_hs                            
#>   student_t(3, 0, 19.3)     sigma                               0   
#>        source
#>       default
#>       default
#>  (vectorized)
#>       default

Quindi li sostituiamo con prior debolmente informativi ma personalizzati:

prior_gaussian <- 
  prior(normal(90, 20), class = "b", coef = "Intercept") +
  prior(normal(0, 15), class = "b", coef = "mom_hs") +
  prior(cauchy(0, 20), class = "sigma")

Questi prior dicono, ad esempio, che:

  • la media del QI (intercetta) è verosimilmente attorno a 90, con ampio margine di incertezza;
  • l’effetto dell’istruzione materna può essere vicino a zero, ma è permessa una variazione ampia;
  • la deviazione standard del QI è positiva e può essere anche piuttosto grande.

67.8.2 Fit del modello con i prior specificati

Stimiamo il modello con i dati, ma usando i nostri prior:

fit_3 <- brm(
  bf(kid_score ~ 1 + mom_hs, center = FALSE), 
  data = kidiq,
  prior = prior_gaussian,
  family = gaussian(),
  backend = "cmdstanr", 
  silent = 0
)

Possiamo poi riassumere le stime a posteriori dei parametri \(\alpha\) e \(\beta\):

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  78.0  2.05    0.0534  0.0414
#> 2 b_mom_hs     11.2  2.29    0.0593  0.0469

67.8.3 Campionamento dalla distribuzione a priori

Per effettuare un prior predictive check, chiediamo a brms di ignorare i dati e di campionare solo dalla distribuzione a priori:

fit_4 <- brm(
  bf(kid_score ~ 1 + mom_hs, center = FALSE), 
  data = kidiq,
  prior = prior_gaussian,
  family = gaussian(),
  backend = "cmdstanr", 
  silent = 0,
  sample_prior = "only"
)

Se esaminiamo il riepilogo di fit_4, noteremo che le “stime” ottenute non derivano dai dati ma riflettono esattamente i prior:

summary(fit_4)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: kid_score ~ 1 + mom_hs 
#>    Data: kidiq (Number of observations: 434) 
#>   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    89.89     20.01    50.65   128.40 1.00     3689     2858
#> mom_hs       -0.04     14.73   -27.63    28.42 1.00     3613     2922
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma   240.42   4297.30     0.66   578.48 1.00     3095     2039
#> 
#> 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).

67.8.4 Visualizzazione del prior predictive check

Utilizziamo pp_check() per confrontare le predizioni generate dai prior con i dati reali:

pp_check(fit_4, ndraws = 100) + xlim(10, 180)

Il grafico mostra la distribuzione dei punteggi simulati (in blu) ottenuti usando solo i prior, senza considerare i dati osservati. L’obiettivo è valutare se i prior producono dati simulati compatibili con la scala reale del fenomeno.

Nel nostro esempio, la distribuzione predittiva a priori è più ampia rispetto ai dati osservati, ma rimane nello stesso ordine di grandezza. Questo è un buon segno: significa che i prior sono abbastanza flessibili da coprire valori plausibili, ma non così ampi da generare predizioni irrealistiche (come punteggi di QI negativi o >200 con alta probabilità).

67.8.5 Perché fare prior predictive checks?

  • Permettono di visualizzare le implicazioni dei prior prima di osservare i dati.
  • Aiutano a diagnosticare prior troppo vaghi, troppo stretti o implausibili.
  • Sono particolarmente utili quando si lavora con dati scarsi, perché in quel caso i prior possono influenzare fortemente le stime.

⚠️ Un prior che sembra “ragionevole” sulla carta può generare predizioni completamente assurde quando combinato con la struttura del modello. Il controllo predittivo a priori è quindi un test visivo fondamentale per valutare la coerenza del modello.

67.8.6 Esempio con Prior Informativi

Un ulteriore approfondimento dell’approccio bayesiano consiste nell’uso di prior informativi. Ad esempio, se studi precedenti indicano che le differenze tra i due gruppi tendono ad essere intorno a 10 punti, possiamo specificare un prior normale centrato su 10 con una deviazione standard moderata. Questo prior può influenzare i risultati quando i dati sono scarsi o poco informativi, ma con un campione ampio come in questo caso, i dati tendono a prevalere sui prior.

fit_5 <- brm(
  kid_score ~ mom_hs,
  data = kidiq,
  prior = c(set_prior("normal(10, 5)", class = "b", coef = "mom_hs")),
  backend = "cmdstanr",
  silent = 0
)
summary(fit_5)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: kid_score ~ mom_hs 
#>    Data: kidiq (Number of observations: 434) 
#>   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    77.81      1.92    74.07    81.58 1.00     4427     3120
#> mom_hs       11.43      2.12     7.37    15.60 1.00     4265     3078
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma    19.88      0.66    18.62    21.20 1.00     4270     3219
#> 
#> 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).

Con prior informativi, i risultati rifletteranno sia le evidenze dei dati sia il nostro stato di conoscenza precedente, arricchendo ulteriormente l’inferenza.

67.8.7 Test di Ipotesi Bayesiano con hypothesis()

Il comando hypothesis() in brms consente di verificare ipotesi specifiche direttamente sulla distribuzione a posteriori. Ad esempio:

hypothesis(fit_1, "mom_hs > 5")
#> Hypothesis Tests for class b:
#>         Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
#> 1 (mom_hs)-(5) > 0     6.78      2.31     3.02    10.54        399
#>   Post.Prob Star
#> 1         1    *
#> ---
#> 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
#> '*': For one-sided hypotheses, the posterior probability exceeds 95%;
#> for two-sided hypotheses, the value tested against lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.

Confrontando questo approccio con il test frequentista, emergono differenze sostanziali nel modo in cui vengono comunicate e interpretate le evidenze:

  • Il test bayesiano non restituisce un p-value, ma una probabilità interpretabile in termini diretti. Ad esempio, l’output di hypothesis(fit_1, "mom_hs > 5") ci dice che la probabilità a posteriori che la differenza tra i gruppi sia maggiore di 5 punti è circa 100%. Questa interpretazione è intuitiva e facilmente comunicabile.

  • L’output include anche un intervallo di credibilità (ad esempio, 90%), che rappresenta l’intervallo entro cui si trova il parametro con una certa probabilità a posteriori. A differenza dell’intervallo di confidenza frequentista, qui possiamo affermare, ad esempio, che c’è una probabilità del 90% che la vera differenza tra gruppi sia compresa tra 3.02 e 10.54 punti.

In sintesi, mentre il test frequentista quantifica quanto sarebbe raro osservare un certo risultato se l’ipotesi nulla fosse vera, il test bayesiano stima direttamente quanto è plausibile che un certo effetto superi una soglia di interesse, alla luce dei dati osservati e dei prior specificati.

67.9 Riflessioni Conclusive

In questo capitolo abbiamo introdotto l’uso del pacchetto brms come strumento per la modellazione bayesiana, illustrandone la sintassi intuitiva e la grande flessibilità. Abbiamo visto come, a partire dalla specifica del modello fino all’interpretazione dei risultati, sia possibile tradurre domande di ricerca complesse in analisi statistiche rigorose e trasparenti.

Uno degli aspetti centrali emersi riguarda i vantaggi concettuali dell’approccio bayesiano rispetto a quello frequentista. In particolare, l’inferenza bayesiana:

  • consente di quantificare direttamente l’incertezza sui parametri,
  • permette di calcolare probabilità su quantità di interesse, come la probabilità che una differenza superi una soglia specifica,
  • e offre la possibilità di incorporare conoscenze pregresse in modo formale attraverso i prior.

Questo approccio si dimostra particolarmente utile in contesti in cui l’ipotesi nulla (ad esempio, che due gruppi abbiano medie esattamente uguali) è poco realistica o poco informativa. A differenza dei test tradizionali, come il t-test di Student, che si fondano su una logica dicotomica e su soglie arbitrarie, l’approccio bayesiano consente di esplorare la gradualità dell’evidenza, concentrandosi su ciò che i dati effettivamente suggeriscono.

Abbiamo inoltre visto che:

  • la regressione bayesiana offre un modo naturale per stimare la differenza tra gruppi, e lo fa producendo una distribuzione a posteriori del parametro di interesse, da cui è possibile derivare probabilità e intervalli di credibilità interpretabili in modo diretto;
  • i prior predictive checks aiutano a valutare se le ipotesi iniziali sui parametri (i prior) generano predizioni compatibili con i dati attesi, migliorando la robustezza del modello;
  • la distribuzione a posteriori consente di esplorare non solo se una differenza esiste, ma quanto è plausibile e con quale grado di incertezza.

In sintesi, l’approccio bayesiano si rivela particolarmente efficace quando si desidera:

  1. evitare assunzioni irrealistiche legate all’ipotesi nulla,
  2. formulare inferenze probabilistiche interpretabili in modo intuitivo,
  3. adattare il modello al contesto della ricerca attraverso l’uso di informazioni a priori,
  4. affrontare modelli complessi in modo trasparente, propagando correttamente l’incertezza.

Con l’aumentare della complessità dei dati e delle domande di ricerca, il pensiero bayesiano offre una struttura concettuale più coerente e flessibile, capace di rispondere a interrogativi sfumati e contestualizzati, invece di limitarsi a un “sì o no” basato su un valore soglia.

In conclusione, brms non è solo un pacchetto per eseguire analisi bayesiane, ma uno strumento che favorisce un cambiamento nella prospettiva analitica: incoraggia a pensare in termini di distribuzioni, incertezza e credenze aggiornabili. Questo capitolo costituisce una base per applicare i modelli bayesiani a scenari reali, offrendo non solo strumenti tecnici, ma anche una mentalità più critica, riflessiva e coerente con le sfide della ricerca quantitativa contemporanea.

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] insight_1.2.0     bayestestR_0.15.3 brms_2.22.0       Rcpp_1.0.14      
#>  [5] posterior_1.6.1   cmdstanr_0.9.0    thematic_0.1.6    MetBrewer_0.2.0  
#>  [9] ggokabeito_0.1.0  see_0.11.0        gridExtra_2.3     patchwork_1.3.0  
#> [13] bayesplot_1.12.0  psych_2.5.3       scales_1.4.0      markdown_2.0     
#> [17] knitr_1.50        lubridate_1.9.4   forcats_1.0.0     stringr_1.5.1    
#> [21] dplyr_1.1.4       purrr_1.0.4       readr_2.1.5       tidyr_1.3.1      
#> [25] tibble_3.2.1      ggplot2_3.5.2     tidyverse_2.0.0   rio_1.2.3        
#> [29] 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] R.utils_2.13.0       fastmap_1.2.0        tensorA_0.36.2.1    
#>  [7] pacman_0.5.1         digest_0.6.37        timechange_0.3.0    
#> [10] estimability_1.5.1   lifecycle_1.0.4      StanHeaders_2.32.10 
#> [13] processx_3.8.6       magrittr_2.0.3       compiler_4.5.0      
#> [16] rlang_1.1.6          tools_4.5.0          utf8_1.2.5          
#> [19] yaml_2.3.10          data.table_1.17.0    labeling_0.4.3      
#> [22] bridgesampling_1.1-2 htmlwidgets_1.6.4    curl_6.2.2          
#> [25] pkgbuild_1.4.7       mnormt_2.1.1         plyr_1.8.9          
#> [28] RColorBrewer_1.1-3   abind_1.4-8          withr_3.0.2         
#> [31] R.oo_1.27.1          datawizard_1.0.2     stats4_4.5.0        
#> [34] grid_4.5.0           colorspace_2.1-1     inline_0.3.21       
#> [37] xtable_1.8-4         emmeans_1.11.1       cli_3.6.5           
#> [40] mvtnorm_1.3-3        rmarkdown_2.29       generics_0.1.3      
#> [43] RcppParallel_5.1.10  rstudioapi_0.17.1    reshape2_1.4.4      
#> [46] tzdb_0.5.0           rstan_2.32.7         parallel_4.5.0      
#> [49] matrixStats_1.5.0    vctrs_0.6.5          V8_6.0.3            
#> [52] Matrix_1.7-3         jsonlite_2.0.0       hms_1.1.3           
#> [55] glue_1.8.0           codetools_0.2-20     ps_1.9.1            
#> [58] distributional_0.5.0 stringi_1.8.7        gtable_0.3.6        
#> [61] QuickJSR_1.7.0       pillar_1.10.2        htmltools_0.5.8.1   
#> [64] Brobdingnag_1.2-9    R6_2.6.1             rprojroot_2.0.4     
#> [67] evaluate_1.0.3       lattice_0.22-7       haven_2.5.4         
#> [70] R.methodsS3_1.8.2    backports_1.5.0      rstantools_2.4.0    
#> [73] coda_0.19-4.1        nlme_3.1-168         checkmate_2.3.2     
#> [76] xfun_0.52            pkgconfig_2.0.3

Bibliografia

Goodman, S. N. (1999). Toward evidence-based medical statistics. 1: The P value fallacy. Annals of internal medicine, 130(12), 995–1004.
Johnson, D. H. (1999). The insignificance of statistical significance testing. The journal of wildlife management, 763–772.
Kruschke, J. K. (2013). Bayesian estimation supersedes the t test. Journal of Experimental Psychology: General, 142(2), 573–603.