Inferenza controfattuale

Inferenza controfattuale#

Il ragionamento causale e il pensiero controfattuale sono argomenti davvero interessanti, ma anche complessi. Tuttavia, possiamo progredire nella comprensione di queste idee attraverso esempi relativamente semplici. Questo notebook si concentra sui concetti e sull’implementazione pratica del ragionamento causale Bayesiano utilizzando PyMC.

Per fare questo useremo un esempio importante ma triste relativo al calcolo delle morti in eccesso dovute al COVID-19. A tal proposito, le idee presentate in questo notebook si sovrappongono fortemente a quelle proposte da Brodersen et al. [BGK+15]. In pratica, cercheremo di stimare il numero di ‘morti in eccesso’ dall’inizio della pandemia di COVID-19, utilizzando dati provenienti dall’Inghilterra e dal Galles. Le morti in eccesso sono definite come:

\[ \text{Morti in eccesso} = \text{Morti osservate} - \text{Morti attese (controfattuali)} \]

Fare un’affermazione sulle morti in eccesso richiede un ragionamento causale e controfattuale. Mentre il numero di morti segnalate non è altro che una misura (sicuramente rumorosa) di un fatto realmente osservabile nel mondo, le morti attese non sono misurabili perché non si realizzano mai nella nostra linea temporale. In altre parole, l’idea delle morti attese è un esperimento controfattuale in cui ci poniamo la domanda “Cosa sarebbe/succederà se?”.

Questo notebook rappresenta una versione semplificata del notebook originale disponibile sul sito ufficiale di PyMC. La differenza fondamentale risiede nel numero di variabili predittive utilizzate. Nel notebook originale, il modello prende in considerazione tre predittori per il numero di morti nell’era pre-covid: l’anno, il mese e la temperatura. In questo tutorial, per ragioni di semplicità, la variabile del mese è stata omessa. Nonostante questa semplificazione, i risultati rimangono in gran parte analoghi a quelli presentati nel tutorial originale.

Preparazione del Notebook#

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import seaborn as sns
import scipy as sc
import statistics as st
import pymc as pm
import pymc.sampling_jax
import pytensor.tensor as pt
import xarray as xr
import arviz as az
import warnings

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=Warning)

/Users/corrado/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

%config InlineBackend.figure_format = 'retina'
RANDOM_SEED = 12345
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")
sns.set_theme(palette="colorblind")

Utilizzeremo qui una serie di funzioni che consento di genereare i grafici delle quantità di interesse con maggiore facilità. Per i nostri scopi, non è importante capire come sono state scritte tali funzioni.

Show code cell content Hide code cell content

def ZeroSumNormal(name, *, sigma=None, active_dims=None, dims, model=None):
    model = pm.modelcontext(model=model)

    if isinstance(dims, str):
        dims = [dims]

    if isinstance(active_dims, str):
        active_dims = [active_dims]

    if active_dims is None:
        active_dims = dims[-1]

    def extend_axis(value, axis):
        n_out = value.shape[axis] + 1
        sum_vals = value.sum(axis, keepdims=True)
        norm = sum_vals / (pt.sqrt(n_out) + n_out)
        fill_val = norm - sum_vals / pt.sqrt(n_out)
        out = pt.concatenate([value, fill_val], axis=axis)
        return out - norm

    dims_reduced = []
    active_axes = []
    for i, dim in enumerate(dims):
        if dim in active_dims:
            active_axes.append(i)
            dim_name = f"{dim}_reduced"
            if name not in model.coords:
                model.add_coord(dim_name, length=len(model.coords[dim]) - 1, mutable=False)
            dims_reduced.append(dim_name)
        else:
            dims_reduced.append(dim)

    raw = pm.Normal(f"{name}_raw", sigma=sigma, dims=dims_reduced)
    for axis in active_axes:
        raw = extend_axis(raw, axis)
    return pm.Deterministic(name, raw, dims=dims)


def format_x_axis(ax, minor=False):
    # major ticks
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y %b"))
    ax.xaxis.set_major_locator(mdates.YearLocator())
    ax.grid(which="major", linestyle="-", axis="x")
    # minor ticks
    if minor:
        ax.xaxis.set_minor_formatter(mdates.DateFormatter("%Y %b"))
        ax.xaxis.set_minor_locator(mdates.MonthLocator())
        ax.grid(which="minor", linestyle=":", axis="x")
    # rotate labels
    for label in ax.get_xticklabels(which="both"):
        label.set(rotation=70, horizontalalignment="right")


def plot_xY(x, Y, ax):
    quantiles = Y.quantile((0.025, 0.25, 0.5, 0.75, 0.975), dim=("chain", "draw")).transpose()

    az.plot_hdi(
        x,
        hdi_data=quantiles.sel(quantile=[0.025, 0.975]),
        fill_kwargs={"alpha": 0.25},
        smooth=False,
        ax=ax,
    )
    az.plot_hdi(
        x,
        hdi_data=quantiles.sel(quantile=[0.25, 0.75]),
        fill_kwargs={"alpha": 0.5},
        smooth=False,
        ax=ax,
    )
    ax.plot(x, quantiles.sel(quantile=0.5), color="C1", lw=3)


# default figure sizes
figsize = (10, 5)

Calcolo delle morti in eccesso dovute al COVID-19#

Esamineremo i dati resi disponibili sul sito ufficiale di PyMC. Sebbene sia possibile condurre l’analisi utilizzando dati italiani, opteremo per i dati forniti sul sito per motivi di semplicità.

df = pd.read_csv("../data/deaths_and_temps_england_wales.csv")
df.head()

	date	temp	deaths	year	month	t	pre
0	2006-01-01	3.8	49124	2006	1	0	True
1	2006-02-01	3.4	42664	2006	2	1	True
2	2006-03-01	3.9	49207	2006	3	2	True
3	2006-04-01	7.4	40645	2006	4	3	True
4	2006-05-01	10.7	42425	2006	5	4	True

df["date"] = pd.to_datetime(df["date"])
df = df.set_index("date")

Lo scopo delle due righe precedenti è convertire una colonna del DataFrame df in un formato data/tempo e impostarla come indice del DataFrame. Vediamo in dettaglio cosa fanno:

df["date"] = pd.to_datetime(df["date"]): Questa riga prende la colonna “date” del DataFrame df e la converte in un formato data/tempo utilizzando la funzione pd.to_datetime(). Questa funzione è molto utile quando si lavora con serie temporali, poiché facilita la manipolazione e l’analisi dei dati. Dopo questa operazione, la colonna “date” sarà di tipo datetime64.
df = df.set_index("date"): Questa riga imposta la colonna “date” come indice del DataFrame. Gli indici in un DataFrame pandas servono per etichettare le righe, facilitando così il filtraggio, la selezione e altre operazioni. Impostare una colonna di date come indice è particolarmente utile quando si lavora con serie temporali, in quanto consente di eseguire facilmente operazioni come il riordinamento dei dati in base alla data o l’esecuzione di calcoli che coinvolgono intervalli di tempo.

Ora dividiamo il DataFrame df in due DataFrame distinti, corrispondenti al periodo pre e post insorgenza del COVID-19.

pre = df[df.index < "2020"]
post = df[df.index >= "2020"]

pre.head()

	temp	deaths	year	month	t	pre
date
2006-01-01	3.8	49124	2006	1	0	True
2006-02-01	3.4	42664	2006	2	1	True
2006-03-01	3.9	49207	2006	3	2	True
2006-04-01	7.4	40645	2006	4	3	True
2006-05-01	10.7	42425	2006	5	4	True

La variabile t è un indice che identifica la sequenza temporale delle misurazioni: 12 mesi per 14 anni.

14*12

pre["t"]

date
2006-01-01      0
2006-02-01      1
2006-03-01      2
2006-04-01      3
2006-05-01      4
             ... 
2019-08-01    163
2019-09-01    164
2019-10-01    165
2019-11-01    166
2019-12-01    167
Name: t, Length: 168, dtype: int64

Generiamo un grafico con il numero di morti in funzione del tempo, distinguendo tra prima e dopo l’insorgenza del COVID-19.

ax = sns.lineplot(data=df, x="date", y="deaths", hue="pre")
format_x_axis(ax)

../_images/976d03e49d06d3f2d91408fad58e5c904d857c3b2fda7414b285a47624a48c5e.png

Esaminiamo i dati relativi al periodo pre-COVID-19. La serie temporale evidenzia una tendenza generale all’incremento del numero di decessi nel corso del tempo, oltre a un marcato effetto stagionale. Quest’ultimo è attribuibile alle fluttuazioni della temperatura.

annual_deaths = pd.DataFrame(pre.groupby("year")["deaths"].sum()).reset_index()
sns.regplot(x="year", y="deaths", data=annual_deaths);

../_images/7707ecd0a50bba7261981b2b0564d0fd767968b2b79cf956680f69871133f75c.png

fig, ax = plt.subplots(1, 2, figsize=figsize)
sns.regplot(x="temp", y="deaths", data=pre, scatter_kws={"s": 40}, order=1, ax=ax[0])
ax[0].set(title="Linear fit (pre COVID-19 data)")
sns.regplot(x="temp", y="deaths", data=pre, scatter_kws={"s": 40}, order=2, ax=ax[1])
ax[1].set(title="Quadratic fit (pre COVID-19 data)");

../_images/48b070cc14a8fe020bc354c1c077aa85bd351af03c45e947c52dff36886d20c4.png

Per una descrizione più accurata, l’effetto della temperatura potrebbe essere modellato attraverso una funzione quadratica. Tuttavia, per semplicità, opteremo per una rappresentazione lineare. Nei dati relativi al periodo pre-COVID-19, un aumento di un grado nella temperatura è associato a una diminuzione media di 764 decessi.

slope, intercept = np.polyfit(pre["temp"], pre["deaths"], 1)
print(f"{slope:.0f} deaths/degree")

-764 deaths/degree

Per il periodo pre-COVID-19, proponiamo dunque un modello di regressione multipla che considera due predittori: l’anno e la temperatura. L’anno è incluso nel modello per tener conto della tendenza lineare all’aumento del numero di morti nel tempo, mentre la temperatura è considerata in quanto fattore significativo.

Il modello matematico è espresso come segue:

\[ \text{deaths}_t \sim \text{TruncatedNormal}(\mu_t, \sigma), \quad \mu_t = \alpha + \beta_1 \cdot \text{time}_t + \beta_2 \cdot \text{temp}_t \]

In questa formula, \(\text{deaths}_t\) rappresenta il numero di morti al tempo \(t\) e viene modellato come una variabile casuale con distribuzione TruncatedNormal avente media \(\mu_t\) e deviazione standard \(\sigma\). La media \(\mu_t\) è una funzione lineare del tempo (\(\text{time}_t\)) e della temperatura (\(\text{temp}_t\)), dove \(\alpha\) è il termine di intercetta e \(\beta_1\) e \(\beta_2\) sono i coefficienti di regressione.

In un’ottica bayesiana, ciascun parametro del modello è trattato come una variabile casuale dotata di una distribuzione a priori. Utilizzeremo le seguenti distribuzioni a priori:

\(\alpha\) (Intercept): distribuzione a priori Normale con media 40.000 e deviazione standard 10.000.
\(\beta_1\) (Coefficiente temporale): distribuzione a priori TruncatedNormal con media 0 e deviazione standard 50, limitata a valori positivi.
\(\beta_2\) (Coefficiente per la temperatura): distribuzione a priori Normale con media 0 e deviazione standard 200.
\(\sigma\): distribuzione a priori HalfNormal con deviazione standard 2.000.

L’obiettivo è quindi quello di stimare la distribuzione a posteriori di questi parametri alla luce dei dati osservati, permettendo così inferenze probabilistiche sul numero di morti in relazione al tempo e alla temperatura.

with pm.Model() as model:
    # observed predictors and outcome
    time = pm.MutableData("time", pre["t"].to_numpy(), dims="t")
    temp = pm.MutableData("temp", pre["temp"].to_numpy(), dims="t")
    deaths = pm.MutableData("deaths", pre["deaths"].to_numpy(), dims="t")

    # priors
    alpha = pm.Normal("intercept", 40_000, 10_000)
    beta_1 = pm.TruncatedNormal("linear trend", 0, 50, lower=0)
    beta_2 = pm.Normal("temp coeff", 0, 200)

    # the actual linear model
    mu = pm.Deterministic(
        "mu",
        alpha + (beta_1 * time) + (beta_2 * temp),
        dims="t",
    )
    sigma = pm.HalfNormal("sigma", 2_000)
    # likelihood
    pm.TruncatedNormal("obs", mu=mu, sigma=sigma, lower=0, observed=deaths, dims="t")

L’utilizzo di una distribuzione TruncatedNormal con un limite inferiore di zero per il numero di morti è una scelta modellistica che riflette alcune considerazioni pratiche e teoriche. Il numero di morti non può essere negativo, quindi un limite inferiore di zero è una scelta naturale per rispettare questa restrizione. Inoltre, le distribuzioni troncate sono spesso più facili da gestire computazionalmente rispetto ad altre distribuzioni asimmetriche che potrebbero essere usate per modellare conteggi come il numero di morti (ad esempio, la distribuzione di Poisson o la distribuzione binomiale negativa). Troncare una distribuzione normale conserva molte delle sue proprietà matematiche desiderabili. Per queste ragioni, la distribuzione TruncatedNormal con un limite inferiore di zero è una scelta ragionevole quando si modellano variabili che non possono assumere valori negativi, come nel caso del numero di morti.

Eseguiamo il campionamento usando i dati pre-COVID-19.

with model:
    idata = pm.sample(random_seed=rng)

Show code cell output Hide code cell output

Auto-assigning NUTS sampler...

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [intercept, linear trend, temp coeff, sigma]

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[15], line 2
with model:
----> 2     idata = pm.sample(random_seed=rng)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/mcmc.py:802, in sample(draws, tune, chains, cores, random_seed, progressbar, step, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, model, **kwargs)
_print_step_hierarchy(step)
try:
--> 802     _mp_sample(**sample_args, **parallel_args)
except pickle.PickleError:
   _log.warning("Could not pickle model, sampling singlethreaded.")

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/mcmc.py:1177, in _mp_sample(draws, tune, step, chains, cores, random_seed, start, progressbar, traces, model, callback, mp_ctx, **kwargs)
# We did draws += tune in pm.sample
draws -= tune
-> 1177 sampler = ps.ParallelSampler(
   draws=draws,
   tune=tune,
   chains=chains,
   cores=cores,
   seeds=random_seed,
   start_points=start,
   step_method=step,
   progressbar=progressbar,
   mp_ctx=mp_ctx,
)
try:
   try:

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/parallel.py:402, in ParallelSampler.__init__(self, draws, tune, chains, cores, seeds, start_points, step_method, progressbar, mp_ctx)
if mp_ctx.get_start_method() != "fork":
   step_method_pickled = cloudpickle.dumps(step_method, protocol=-1)
--> 402 self._samplers = [
   ProcessAdapter(
       draws,
       tune,
       step_method,
       step_method_pickled,
       chain,
       seed,
       start,
       mp_ctx,
   )
   for chain, seed, start in zip(range(chains), seeds, start_points)
]
self._inactive = self._samplers.copy()
self._finished: List[ProcessAdapter] = []

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/parallel.py:403, in <listcomp>(.0)
if mp_ctx.get_start_method() != "fork":
   step_method_pickled = cloudpickle.dumps(step_method, protocol=-1)
self._samplers = [
--> 403     ProcessAdapter(
       draws,
       tune,
       step_method,
       step_method_pickled,
       chain,
       seed,
       start,
       mp_ctx,
   )
   for chain, seed, start in zip(range(chains), seeds, start_points)
]
self._inactive = self._samplers.copy()
self._finished: List[ProcessAdapter] = []

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/parallel.py:259, in ProcessAdapter.__init__(self, draws, tune, step_method, step_method_pickled, chain, seed, start, mp_ctx)
   step_method_send = step_method
self._process = mp_ctx.Process(
   daemon=True,
   name=process_name,
   (...)
   ),
)
--> 259 self._process.start()
# Close the remote pipe, so that we get notified if the other
# end is closed.
remote_conn.close()

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/multiprocessing/process.py:121, in BaseProcess.start(self)
assert not _current_process._config.get('daemon'), \
      'daemonic processes are not allowed to have children'
_cleanup()
--> 121 self._popen = self._Popen(self)
self._sentinel = self._popen.sentinel
# Avoid a refcycle if the target function holds an indirect
# reference to the process object (see bpo-30775)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/multiprocessing/context.py:300, in ForkServerProcess._Popen(process_obj)
@staticmethod
def _Popen(process_obj):
   from .popen_forkserver import Popen
--> 300     return Popen(process_obj)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/multiprocessing/popen_forkserver.py:35, in Popen.__init__(self, process_obj)
def __init__(self, process_obj):
   self._fds = []
---> 35     super().__init__(process_obj)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/multiprocessing/popen_fork.py:19, in Popen.__init__(self, process_obj)
self.returncode = None
self.finalizer = None
---> 19 self._launch(process_obj)

File ~/opt/anaconda3/envs/pymc_env/lib/python3.11/multiprocessing/popen_forkserver.py:58, in Popen._launch(self, process_obj)
self.finalizer = util.Finalize(self, util.close_fds,
                              (_parent_w, self.sentinel))
with open(w, 'wb', closefd=True) as f:
---> 58     f.write(buf.getbuffer())
self.pid = forkserver.read_signed(self.sentinel)

KeyboardInterrupt: 

Esaminiamo ora l’adeguatezza del modello confrontando le predizioni del modello con i dati osservati.

with model:
    idata.extend(pm.sample_posterior_predictive(idata, random_seed=RANDOM_SEED))

Sampling: [obs]

100.00% [4000/4000 00:00<00:00]

fig, ax = plt.subplots(figsize=figsize)

az.plot_hdi(pre.index, idata.posterior_predictive["obs"], hdi_prob=0.5, smooth=False)
az.plot_hdi(pre.index, idata.posterior_predictive["obs"], hdi_prob=0.95, smooth=False)
ax.plot(pre.index, pre["deaths"], label="observed")
format_x_axis(ax)
ax.set(title="Posterior predictive distribution in the pre COVID-19 era")
plt.legend();

../_images/99e035f7fd15865a11cae7452a317cd1ea951eaf7e1746596cdc5605c8bf2f10.png

Il modello, pur non essendo esente da limitazioni, dimostra un’alta capacità predittiva del numero di decessi in relazione al tempo, avvalendosi unicamente del tempo e della temperatura come variabili indipendenti.

Ora focalizziamoci sui dati osservati, aggregati esclusivamente in base al fattore temporale e trascurando l’effetto della temperatura, per metterli a confronto con le stime fornite dal modello.

temp = idata.posterior["mu"].mean(dim=["chain", "draw"]).to_dataframe()
pre = pre.assign(deaths_predicted=temp["mu"].values)

fig, ax = plt.subplots(1, 2, figsize=figsize, sharey=True)
sns.lineplot(data=pre, x="year", y="deaths", ax=ax[0], lw=3)
ax[0].set(title="Observed")
sns.lineplot(data=pre, x="year", y="deaths_predicted", ax=ax[1], lw=3)
ax[1].set(title="Model predicted mean");

../_images/10d8082bb44e3d369a6ea01ade564d2ca4efef4cd8a09e4d9db6d92a149056e9.png

Possiamo notare una buona corrispondenza tra i dati osservati e le predizioni del modello.

Inferenza controfattuale#

Ora utilizzeremo il metodo sample_posterior_predictive per generare campioni dalla distribuzione predittiva a posteriori del modello, dato un nuovo set di dati sui predittori (“time” e “temp” in questo caso, per gli anni 2006-2020). In altre parole, stiamo cercando di capire quale sarebbe la distribuzione del numero di morti in anni futuri (2021-2023) dati i nuovi valori per il tempo e la temperatura, assumendo che il “meccanismo generatore dei dati” non cambi nel tempo – ovvero assumiamo che il Covid-19 non si sia verificato, insieme a tutte le altre cose che lo hanno accompagnato, come la riduzione dei servizi sanitari che, normalmente, sarebbero stati presenti.

Questo è un esempio di “estrapolazione” nel senso che stiamo utilizzando il modello addestrato su un certo set di dati per fare previsioni su un diverso set di dati che non è stato utilizzato durante l’addestramento del modello. Quindi, in termini semplici, stiamo cercando di prevedere come il numero di morti cambierà in funzione del tempo e della temperatura in anni futuri (2021-2023), sotto l’ipotesi che la relazione tra queste variabili e il numero di morti rimanga costante.

with model:
    pm.set_data(
        {
            "time": post["t"].to_numpy(),
            "temp": post["temp"].to_numpy(),
        }
    )
    counterfactual = pm.sample_posterior_predictive(
        idata, var_names=["obs"], random_seed=rng
    )

Sampling: [obs]

100.00% [4000/4000 00:00<00:00]

La figura illustra le previsioni sul numero di morti per gli anni 2021-2023, basate esclusivamente sui dati relativi al tempo e alla temperatura. Queste stime presuppongono che l’impatto di queste variabili sia analogo a quello osservato nel periodo 2006-2020.

fig, ax = plt.subplots(figsize=figsize)

plot_xY(post.index, counterfactual.posterior_predictive["obs"], ax)
format_x_axis(ax, minor=True)
ax.plot(post.index, post["deaths"], label="reported deaths")
ax.set(title="Counterfactual: Posterior predictive forecast of deaths if COVID-19 had not appeared")
plt.legend();

../_images/d74b4fb9407d04c124b2e59bf463fd60cab88c08c56b708baec6b0ec9f58eb4b.png

Alla figura è sovrapposta una spezzata (la linea blu) che, nel periodo post-covid, descrive il numero effettivo di morti.

Ora abbiamo gli elementi necessari per calcolare il numero di morti in eccesso. In particolare, il numero di morti segnalato e la previsione controfattuale Bayesiana su quanti sarebbero morti se nulla fosse cambiato dal periodo pre a quello post-COVID-19.

Inferenza controfattuale

Contents

Inferenza controfattuale#

Preparazione del Notebook#

Calcolo delle morti in eccesso dovute al COVID-19#

Inferenza controfattuale#

Numero di morti in eccesso dall’onset del Covid-19#