R Academy · Lezione

Controlli predittivi posteriori

Valuti l'adattamento del modello confrontando le distribuzioni dei dati simulati e osservati

Lezione 4 di 413 passaggi

Controlli predittivi posteriori è una lezione R Academy gratuita su CoddyKit. Questa è la lezione 4 di 4. Puoi leggere la lezione completa qui gratuitamente — poi esercitati direttamente nel browser con un editor di codice integrato e un tutor IA disponibile 24/7. Fa parte del percorso di apprendimento R Academy, e i tuoi progressi si sincronizzano tra il web e l'app CoddyKit. Il corso R Academy include 4 lezioni in totale.

Che cosa sono i controlli predittivi a posteriori?

Dopo aver adattato un modello bayesiano, è necessario chiedersi: questo modello genera dati simili a quelli osservati? I controlli predittivi a posteriori (PPC) rispondono simulando dataset replicati yrep dalla distribuzione a posteriori e confrontandoli graficamente con i dati osservati y.

Estrarre i campioni a posteriori

extract(fit, pars='mu') restituisce una lista denominata; $mu è un vettore numerico di tutti i campioni successivi al warmup per quel parametro. Con 4 catene e 1000 iterazioni dopo il warmup si ottengono 4000 campioni.

# library(rstan)
# mu_samples    <- extract(fit, pars = 'mu')$mu
# sigma_samples <- extract(fit, pars = 'sigma')$sigma
#
# cat('Samples drawn:', length(mu_samples), '
')
# cat('Posterior mean of mu:', mean(mu_samples), '
')
# cat('90% CI:', quantile(mu_samples, c(0.05, 0.95)), '
')

Generare yrep dai campioni a posteriori

Per ogni estrazione a posteriori (mu_s, sigma_s), simuli un dataset replicato delle stesse dimensioni dei dati originali. Memorizzi questi dataset in una matrice yrep, in cui ogni riga rappresenta un dataset simulato.

# y <- c(2.1, 1.8, 2.4, 1.9, 2.3, 2.0, 1.7, 2.2)
# n_obs <- length(y)
# S     <- length(mu_samples)  # 4000 posterior draws
#
# yrep <- matrix(NA, nrow = S, ncol = n_obs)
# for (s in seq_len(S)) {
#   yrep[s, ] <- rnorm(n_obs, mean = mu_samples[s], sd = sigma_samples[s])
# }
# dim(yrep)  # [4000, 8]

ppc_dens_overlay() — Confronto delle densità

bayesplot::ppc_dens_overlay(y, yrep[1:50,]) sovrappone la densità kernel dei dati osservati (linea scura) alle densità di 50 dataset simulati scelti casualmente (linee chiare). Un buon adattamento del modello significa che la linea scura si trova all'interno dell'insieme delle linee chiare.

# library(bayesplot)
#
# ppc_dens_overlay(y, yrep[1:50, ])
#
# Interpretation:
# - Dark line (y_obs) surrounded by light lines (yrep): good fit
# - Dark line systematically outside the cloud: model misfit
# - Light lines much wider than dark: overdispersed model
# - Light lines much narrower than dark: underdispersed model

ppc_stat() — Controllo della statistica test

ppc_stat(y, yrep, stat = 'mean') mostra un istogramma della statistica test (ad esempio, la media) calcolata su ciascun dataset simulato, con una linea verticale in corrispondenza della statistica osservata. Se il valore osservato ricade nella parte centrale dell'istogramma, il modello rappresenta adeguatamente quell'aspetto dei dati.

# library(bayesplot)
#
# ppc_stat(y, yrep, stat = 'mean')   # does model capture the mean?
# ppc_stat(y, yrep, stat = 'sd')     # does model capture spread?
# ppc_stat(y, yrep, stat = 'max')    # does model capture extremes?
#
# If observed stat is in the tail of the histogram,
# the model fails to reproduce that statistic.

p-value bayesiano

Il p-value bayesiano (p-value predittivo a posteriori) è la proporzione di dataset simulati la cui statistica test è più estrema del valore osservato. Valori vicini a 0.5 indicano una buona calibrazione; valori vicini a 0 o 1 indicano un cattivo adattamento del modello per quella statistica.

# Bayesian p-value for the mean:
# obs_mean <- mean(y)
# rep_means <- apply(yrep, 1, mean)
# pval <- mean(rep_means >= obs_mean)
# cat('Bayesian p-value (mean):', round(pval, 3), '
')
# # 0.5 is perfect; < 0.05 or > 0.95 suggests misfit

Altre funzioni PPC di bayesplot

bayesplot offre numerose visualizzazioni PPC oltre alla sovrapposizione delle densità:

ppc_hist(y, yrep[1:8,]) — griglia di istogrammi
ppc_scatter_avg(y, yrep) — grafico a dispersione dei dati osservati rispetto alla media di yrep
ppc_intervals(y, yrep) — intervalli di incertezza attorno a ciascuna osservazione
ppc_rootogram(y, yrep) — per dati di conteggio

# library(bayesplot)
#
# # Grid of 8 simulated histograms vs the observed
# ppc_hist(y, yrep[1:8, ])
#
# # Scatter: y_obs (x) vs mean of yrep (y) — should hug diagonal
# ppc_scatter_avg(y, yrep)
#
# # 50% and 90% posterior predictive intervals around each y_i
# ppc_intervals(y, yrep)

Interpretare i grafici PPC — Cattivo adattamento del modello

Pattern comuni di cattivo adattamento e relative cause:

yrep troppo ampia — la prior è troppo diffusa o il modello presenta una dispersione eccessiva
yrep spostata — famiglia di verosimiglianza errata (ad esempio, normale per dati asimmetrici)
yrep non riproduce la multimodalità — è necessario un modello a miscela
yrep non riproduce i valori estremi — è necessaria una distribuzione con code pesanti

# Example: if data has a long right tail but yrep does not,
# consider switching:
# y ~ normal(mu, sigma)  =>  y ~ student_t(nu, mu, sigma)
#
# Or for count data:
# y ~ poisson(lambda)  =>  y ~ neg_binomial_2(mu, phi)  (overdispersion)
cat('PPCs guide model improvement by revealing specific failure modes
')

PPC nel blocco Stan model

È possibile generare yrep direttamente in Stan usando il blocco generated quantities. In questo modo si evita di estrarre nuovamente i parametri in R e il risultato è equivalente dal punto di vista computazionale.

# Stan model with generated quantities:
# '
# generated quantities {
#   array[N] real y_rep;
#   for (n in 1:N) {
#     y_rep[n] = normal_rng(mu, sigma);
#   }
# }
# '
# Then extract in R:
# yrep <- extract(fit, pars = 'y_rep')$y_rep  # [S, N] matrix

Convalida incrociata leave-one-out

Oltre ai PPC, loo::loo(fit) calcola la convalida incrociata leave-one-out per confrontare modelli concorrenti. Si preferisce il modello con ELPD (densità predittiva logaritmica attesa) più alto. Usi loo::loo_compare(loo1, loo2) per classificare i modelli.

# library(loo)
# loo1 <- loo(fit1)  # normal model
# loo2 <- loo(fit2)  # student-t model
#
# comparison <- loo_compare(loo1, loo2)
# print(comparison)
#
# Model with elpd_diff > 0 is preferred
# se_diff > |elpd_diff| means difference is not reliable

Buone pratiche per i PPC

Segua queste pratiche per eseguire controlli predittivi a posteriori rigorosi:

Inizi sempre con ppc_dens_overlay() come controllo globale di plausibilità
Prosegua con statistiche specifiche del dominio (ppc_stat()) pertinenti agli obiettivi dell'analisi
Usi almeno 50 estrazioni yrep per i controlli visivi; tutte le 4000 per i p-value bayesiani
I PPC non superati guidano il miglioramento del modello: sono diagnostici, non rappresentano un fallimento

# Workflow:
# 1. Fit model -> extract() -> generate yrep matrix
# 2. ppc_dens_overlay(y, yrep[1:50,])  -- visual global check
# 3. ppc_stat(y, yrep, stat='mean')    -- check mean
# 4. ppc_stat(y, yrep, stat='sd')      -- check spread
# 5. ppc_stat(y, yrep, stat='max')     -- check tails
# 6. If misfit found -> revise model -> refit -> re-check

Verifica rapida: interpretazione del p-value bayesiano

Che cosa indica per il modello un p-value bayesiano pari a 0.03 per la statistica del massimo?

Rispondi nella lezione interattiva

Riepilogo dei controlli predittivi a posteriori

Flusso di lavoro PPC in RStan e bayesplot:

Estrarre i campioni: extract(fit, pars='mu')$mu
Generare yrep: eseguire un ciclo sulle estrazioni a posteriori chiamando rnorm(n, mu_s, sigma_s)
Controllo globale: ppc_dens_overlay(y, yrep[1:50,])
Controlli delle statistiche: ppc_stat(y, yrep, stat='mean')
p-value bayesiano: mean(apply(yrep,1,stat) >= stat(y)) — un valore vicino a 0.5 è positivo
Per confrontare i modelli, usi loo::loo_compare()

Gratis per iniziare

Impara R con un tutor IA — gratis

Scrivi ed esegui vero codice nel tuo browser, ricevi aiuto istantaneo da un tutor IA disponibile 24/7, e riprendi da dove hai lasciato sul web o nell'app.

Corsi: 43
Lezioni: 159

Inizia a imparare gratis Visualizza il corso completo

Domande Frequenti

La lezione «Controlli predittivi posteriori» è gratuita?

Sì — il testo completo di «Controlli predittivi posteriori» è gratuito qui sul web. Per esercitarvi in modo interattivo (un editor di codice integrato e un tutor IA 24/7) e sbloccare il resto del corso R Academy, passa a CoddyKit PRO. Il corso R Academy include 4 lezioni in totale.

Cosa imparerò in «Controlli predittivi posteriori»?

Valuti l'adattamento del modello confrontando le distribuzioni dei dati simulati e osservati Eserciti R Academy con codice pratico che esegui direttamente nel browser, e un tutor IA 24/7 risponde alle tue domande mentre lavori sulla lezione.

Ho bisogno di esperienza per iniziare R Academy?

Non è richiesta alcuna esperienza precedente. R Academy su CoddyKit è strutturato per principianti e studenti avanzati, quindi puoi iniziare da qui o dall'inizio e procedere al tuo ritmo. Questa è la lezione 4 di 4.

Quanto tempo richiede la lezione «Controlli predittivi posteriori»?

La maggior parte delle lezioni CoddyKit richiede circa 5–10 minuti. Ogni lezione è breve e interattiva, quindi fai progressi costanti e riprendi esattamente da dove hai lasciato su web e app.

Posso scrivere ed eseguire codice in questa lezione R Academy?

Sì. Ogni lezione R Academy include un editor di codice integrato, quindi scrivi ed esegui codice reale direttamente nel tuo browser e ricevi feedback istantaneo dall'IA — nessuna configurazione locale necessaria.

Tutte le lezioni di questo corso

← Torna a R Academy