56 Distribuzione predittiva a posteriori

In questo capitolo imparerai a:

utilizzare la distribuzione predittiva a posteriori per fare previsioni sui dati futuri, incorporando sia l’incertezza sui parametri del modello sia la variabilità intrinseca del processo generativo.
comprendere il ruolo della distribuzione predittiva a posteriori nel verificare la coerenza del modello bayesiano rispetto ai dati osservati.
applicare il concetto di distribuzione predittiva a posteriori al caso beta-binomiale per un’analisi pratica e intuitiva.

Prerequisiti

Leggere il capitolo Posterior Inference & Prediction di Bayes Rules!.
Consultare Bayesian statistics and modelling (Schoot et al., 2021).

Preparazione del Notebook

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

56.1 Introduzione

Nel contesto dell’inferenza bayesiana, uno degli obiettivi principali è non solo stimare i parametri di un modello (ad esempio, la probabilità \(p\) di successo in un esperimento binomiale) ma anche fare previsioni su dati futuri basandosi su ciò che abbiamo osservato. La distribuzione predittiva a posteriori risponde proprio a questa esigenza, combinando:

La nostra incertezza sui parametri descritta dalla distribuzione a posteriori.
La variabilità intrinseca del processo che genera i dati futuri.

In termini semplici, la distribuzione predittiva a posteriori ci dice quali risultati futuri sono plausibili, dato ciò che sappiamo dai dati osservati e dal modello.

56.2 Definizione Formale

Supponiamo di avere un insieme di dati osservati \(y = \{y_1, y_2, \ldots, y_n\}\), generati da un modello che dipende da un parametro sconosciuto \(\theta\). Il parametro \(\theta\) può rappresentare qualsiasi caratteristica del modello, come una probabilità di successo, una media, o un coefficiente in un modello di regressione. Inizialmente, la nostra conoscenza su \(\theta\) è rappresentata dalla distribuzione a priori \(p(\theta)\), che riflette ciò che sappiamo (o non sappiamo) su \(\theta\) prima di osservare i dati.

Dopo aver osservato i dati \(y\), possiamo aggiornare la nostra conoscenza su \(\theta\) utilizzando la formula di Bayes per calcolare la distribuzione a posteriori \(p(\theta \mid y)\):

\[ p(\theta \mid y) = \frac{p(y \mid \theta) p(\theta)}{p(y)}, \]

dove:

\(p(\theta \mid y)\): la distribuzione a posteriori rappresenta la nostra conoscenza aggiornata su \(\theta\) dopo aver osservato i dati.
\(p(y \mid \theta)\): la verosimiglianza è la probabilità di osservare i dati dati i parametri del modello.
\(p(\theta)\): la distribuzione a priori rappresenta la conoscenza iniziale su \(\theta\).
\(p(y)\): l’evidenza è la probabilità totale dei dati osservati, calcolata come:

\[ p(y) = \int p(y \mid \theta) p(\theta) \, d\theta. \]

56.2.1 Previsione di Nuovi Dati

Quando vogliamo prevedere un nuovo dato, indicato con \(\tilde{y}\), la nostra attenzione si sposta sulla distribuzione predittiva a posteriori, \(p(\tilde{y} \mid y)\).

56.2.1.1 Cosa rappresenta \(\tilde{y}\)?

\(\tilde{y}\) rappresenta un dato futuro o non osservato. Ad esempio, se i dati \(y\) rappresentano il numero di successi osservati in una serie di lanci di una moneta, \(\tilde{y}\) potrebbe rappresentare il numero di successi in una nuova serie di lanci.
L’obiettivo è stimare \(p(\tilde{y} \mid y)\), cioè la probabilità del nuovo dato \(\tilde{y}\) dato ciò che abbiamo osservato in \(y\).

56.2.1.2 Cosa rappresenta \(p(\tilde{y} \mid \theta)\)?

\(p(\tilde{y} \mid \theta)\) è la probabilità del nuovo dato \(\tilde{y}\) dato un particolare valore del parametro \(\theta\).
Ad esempio, in un modello binomiale, \(p(\tilde{y} \mid \theta)\) corrisponde alla probabilità di ottenere \(\tilde{y}\) successi su \(n_{\text{new}}\) prove, data la probabilità di successo \(\theta\).

56.2.1.3 Combinazione di \(p(\tilde{y} \mid \theta)\) con \(p(\theta \mid y)\)

Poiché non conosciamo esattamente \(\theta\), dobbiamo considerare tutte le possibili ipotesi su \(\theta\), pesandole in base alla loro probabilità a posteriori \(p(\theta \mid y)\). Questo porta alla formula per la distribuzione predittiva a posteriori:

\[ p(\tilde{y} \mid y) = \int p(\tilde{y} \mid \theta) p(\theta \mid y) \, d\theta. \]

56.2.1.4 Interpretazione

La distribuzione predittiva a posteriori, \(p(\tilde{y} \mid y)\), rappresenta la nostra miglior stima della probabilità del nuovo dato \(\tilde{y}\), tenendo conto sia dei dati osservati \(y\) sia dell’incertezza su \(\theta\).

56.2.2 Caso Discreto

Se il parametro \(\theta\) assume un numero finito di valori, l’integrale si semplifica in una somma:

\[ p(\tilde{y} \mid y) = \sum_{\theta} p(\tilde{y} \mid \theta) p(\theta \mid y). \]

Questo approccio è utile nei modelli discreti o quando si approssimano i parametri con un numero finito di valori campionati.

56.3 Il Caso Beta-Binomiale

Consideriamo un classico esperimento binomiale: lanciare una moneta \(n\) volte e osservare il numero di successi \(y\) (ad esempio, il numero di “teste”). In un contesto bayesiano, il processo di analisi si articola in tre fasi:

Distribuzione a Priori: Prima di osservare i dati, formuliamo una distribuzione a priori sulla probabilità \(p\) di successo, che riflette le nostre conoscenze iniziali (o la loro assenza). Una scelta comune è la distribuzione Beta(\(\alpha, \beta\)), perché è flessibile e ben definita per probabilità. Per esempio:
- \(\alpha\) rappresenta il numero “fittizio” di successi osservati.
- \(\beta\) rappresenta il numero “fittizio” di insuccessi osservati.
Distribuzione a Posteriori: Dopo aver osservato \(y\) successi su \(n\) prove, aggiorniamo la distribuzione a priori combinandola con i dati osservati, ottenendo la distribuzione a posteriori di \(p\):

\[ \alpha_{\text{post}} = \alpha_{\text{prior}} + y, \quad \beta_{\text{post}} = \beta_{\text{prior}} + (n - y). \]

Questa distribuzione rappresenta ciò che sappiamo di \(p\) dopo aver osservato i dati.
Distribuzione Predittiva a Posteriori: Per prevedere il numero di successi futuri \(y_{\text{new}}\) in un nuovo esperimento con \(n_{\text{new}}\) prove, combiniamo la distribuzione a posteriori di \(p\) con la variabilità intrinseca del processo binomiale. In pratica:
- Campioniamo \(p\) dalla distribuzione a posteriori (\(p \sim \text{Beta}(\alpha_{\text{post}}, \beta_{\text{post}})\)).
- Usiamo ciascun campione di \(p\) per simulare nuovi dati (\(y_{\text{new}} \sim \text{Binomiale}(n_{\text{new}}, p)\)).

Questo processo tiene conto sia dell’incertezza sui parametri (\(p\)) sia della variabilità intrinseca nei dati futuri.

56.4 Simulazione della Distribuzione Predittiva a Posteriori

56.4.1 Impostazione dei Parametri

Supponiamo di aver osservato \(y = 70\) successi su \(n = 100\) prove.
Utilizziamo una distribuzione a priori Beta(\(2, 2\)), che rappresenta una conoscenza iniziale debolmente informativa, con una leggera preferenza per \(p \approx 0.5\).

56.4.2 Calcolo della Distribuzione a Posteriori

Aggiorniamo la distribuzione a priori con i dati osservati. I parametri aggiornati sono:

\[ \alpha_{\text{post}} = \alpha_{\text{prior}} + y = 2 + 70 = 72, \quad \beta_{\text{post}} = \beta_{\text{prior}} + (n - y) = 2 + 30 = 32. \]
La distribuzione a posteriori è quindi \(p \sim \text{Beta}(72, 32)\), che descrive la probabilità aggiornata di successo basata sui dati.

56.4.3 Simulazione dei Dati Futuri

Generiamo \(n_{\text{sim}} = 1000\) campioni di \(p\) dalla distribuzione a posteriori Beta(72, 32).
Per ogni campione di \(p\), simuliamo \(y_{\text{new}}\), il numero di successi futuri su \(n_{\text{new}} = 10\) prove, utilizzando la distribuzione binomiale:

\[ y_{\text{new}} \sim \text{Binomiale}(n_{\text{new}} = 10, p). \]
Infine, convertiamo \(y_{\text{new}}\) in proporzioni predette:

\[ \text{Proporzione predetta} = \frac{y_{\text{new}}}{n_{\text{new}}}. \]

# Impostazione del seed per riproducibilità
set.seed(123)

# Parametri osservati
y <- 70
n <- 100

# Parametri a priori
alpha_prior <- 2
beta_prior <- 2

# Calcolo dei parametri a posteriori
alpha_post <- alpha_prior + y
beta_post <- beta_prior + (n - y)

# Simulazione della distribuzione a posteriori
p_samples <- rbeta(1000, alpha_post, beta_post)

# Simulazione di nuovi dati per n_new = 10
y_preds <- sapply(p_samples, function(p) {
  rbinom(1, 10, p)
})

# Calcolo delle proporzioni predette
prop_preds <- y_preds / 10

56.4.4 Risultati Attesi

Distribuzione a Posteriori di \(p\):
- Centrata attorno a \(0.7\), con una varianza ridotta grazie alla dimensione del campione \(n = 100\).
Distribuzione Predittiva per \(n_{\text{new}} = 10\):
- Variabilità più ampia rispetto a \(p\), dovuta al numero ridotto di prove (\(n_{\text{new}} = 10\)).
- Riflette sia l’incertezza su \(p\) sia la variabilità intrinseca dei dati futuri.

# Visualizzazione
par(mfrow = c(1, 2))
curve(dbeta(x, alpha_prior, beta_prior), from = 0, to = 1, 
      main = "Distribuzione a Priori", col = "blue", lwd = 2)
curve(dbeta(x, alpha_post, beta_post), from = 0, to = 1, 
      main = "Distribuzione a Posteriori", col = "red", lwd = 2)

hist(prop_preds, breaks = 20, col = "lightblue", freq = FALSE,
     main = "Distribuzione Predittiva (n_new = 10)",
     xlab = "Proporzione di Successi")
abline(v = y / n, col = "blue", lwd = 3, lty = 2)

56.4.5 Spiegazione del Codice

Simulazione di nuovi dati per \(n_{\text{new}} = 10\):

y_preds <- sapply(p_samples, function(p) {
  rbinom(1, 10, p)
})

p_samples contiene \(1000\) valori di \(p\) simulati dalla distribuzione Beta(72, 32).
Per ciascun \(p\), rbinom(1, 10, p) genera un valore di \(y_{\text{new}}\), simulando il numero di successi futuri su \(n_{\text{new}} = 10\) prove.

Risultato: un vettore di 1000 valori di \(y_{\text{new}}\), uno per ciascun campione di \(p\).

Calcolo delle proporzioni predette:

prop_preds <- y_preds / 10

Divide ciascun valore di \(y_{\text{new}}\) per \(n_{\text{new}} = 10\), ottenendo le proporzioni di successi predette.

Risultato: un vettore di 1000 proporzioni predette (\(y_{\text{new}} / n_{\text{new}}\)).

56.4.6 Interpretazione

La distribuzione a posteriori di \(p\), Beta(72, 32), è centrata attorno a \(p \approx 0.7\), coerente con i dati osservati (\(y / n = 0.7\)).
Le proporzioni predette mostrano la variabilità combinata dell’incertezza su \(p\) (dalla distribuzione a posteriori) e della variabilità binomiale per un campione futuro di 10 prove.
L’istogramma delle proporzioni predette è più ampio rispetto alla distribuzione a posteriori di \(p\), riflettendo l’incertezza aggiuntiva derivante dal numero ridotto di prove (\(n_{\text{new}} = 10\)).

56.5 Posterior Predictive Check (PP-Check)

La distribuzione predittiva a posteriori ottenuta dalla simulazione (\(n_{\text{new}} = 10\)) può essere utilizzata per effettuare un Posterior Predictive Check (PP-Check). Questo controllo confronta i dati osservati con i dati simulati dal modello per verificare se il modello è in grado di riprodurre caratteristiche rilevanti dei dati osservati.

In questa simulazione, il Posterior Predictive-Check suggerisce che il modello è ben specificato per i dati osservati (\(y = 70\), \(n = 100\)): la proporzione osservata (\(0.7\)) è vicina al centro della distribuzione predittiva. Questo significa che il modello può essere considerato valido per fare previsioni sui dati futuri, almeno nel contesto specificato.

Il fatto che la distribuzione predittiva a posteriori rappresenti accuratamente la proporzione osservata può sembrare ovvio nel caso presente. Questo avviene perché stiamo lavorando con un modello semplice e ben specificato, utilizzato qui con l’intento di chiarire la logica del concetto di distribuzione predittiva a posteriori. Tuttavia, nei modelli più complessi, tipici delle indagini psicologiche, la corrispondenza tra i dati osservati e quelli predetti dal modello non può mai essere data per scontata. È essenziale verificarla attraverso la distribuzione predittiva a posteriori.

Se i dati osservati non sono ben rappresentati dalla distribuzione predittiva a posteriori, ciò indica che il modello non è adeguato a spiegare i dati e necessita di revisione. Questo processo di verifica non solo garantisce la coerenza del modello con i dati disponibili, ma consente anche di identificare eventuali aree problematiche nella specificazione del modello, come prior inappropriati o assunzioni non realistiche. In definitiva, il confronto tra i dati osservati e quelli predetti è un passo fondamentale per la validazione di modelli complessi in contesti psicologici e scientifici.

56.6 Riflessioni Conclusive

La distribuzione predittiva a posteriori è uno strumento centrale nell’inferenza bayesiana, poiché consente di fare previsioni sui dati futuri integrando l’incertezza sui parametri del modello con la variabilità intrinseca del processo generativo. Questa capacità va oltre la semplice stima dei parametri, permettendo di confrontare le previsioni del modello con i dati reali, un passaggio fondamentale per verificare la coerenza e l’utilità del modello stesso.

Nel flusso di lavoro bayesiano, la distribuzione predittiva a posteriori svolge un ruolo chiave nella valutazione del modello. Ad esempio, consente di effettuare controlli predittivi a posteriori per identificare discrepanze tra i dati osservati e quelli previsti. Tali controlli aiutano a diagnosticare problemi di specificazione del modello, a valutare l’adeguatezza delle scelte a priori e a guidare eventuali revisioni del modello.

Inoltre, il caso beta-binomiale utilizzato in questo capitolo rappresenta un esempio intuitivo e potente: evidenzia come l’incertezza sui parametri possa essere tradotta in previsioni probabilistiche robuste, senza la necessità di fare assunzioni rigide o non realistiche. Questo approccio non solo formalizza l’incertezza in modo rigoroso, ma permette anche di comunicare le previsioni in modo trasparente e interpretabile, caratteristiche essenziali in ambito decisionale e scientifico.

In sintesi, la distribuzione predittiva a posteriori è un elemento fondamentale della modellazione bayesiana, che lega l’inferenza paramatrica alla previsione empirica, contribuendo a rendere l’intero processo inferenziale più affidabile, interpretabile e applicabile a scenari complessi.

Esercizi

Problemi

Consideriamo i dati della SWLS somministrata a un campione di studenti, ottenendo per ciascuno uno score complessivo. Per semplicità, vogliamo “dichiarare positivo” lo studente se il punteggio SWLS supera una determinata soglia (ad esempio, 20 su 35). In questo modo otteniamo una variabile dicotomica (0/1), che useremo come “successo” in un modello binomiale.

Dati e conteggio dei successi
- Carica il dataset con le risposte SWLS.
- Costruisci la variabile binaria (ad esempio SWLS_dich) che vale 1 se lo score ≥ 20, e 0 altrimenti.
- Calcola il numero di successi (numero di persone che superano la soglia) e il numero totale di osservazioni (N).
Modello beta-binomiale (approccio manuale via simulazione)
- Specifica una distribuzione Beta(a, b) come prior per la probabilità di successo \(p\). Scegli una coppia \((a, b)\) relativamente poco informativa, ad esempio (2,2) o (1,1).
- Osservando \(y\) successi su \(n\) soggetti, aggiorna i parametri a posteriori: \[ a_{\text{post}} = a + y, \quad b_{\text{post}} = b + (n - y). \]
- Simula un gran numero di campioni di \(p\) dalla distribuzione Beta\(\bigl(a_{\text{post}},\, b_{\text{post}}\bigr)\).
- Per ciascun campione di \(p\), genera un valore \(\tilde{y}\) da una Binomiale\(\bigl(n_{\text{new}}, p\bigr)\), dove \(n_{\text{new}}\) è la dimensione di un ipotetico nuovo campione (che puoi scegliere, ad esempio, uguale a \(n\) oppure un valore diverso). Otterrai così una posterior predictive distribution per \(\tilde{y}\).
- Infine, calcola statistiche descrittive (media, varianza, intervalli) e/o disegna un istogramma di \(\tilde{y}\) o della proporzione \(\tilde{y}/n_{\text{new}}\).

Replicare con brms

Usa il pacchetto brms per costruire un modello binomiale. Per esempio:

library(brms)

# Crea un data frame con la variabile dicotomica
df_binom <- data.frame(
  successes = y,    # conteggio dei successi
  failures  = n - y
)

# Modello binomiale con prior Beta(a,b) approssimato tramite logit
fit_brms <- brm(
  bf(successes | trials(n) ~ 1), 
  data = df_binom,
  family = binomial(link = "logit"),
  prior = c(
    prior(beta(2, 2), class = "Intercept", dpar = "mu") 
    # NOTA: la specifica di una "beta(2,2)" diretta sull'intercetta
    # è un'approssimazione, tipicamente serve passare a una scala logit.
    # In brms, di solito si usa prior su scale normali dell'intercetta.
  ),
  seed = 123
)

(Le specifiche del prior potrebbero richiedere una formulazione differente se vuoi rispettare esattamente la corrispondenza con Beta(a,b). In ogni caso, l’idea è mostrare come definire un prior e costruire un modello binomiale con brms.)

Verifica la convergenza e poi estrai la posterior predictive distribution con le funzioni di brms:
```
pp_check(fit_brms, nsamples = 100)
```
Questo ti mostrerà come i dati predetti dal modello (in termini di binomiale) si confrontano con i dati osservati.

Confronto e interpretazione
- Metti a confronto i risultati della simulazione “manuale” (Beta-Binomial) e quelli ottenuti con brms. Noterai che le distribuzioni predittive dovrebbero essere coerenti, se hai impostato un prior per brms simile a quello del modello Beta-Binomiale.
- Discuti brevemente se la distribuzione predittiva a posteriori acquisita è plausibile rispetto ai dati osservati. Ad esempio, la probabilità di osservare \(\tilde{y}\) simile a \(y\) dovrebbe essere relativamente alta se il modello è appropriato.
- Se vuoi, puoi cambiare \(n_{\text{new}}\) (es. previsione su 200 soggetti futuri) per vedere come la variabilità della previsione si “ridimensiona” o cresce a seconda della taglia del campione.

Soluzioni

Costruzione del dataset

Se la SWLS varia tra 5 e 35, e la soglia è 20, puoi fare:

df$SWLS_dich <- ifelse(df$SWLS_score >= 20, 1, 0)
y <- sum(df$SWLS_dich)
n <- nrow(df)

Approccio Beta-Binomial manuale

Prior: \((a, b) = (2, 2)\)
Posterior: \((a_{\text{post}}, b_{\text{post}}) = (2 + y,\, 2 + n - y)\).

Generazione dei campioni:

N_sim <- 2000
p_post <- rbeta(N_sim, a_post, b_post)
y_pred <- rbinom(N_sim, size = n_new, prob = p_post)

# Se preferisci la proporzione futura:
prop_pred <- y_pred / n_new

Statistiche:

mean_p <- mean(p_post) # media a posteriori di p
quantile_p <- quantile(p_post, c(0.025, 0.975))  

mean_prop_pred <- mean(prop_pred)
quantile_prop_pred <- quantile(prop_pred, c(0.025, 0.975))

Grafici (istogramma e densità):

hist(prop_pred, freq=FALSE, col='lightblue',
     main='Posterior Predictive Distribution: prop. di successi')

Modello con brms

Usa la sintassi di una binomiale con offset o con trials(n).

Specifica un prior che approssimi Beta(2,2) sullo scale logit, ad esempio:

# Beta(2,2) ha media ~ 0.5, varianza relativamente ampia.
# Approssimandola su scala logit ~ normal(0, 2.2) 
# (valore indicativo: la normal(0, 2) su logit copre un intervallo ampio).

prior_approx <- prior(normal(0, 2), class = "Intercept")

Esegui pp_check(fit_brms) e interpreta.

Interpretazione
- Se la soglia scelta per la SWLS cattura un “buon livello di soddisfazione”, potresti aspettarti una certa % di successi.
- Se i dati futuri simulati sono coerenti con i dati reali — ad esempio, la media di \(\tilde{y}\) è vicina a \(y\) — allora il modello sembra descrivere bene la realtà. Altrimenti, potresti rivedere la soglia o la specifica del prior.

L’elemento chiave è che la distribuzione predittiva a posteriori (posterior predictive distribution) non si limita a considerare un solo valore di \(p\), bensì campiona molteplici valori plausibili (dalla posterior), e per ciascuno simula un potenziale outcome. Così facendo, si riflette pienamente l’incertezza residua sul parametro e l’aleatorietà del processo binomiale.

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] thematic_0.1.6   MetBrewer_0.2.0  ggokabeito_0.1.0 see_0.10.0      
#>  [5] gridExtra_2.3    patchwork_1.3.0  bayesplot_1.11.1 psych_2.4.12    
#>  [9] scales_1.3.0     markdown_1.13    knitr_1.49       lubridate_1.9.4 
#> [13] forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4      purrr_1.0.4     
#> [17] readr_2.1.5      tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.1   
#> [21] tidyverse_2.0.0  rio_1.2.3        here_1.0.1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] generics_0.1.3    stringi_1.8.4     lattice_0.22-6    hms_1.1.3        
#>  [5] digest_0.6.37     magrittr_2.0.3    evaluate_1.0.3    grid_4.4.2       
#>  [9] timechange_0.3.0  fastmap_1.2.0     rprojroot_2.0.4   jsonlite_1.9.1   
#> [13] mnormt_2.1.1      cli_3.6.4         rlang_1.1.5       munsell_0.5.1    
#> [17] withr_3.0.2       tools_4.4.2       parallel_4.4.2    tzdb_0.4.0       
#> [21] colorspace_2.1-1  pacman_0.5.1      vctrs_0.6.5       R6_2.6.1         
#> [25] lifecycle_1.0.4   htmlwidgets_1.6.4 pkgconfig_2.0.3   pillar_1.10.1    
#> [29] gtable_0.3.6      glue_1.8.0        xfun_0.51         tidyselect_1.2.1 
#> [33] rstudioapi_0.17.1 farver_2.1.2      htmltools_0.5.8.1 nlme_3.1-167     
#> [37] rmarkdown_2.29    compiler_4.4.2

Bibliografia

Schoot, R. van de, Depaoli, S., King, R., Kramer, B., Märtens, K., Tadesse, M. G., Vannucci, M., Gelman, A., Veen, D., Willemsen, J., et al. (2021). Bayesian statistics and modelling. Nature Reviews Methods Primers, 1(1), 1.