Inferenza bayesiana con MCMC

33. Inferenza bayesiana con MCMC#

In questo capitolo utilizzeremo PyMC, un pacchetto Python per la modellazione statistica bayesiana. Rispetto al capitolo precedente, in cui abbiamo utilizzato le funzioni di scipy.stats, in questo caso useremo i linguaggi di programmazione probabilistici (PPL) per scrivere il nostro campionatore. I PPL consentono di definire i modelli bayesiani utilizzando il codice del PPL e di eseguire l’inferenza bayesiana in modo abbastanza automatizzato, permettendo agli utenti di concentrarsi sulla costruzione dei modelli e meno sui dettagli matematici e computazionali.

L’algoritmo Metropolis, che abbiamo visto precedentemente, consente di generare campioni da distribuzioni di probabilità creando una catena di Markov che ha come distribuzione di equilibrio (o stazionaria) la distribuzione desiderata. Tuttavia, per modelli complessi, questo algoritmo risulta poco efficiente, richiedendo molto tempo per raggiungere una distribuzione stazionaria. Sono stati sviluppati algoritmi Monte Carlo a catena di Markov (MCMC) più efficienti, come il campionatore No-U-Turn (NUTS), gli algoritmi Metropolis-Hastings, Gibbs Sampler e Hamiltonian Monte Carlo, tutti implementati in vari framework per la programmazione probabilistica (PP), noti come “Universal Inference Engines”. PyMC e Stan [CGH+17] sono i due PPL più popolari. Con PyMC è possibile definire modelli complessi, anche con molte migliaia di parametri, utilizzando una sintassi leggibile e intuitiva. In questo capitolo utilizzeremo PyMC per eseguire l’inferenza su una o due proporzioni binomiali utilizzando il campionamento Monte Carlo a catena di Markov.

Iniziamo a caricare i pacchetti necessari.

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import pymc as pm
import scipy.stats as st
from scipy.stats import beta
# PyMC generates a FutureWarning we don't need to deal with yet
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)

print(f"Running on PyMC v{pm.__version__}")

Running on PyMC v5.5.0

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
# Initialize random number generator
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)

plt.style.use("bmh")
plt.rcParams["figure.figsize"] = [10, 6]
plt.rcParams["figure.dpi"] = 100
plt.rcParams["figure.facecolor"] = "white"

sns.set_theme(palette="colorblind")

%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = "svg"

Un esempio pratico di utilizzo del linguaggio probabilistico PyMC consiste nell’eseguire una semplice operazione aritmetica. Iniziamo a sommare due numeri interi in Python.

# adding 2 integers in Python
a = 2
b = 3
c = a + b
print(c)

Facciamo ora la stessa cosa usando PyMC.

# adding 2 random variables in Pymc
with pm.Model() as example:
    a = pm.Normal("a", 2, 0.5)
    b = pm.Normal("b", 3, 0.2)
    c = pm.Deterministic("c", a + b)
    trace_1 = pm.sample()

Auto-assigning NUTS sampler...

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [a, b]

100.00% [8000/8000 00:01<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 20 seconds.

Il risultato che otteniamo è una distribuzione di probabilità di \(c\). Dalla figura vediamo che la moda della distribuzione di \(c\) è uguale al risultato atteso.

az.plot_trace(trace_1)
plt.tight_layout()

_images/736499fa4f298365d2151ba935f3be257b04fbde842cef849cde9a71f79aa351.svg

33.1. Il presidente Trump e l’idrossiclorochina#

Esaminiamo ora un’inferenza, non il caso di una situazione deterministica come quella vista in precedenza. Per fare un esempio concreto, considereremo un set di dati reali che può essere analizzato usando lo schema Beta-Binomiale. Cito dal Washington Post del 7 aprile 2020:

One of the most bizarre and disturbing aspects of President Trump’s nightly press briefings on the coronavirus pandemic is when he turns into a drug salesman. Like a cable TV pitchman hawking ‘male enhancement’ pills, Trump regularly extols the virtues of taking hydroxychloroquine, a drug used to treat malaria and lupus, as a potential ‘game changer’ that just might cure Covid-19.

Esaminiamo le evidenze a supporto dell’ipotesi che l’idrossiclorochina possa essere utile per la cura del Covid-19, ovvero le evidenze che erano disponibili nel momento in cui il Donald Trump ha fatto le affermazioni riportate sopra (in seguito, quest’idea è stata screditata). Tali evidenze sono fornite da uno studio di Gautret et al. [GLP+20]. Il disegno sperimentale di Gautret et al. [GLP+20] comprende, tra le altre cose, il confronto tra una condizione sperimentale e una condizione di controllo. Il confronto importante è tra la proporzione di paziente positivi al virus SARS-CoV-2 nel gruppo sperimentale (a cui è stata somministrata l’idrossiclorochina; 6 su 14) e la proporzione di paziente positivi nel gruppo di controllo (a cui non è stata somministrata l’idrossiclorochina; ovvero 14 su 16). Obiettivo di questo capitolo è mostrare come si possa fare inferenza sui dati di Gautret et al. [GLP+20] usando PyMC. Per semplicità, iniziamo considerando solo il gruppo di controllo.

33.2. Una proporzione#

Sulla base di ciò che è stato detto nel Distribuzioni coniugate, sappiamo che, quando i dati sono rappresentati da una proporzione \(\theta\), e quando utilizziamo una distribuzione a priori Beta per \(\theta\), la distribuzione a posteriori di \(\theta\) è specificata dallo schema beta-binomiale. Se scegliamo, ad esempio, una \(Beta(2, 2)\) quale distribuzione a priori per \(\theta\), il modello diventa:

(33.1)#\[\begin{split} \begin{align} y &\sim Bin(n, \theta) \notag\\ \theta &\sim Beta(2, 2) \end{align} \end{split}\]

dove la prima riga definisce la funzione di verosimiglianza e la seconda riga definisce la distribuzione a priori.

In questo caso specifico, è possibile ottenere una soluzione analitica in forma chiusa per stimare la distribuzione a posteriori e, quindi, non è necessario utilizzare MCMC. Nel caso di \(y\) = 14 e \(n\) = 16, ad esempio, otteniamo una Beta(16, 4). Esaminiamo graficamente le tre funzioni.

n_gridpoints = 1000

alpha_prior = 2
beta_prior = 2

y = 14
ntrials = 16

alpha_post = alpha_prior + y
beta_post = beta_prior + ntrials - y

theta = np.linspace(0, 1, n_gridpoints)

_, ax = plt.subplots()

# prior
fx_prior = st.beta.pdf(theta, alpha_prior, beta_prior)
_ = ax.plot(theta, fx_prior)

# likelihood
like_unstandardized = st.binom.pmf(y, ntrials, theta) 
arbitrary_constant = 18
like = like_unstandardized * arbitrary_constant
_ = ax.plot(theta, like)

# posterior
fx_post = st.beta.pdf(theta, alpha_post, beta_post)
_ = ax.plot(theta, fx_post, color="black")

ax.set_xlabel("$\\theta$")
ax.set_ylabel("Densità")

Text(0, 0.5, 'Densità')

_images/6ffcb4a1efbced3dd528f640a5078881bd76c93978dee3a1cd78c4227290c1b0.svg

Anche se in questo caso è possibile ottenere una soluzione analitica per la distribuzione a posteriori, nella maggior parte dei modelli di inferenza bayesiana questo non è possibile, ed è pertanto necessario fare ricorso a tecniche di approssimazione numerica come MCMC. In questo esempio, applicheremo MCMC a un caso in cui conosciamo già la soluzione analitica, in modo da poter confrontare i risultati ottenuti tramite la soluzione analitica e quelli ottenuti tramite l’approssimazione numerica.

33.3. Dedurre una proporzione con PyMC#

Ora eseguiremo la stessa analisi che abbiamo svolto in precedenza utilizzando il metodo numerico Markov Chain Monte Carlo. Supponiamo di avere già installato PyMC. Una volta installato, dobbiamo successivamente importare le librerie necessarie, che includono Matplotlib, Numpy, Scipy, Arviz e lo stesso PyMC.

Vediamo ora come specificare il modello beta-binomiale mediante PyMC. Per svolgere l’analisi mediante PyMC è necessario prima specificare la struttura del modello bayesiano e poi eseguire il campionamento dalla distribuzione a posteriori. Esaminiamo questi due passaggi per l’esempio presente.

33.4. Dati#

I dati sono i seguenti e indicano che la verosimiglianza corrisponderà al modello di probabilità Binomiale.

ntrials = 16
y = 14

33.5. Distribuzione a priori#

Useremo una distribuzione Beta per rappresentare tutti possibili valori della probabilità di successo, che chiameremo \(\theta\). Decidiamo di rappresentare le nostre credenze a priori rispetto a \(\theta\) usando una distribuzione a priori debolmente informativa, ovvero una distribuzione a priori che non introduce alcun bias nella stima della distribuzione a posteriori, ma si limita a escludere le possibilità estreme (\(\theta\) = 0, \(\theta\) = 1) assegnando alle rimanenti possibilità la maggiore incertezza possibile. Ciò può essere ottenuto scegliendo una Beta(\(\alpha\) = 2, \(\beta\) = 2). Definiamo i parametri della distribuzione Beta.

alpha_prior = 2
beta_prior = 2

33.5.1. Specificare il modello#

Per specificare il modello con PyMC si usa context with di Python. Ciò consente di assegnare tutti i parametri, gli argomenti e i valori iniziali a un’istanza pymc.Model (che qui denominiamo bb_model).

bb_model = pm.Model()

with bb_model:
    # Prior
    theta = pm.Beta("theta", alpha=alpha_prior, beta=beta_prior)
    # Likelihood
    obs = pm.Binomial("obs", p=theta, n=ntrials, observed=y)

La prima linea di codice

bb_model = pm.Model()

crea un nuovo oggetto di classe Model che è un contenitore per le variabili casuali del modello. Dopo l’istanziazione del modello, la successiva specificazione delle componenti del modello viene eseguita all’interno di un’istruzione with:.

Questo crea un context manager, con il nostro bb_model come contesto, che include tutte le istruzioni fino alla fine del blocco indentato. Ciò significa che tutti gli oggetti PyMC introdotti nel blocco di codice indentato sotto l’istruzione with: vengono aggiunti al modello. In assenza di questo idioma del context manager, saremmo costretti ad associare manualmente ciascuna delle variabili a bb_model.

Con la chiamata del costruttore pm.Beta si crea una variabile casuale da usare come distribuzione a priori. Il primo argomento è sempre il nome della variabile casuale, che dovrebbe corrispondere al nome della variabile Python assegnata, poiché a volte viene utilizzato per recuperare la variabile dal modello per riassumere l’output. I restanti argomenti richiesti per un oggetto stocastico sono i parametri, in questo caso alpha e beta, a cui assegniamo i valori degli iperparametri del modello. Nel caso presente, usiamo i valori alpha_prior e beta_prior come parametri della distribuzione Beta a priori.

La riga finale del modello definisce obs, la distribuzione campionaria della variabile di esito nel set di dati. Questo è un caso speciale di variabile stocastica chiamata observed stochastic, la quale rappresenta la verosimiglianza dei dati del modello. Tale variabile stocastica è identica alle altre variabili stocastiche, tranne per l’argomento obs, il quale indica che i valori per questa variabile sono stati osservati e non devono essere modificati dal modello. I dati possono essere passati sotto forma di oggetto ndarray o DataFrame. Nel caso presente, per la funzione di verosimiglianza Binomiale specifichiamo il parametro p=theta, il numero di prove n=ntrials e il numero di successi observed=y.

Si noti che, a differenza delle distribuzioni a priori del modello, i parametri per la distribuzione Beta di obs non sono valori fissi, ma piuttosto corrispondono all’oggetto stocastico theta. Questo crea una relazione genitore-figlio tra la verosimiglianza e questa variabile.

33.5.3. Analisi a posteriori#

Esaminiamo l’accettanza.

az.plot_posterior(
    idata,
    group="sample_stats",
    var_names="acceptance_rate",
    hdi_prob="hide",
    kind="hist",
)

<Axes: title={'center': 'acceptance_rate'}>

_images/51739573a90b4e4b620b8d2abc75e0e0681521aef68151962409f618841f6d1a.svg

Per combinare catene e iterazioni, utilizziamo la funzione arviz.extract(). Creiamo un oggetto DataArray di Xarray chiamato post.

post = az.extract(idata)

Le stime a posteriori di \(\theta\) sono accessibili nel modo seguente.

post.mean()

Oppure nel modo seguente:

post["theta"].mean()

<xarray.DataArray 'theta' ()>
array(0.79920394)

Possiamo calcolare la mediana a posteriori di \(\theta\) in questo modo:

post.median()

La deviazione standard della stima a posteriori è

np.std(post["theta"])

<xarray.DataArray 'theta' ()>
array(0.08857597)

Un sommario della distribuzione a posteriori si ottiene con il metodo az.summary().

az.summary(idata, hdi_prob=0.95, round_to=3)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
theta	0.799	0.089	0.629	0.955	0.001	0.001	3427.614	3845.485	1.001

Si ottiene così l’intervallo di credibilità a più alta densità a posteriori (HPD) al 95%. Questo intervallo ci informa sul fatto che, a posteriori, possiamo essere certi al 95%, che il vero valore del parametro \(\theta\) sia contenuto nell’intervallo [0.625, 0.95].

Un risultato praticamente identico si ottiene usando l’algoritmo di Metropolis, anziché NUTS.

az.summary(metropolis_idata, hdi_prob=0.95, round_to=3)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
theta	0.798	0.087	0.626	0.954	0.002	0.001	1802.807	2094.776	1.003

L’output grafico dell’analisi è riportato nel grafico seguente.

p_post = post["theta"]

# Posterior: Beta(alpha + y, beta + n - y)
alpha_post = alpha_prior + y
beta_post = beta_prior + ntrials - y

plt.hist(
    p_post,
    bins=50,
    histtype="step",
    density=True,
    label="A posteriori (MCMC)",
    color="red",
)

# Plot the analytic prior and posterior beta distributions
x = np.linspace(0, 1, 100)
plt.plot(x, beta.pdf(x, alpha_prior, beta_prior), "--", label="A priori", color="blue")
plt.plot(
    x,
    beta.pdf(x, alpha_post, beta_post),
    label="A posteriori (analitico)",
    color="green",
)

# Update the graph labels
plt.legend(title=" ", loc="best")
plt.xlabel("$\\theta$, Proporzione di successi")
plt.ylabel("Densità")

Text(0, 0.5, 'Densità')

_images/b74d3e582f279da8ae7f755731c605532373abc51fc91a1401ed0a6f46ba74e8.svg

In questo esempio, con 8,000 campioni, la convergenza dell’algoritmo NUTS è estremamente buona. L’istogramma segue da vicino la distribuzione a posteriori calcolata analiticamente, come da previsione.

Un tracciato della catena di Markov illustra questa esplorazione rappresentando il valore \(\theta\) sulle ordinate e l’indice progressivo di in ogni iterazione sull’ascissa. Il trace plot è estremamente utile per valutare la convergenza di un algoritmo MCMC e se è necessario escludere un periodo di campioni iniziali (noto come burn-in). Per produrre la traccia chiamiamo semplicemente az.plot_trace() con la variabile idata:

az.plot_trace(idata)
plt.show()

_images/c391bdc39522bc9b61e14b3e45caf14f54ce43984709644db8b1d8fa31abe2ef.svg

La traccia descrive il comportamento longitudinale delle catene di Markov. Possiamo usare il metodo az.plot_trace() per visualizzare sia la traccia sia una stima della densità del kernel (KDE) dell’istogramma delle stime a posteriori, ovvero, dei valori che le catene MCMC visitano lungo il loro percorso, ignorando l’ordine di queste visite. Da notare come la stima di KDE della convinzione a posteriori nella probabilità di efficacia dell’idrossiclorochina riflette sia la convinzione a priori di 𝜎=0.22 che i nostri dati con una correttezza campionaria di 𝜎=0.09. Inoltre possiamo vedere che la procedura di campionamento MCMC è “convergente alla distribuzione” poiché la serie di campionamento sembra stazionaria.

Svolgendo un’analisi bayesiana simile a questa, Gautret et al. [GLP+20] hanno trovato che gli intervalli di credibilità del gruppo di controllo e del gruppo sperimentale non si sovrappongono. Questo fatto viene interpretato dicendo che il parametro \(\theta\) è diverso nei due gruppi. Sulla base di queste evidenza, Gautret et al. [GLP+20] hanno concluso, con un grado di certezza soggettiva del 95%, che nel gruppo sperimentale vi è una probabilità più bassa di risultare positivi al SARS-CoV-2 rispetto al gruppo di controllo. In altri termini, l’analisi statistica condotta da Gautret et al. [GLP+20] suggerisce che l’idrossiclorochina è una terapia efficace per il Covid-19.

33.6. La critica di Hulme (2020)#

Un articolo pubblicato da Hulme et al. [HWD+20] si è posto il problema di rianalizzare i dati di Gautret et al. [GLP+20].1 Tra gli autori di questo articolo figura anche Eric-Jan Wagenmakers, uno psicologo molto conosciuto per i suoi contributi metodologici. Hulme et al. [HWD+20] osservano che, nelle loro analisi statistiche, Gautret et al. [GLP+20] hanno escluso alcuni dati. Nel gruppo sperimentale, infatti, vi erano alcuni pazienti i quali, anziché migliorare, sono in realtà peggiorati. L’analisi statistica di Gautret et al. [GLP+20] ha escluso i dati di questi pazienti. Se consideriamo tutti i pazienti – non solo quelli selezionati da Gautret et al. [GLP+20] – la situazione diventa la seguente:

gruppo sperimentale: 10 positivi su 18;
gruppo di controllo: 14 positivi su 16.

L’analisi dei dati proposta da [HWD+20] richiede l’uso di alcuni strumenti statistici che, in queste dispense, non verranno discussi. Ma possiamo giungere alle stesse conclusioni raggiunte da questi ricercatori anche usando le procedure statistiche descritte nel Paragrafo successivo.

33.8. Commenti e considerazioni finali#

La ricerca di Gautret et al. [GLP+20] include altre informazioni e altre analisi statistiche che non sono state qui considerate. Tuttavia, la semplice analisi statistica che abbiamo qui descritto è stata in grado di replicare le conclusioni a cui sono giunti (per altra via) Hulme et al. [HWD+20].

33.9. Watermark#

%load_ext watermark
%watermark -n -u -v -iv -w

Last updated: Tue Jul 04 2023

Python implementation: CPython
Python version       : 3.11.3
IPython version      : 8.12.0

seaborn   : 0.12.2
arviz     : 0.15.1
matplotlib: 3.7.1
pymc      : 5.5.0
numpy     : 1.24.3
pandas    : 1.5.3
scipy     : 1.10.1

Watermark: 2.3.1

1: Si veda https://osf.io/5dgmx/.