COVID-19 Data Analysis

GLM Poisson

Il Modello Lineare Generalizzato con distribuzione di Poisson e offset per la stima di efficacia.

Max Pierini


Negli aggiornamenti dell'Istituto Superiore di Sanità sulla pandemia di COVID-19 in Italia, e relativa campagna vaccinale, viene riportata una stima dell'efficacia vaccinale per stato vaccinale e classi d'età.

Nelle note si legge che è stato utilizzato un Modello Lineare Generalizzato con distribuzione di Poisson e offset. Come funziona? Ne diamo qui una breve spiegazione riassuntiva con qualche esempio.

Cenni di teoria

Il Modello Lineare Generalizzato (GLM) con distribuzione di Poisson è uno tra i più utilizzati modelli di regressione ed è definito come

\begin{align}\label{first} &\log(\lambda) = \sum_{n=0}^{N} \beta_n X_n + \varepsilon \tag{1.1} \\ &y \sim \mathbf{Poisson}(\lambda) \tag{1.2} \end{align}

dove $y$ è la variabile dipendente osservata distribuita Poisson con parametro $\lambda$ (vedremo tra poco cosa significa), $X_n$ sono le $N$ variabili indipendenti osservate (dette anche covariate) e $\beta_n$ sono i coefficienti di regressione. Il termine $\varepsilon$ rappresenta l'errore che si commette nella stima e possiamo, per l'attuale trattazione, non considerarlo.

Per $n=0$ la variabile $X_0$ è posta per definizione uguale a $1$ se si intende che il modello abbia un'intercetta oppure uguale a $0$ se non si vuole che il modello abbia un'intercetta.

Se, ad esempio, avessimo una sola variabile indipendente il modello con intercetta risulterebbe

\begin{align} &\log(\lambda) = \beta_0 + \beta_1 X_1 \tag{2.1} \\ &y \sim \mathbf{Poisson}(\lambda) \tag{2.2} \end{align}

dove $\beta_0$ è l'intercetta, mentre il modello senza intercetta risulterebbe

\begin{align} &\log(\lambda) = \beta_1 X_1 \tag{3.1} \\ &y \sim \mathbf{Poisson}(\lambda) \tag{3.2} \end{align}

Dalla (\ref{first}), esponenziando termine a termine ed omettendo il termine errore $\varepsilon$, otteniamo che

\begin{equation} \lambda = \exp\left\{\sum_{n=0}^{N} \beta_n X_n\right\} \tag{4} \end{equation}

Come dicevamo, la variabile $y$ è distribuita Poisson, ovvero la sua Funzione di Massa di Probabilità (PMF) è definita come

\begin{equation} p(y;\lambda) = \frac{\lambda^y e^{-\lambda}}{y!} \tag{5} \end{equation}

e restituisce la probabilità che la variabile assuma un valore $y$ dato il parametro $\lambda$ ed ha media e varianza identiche

\begin{equation} \mu = \sigma^2 = \lambda \tag{6} \end{equation}

Avendo $M$ osservazioni delle variabili indipendenti $X_n$ e della variabile dipendente $y$, il modello risulta

\begin{align} &\log(\lambda_m) = \sum_{n=0}^{N} \beta_n X_{n,m} \tag{7.1} \\ &y_{m} \sim \mathbf{Poisson}(\lambda_{m}) \tag{7.2} \\ &p(y_m;\lambda_m) = \frac{\lambda_m^{y_m} e^{-\lambda_m}}{y_m!} \tag{7.3} \\ &p(y_1 ... y_m;\lambda_1 ... \lambda_m) = \prod_{m=1}^{M} \frac{\lambda_m^{y_m} e^{-\lambda_m}}{y_m!} \tag{7.4} \end{align}

ma è conveniente utilizzare la notazione matriciale

\begin{align}\label{mat} &\log(\boldsymbol\Lambda) = \mathbf{X} \boldsymbol\beta \tag{8.1} \\ &\mathbf{y} \sim \mathbf{Poisson}(\boldsymbol\Lambda) \tag{8.2} \end{align}

dove $\mathbf{y}$ è un vettore colonna di $M$ righe contenente tutte le $M$ osservazioni della variabile dipendente $y_m$ con $m$ che va da $1$ a $M$

\begin{equation} \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \cdots \\ y_M \end{bmatrix} \tag{9} \end{equation}

$\boldsymbol\Lambda$ è un vettore colonna di $M$ righe di tutti gli $M$ parametri $\lambda_m$ della distribuzione di Poisson con $m$ che va da $1$ a $M$

\begin{equation} \boldsymbol\Lambda = \begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_M \end{bmatrix} \tag{10} \end{equation}

$\boldsymbol\beta$ è un vettore colonna di $N$ righe di tutti gli $N$ coefficienti di regressione $\beta_n$ con $n$ che va da $1$ a $N$

\begin{equation} \boldsymbol\beta = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_N \end{bmatrix} \tag{11} \end{equation}

e $\mathbf{X}$ è una matrice di $M$ righe ed $N$ colonnne di tutte le $M$ osservazioni delle $N$ variabili indipendenti $X_{m,n}$ con $m$ che va da $1$ a $M$ e $n$ che va da $1$ a $N$

\begin{equation} \mathbf{X} = \begin{bmatrix} \begin{array}{c|cccc} & \textbf{var 0} & \textbf{var 1} & \cdots & \textbf{var N} \\ \hline \textbf{obs 1} & X_{1,0} & X_{1,1} & \cdots & X_{1,N} \\ \textbf{obs 2} & X_{2,0} & X_{2,1} & \cdots & X_{2,N} \\ \vdots & & & \vdots \\ \textbf{obs M} & X_{M,0} & X_{M,1} & \cdots & X_{M,N} \end{array} \end{bmatrix} \tag{12} \end{equation}

ovvero

\begin{equation} \log \begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_M \end{bmatrix} = \begin{bmatrix} X_{1,0} & X_{1,1} & \cdots & X_{1,N} \\ X_{2,0} & X_{2,1} & \cdots & X_{2,N} \\ \vdots & & & \vdots \\ X_{M,0} & X_{M,1} & \cdots & X_{M,N} \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_N \end{bmatrix} \tag{13.1} \end{equation}\begin{equation} \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_M \end{bmatrix} \sim \mathbf{Poisson} \begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_M \end{bmatrix} \tag{13.2} \end{equation}

Per trovare i coefficienti $\boldsymbol\beta_n$ dobbiamo ottimizzare la log-likelihood negativa.

La likelihood $\mathcal{L}(\boldsymbol\Lambda)$ è proporzionale alla PMF $p(\mathbf{y};\boldsymbol\Lambda)$ a meno del denominatore $y!$ che non dipende da $\boldsymbol\beta$ e si può quindi trattare come una costante

\begin{align} \mathcal{L}(\boldsymbol\Lambda) = \prod_{m=1}^{M} \lambda_m^{y_m} e^{-\lambda_m} \tag{14.1} \\ \mathcal{L}(\boldsymbol\Lambda) \propto p(y_1 ... y_m;\lambda_1 ... \lambda_m) \tag{14.2} \end{align}

La log-likelihood $\mathscr{L}(\boldsymbol\Lambda)$ sarà pertanto

\begin{align} \mathscr{L}(\boldsymbol\Lambda) &= \log[ \mathcal{L}(\boldsymbol\Lambda) ] \tag{15.1} \\&= \log\left[ \prod_{m=1}^{M} \lambda_m^{y_m} e^{-\lambda_m} \right] \tag{15.2} \\&= \sum_{m=1}^M \log\left[ \lambda_m^{y_m} e^{-\lambda_m} \right] \tag{15.3} \\&= \sum_{m=1}^M \log\Big( \exp\big\{ y_m \log(\lambda_m) \big\} \Big) + \log\left( e^{-\lambda_m} \right) \tag{15.4} \\&= \sum_{m=1}^M y_m \log(\lambda_m) - \lambda_m \tag{15.5} \end{align}

quindi la log-likelihood negativa

\begin{equation}\label{ell} \ell(\boldsymbol\Lambda) = -\mathscr{L}(\boldsymbol\Lambda) = \sum_{m=1}^M \lambda_m - y_m \log(\lambda_m) \tag{16} \end{equation}

Esponenziando termine a termine la (\ref{mat}) otteniamo

\begin{equation} \boldsymbol\Lambda = e^{\mathbf{X} \boldsymbol\beta } \tag{17} \end{equation}

da cui, sostituendo $\boldsymbol\Lambda$ nella (\ref{ell}) otteniamo la log-likelihood negativa in funzione di $\boldsymbol\beta$

\begin{align} \ell(\boldsymbol\beta, \mathbf{X}) &= \sum_{m=1}^M e^{\mathbf{X} \boldsymbol\beta} - y_m \mathbf{X} \boldsymbol\beta \tag{18} \end{align}

e possiamo ora calcolarne la derivata parziale prima in $\partial \boldsymbol\beta$ ovvero il gradiente rispetto a tutte le componenti $\beta_n$

\begin{align} \frac{\partial \ell}{\partial \boldsymbol\beta} &= \sum_{m=1}^M \mathbf{X}^T e^{\mathbf{X} \boldsymbol\beta} - y_m \mathbf{X}^T \tag{19.1} \\&= \sum_{m=1}^M \mathbf{X}^T \left( e^{ \mathbf{X} \boldsymbol\beta} - y_m\right) \tag{19.2} \\&= \mathbf{X}^T \left( \hat{\boldsymbol\Lambda} - \mathbf{y}\right) \tag{19.3} \\&= \nabla\ell \tag{19.4} \end{align}

dove $\nabla\ell$ è il gradiente calcolato su ciasun coefficiente $\beta_n$ ed è perciò un vettore colonna di $N$ righe.

\begin{equation} \begin{bmatrix} X_{1,0} & \cdots & X_{1,N} \\ X_{2,0} & \cdots & X_{2,N} \\ \vdots & & \vdots \\ X_{M,0} & \cdots & X_{M,N} \end{bmatrix}^T \left( \begin{bmatrix} \hat{\lambda}_1 \\ \hat{\lambda}_2 \\ \vdots \\ \hat{\lambda}_M \end{bmatrix} - \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_M \end{bmatrix} \right) = \begin{bmatrix} \frac{\partial \ell}{\partial \beta_0} \\ \frac{\partial \ell}{\partial \beta_1} \\ \vdots \\ \frac{\partial \ell}{\partial \beta_N} \end{bmatrix} \tag{20} \end{equation}

Ponendo il gradiente uguale a zero otteniamo i punti critici ottimizzando $\boldsymbol\beta$. Questo può essere ottenuto tramite un algoritmo di discesa del gradiente il cui pseudo-codice più semplice può essere definito come

\begin{align*} &\textrm{INIT }\hat{\boldsymbol\beta} \tag{21} \\ &\textrm{INIT }\hat{\boldsymbol\Lambda} = e^{\mathbf{X} \hat{\boldsymbol\beta}} \\ &\textrm{while }\nabla\ell \neq \mathbf{0} \\ &\;\;\;\;\;\; \nabla\ell = \mathbf{X}^T \left( \hat{\boldsymbol\Lambda} - \mathbf{y}\right) \\ &\;\;\;\;\;\;\hat{\boldsymbol\beta} = \hat{\boldsymbol\beta} - \eta \nabla\ell \\ &\;\;\;\;\;\;\hat{\boldsymbol\Lambda} = e^{\mathbf{X} \hat{\boldsymbol\beta}} \end{align*}

dove $\hat{\boldsymbol\beta}$ può essere inizializzato a valori random oppure vicini a zero ed $\eta$ è il tasso di apprendimento (learning rate).

Più precisamente, l'ottimizzazione dei coefficienti $\beta_n$ con algortimo di discesa del gradiente permette di trovare i valori dei coefficienti tali per cui i punti $(X_{m,n}, y_m)$ osservati si trovino "il più possibile" sul grafico della funzione della log-likelihood negativa $\ell(\boldsymbol\Lambda)$.


Nel caso speciale in cui ci sia una sola variabile indipendente senza intercetta sappiamo che

\begin{align} &\log(\lambda) = \beta X \tag{22.1} \\ &y \sim \mathbf{Poisson}(\lambda) \tag{22.2} \\ &p(y; \lambda) = \frac{\lambda^{y} e^{-\lambda}}{y!} \tag{22.3} \end{align}

da cui

\begin{equation} \lambda = e^{\beta X} \tag{23} \end{equation}

e la Likelihood per ciasun $\lambda$ sarà

\begin{equation} \mathcal{L}(\lambda) = \lambda^{y} e^{-\lambda} \tag{24} \end{equation}

pertanto la log-likelihood negativa corrisponde a

\begin{equation} \ell(\lambda) = e^{\beta X} - \beta X y \tag{25} \end{equation}

Il gradiente rispetto a $\beta$ si riduce quindi alla sola derivata parziale prima in $\partial \beta$

\begin{equation} \nabla\ell = \frac{\partial \ell}{\partial \beta} = X e^{\beta X} - X y \tag{26} \end{equation}

Ponendolo uguale a zero e risolvendo per $\beta$

\begin{align} &X e^{\beta X} - X y = 0 \tag{27.1} \\ &X e^{\beta X} = X y \tag{27.2} \\ &e^{\beta X} = y \tag{27.3} \\ &\beta X = \log(y) \tag{27.4} \end{align}

che è esattamente la definizione di GLM con distribuzione di Poisson.

Dunque la log-likelihood negativa $\ell(\lambda)$ ha un punto critico per $\beta = \frac{\log(y)}{X}$; verifichiamo facilmente che $\nabla\ell > 0$ per $\beta > \frac{\log(y)}{X}$ ed è quindi un punto di minimo della log-likelihood negativa $\ell(\lambda)$ rispetto a $\beta$.

Con più parametri $\beta_n$ non è possibile porre "direttamente" $\nabla\ell=0$ e dobbiamo appunto procedere con un algoritmo, come quello di discesa del gradiente.

I questo esempio GeoGebra

  • $\mathrm{X}$ e $\mathrm{Y}$ sono le variabili osservate $X_m$ e $y_m$
  • $\beta$ è il coefficiente reale
  • $\beta_{\mathrm{hat}}$ è il coefficiente stimato $\hat{\beta}$
  • $f(x,y)$ è la log-likelihood negativa $\ell(\lambda)$
  • $\mathrm{I1}$ sono i punti osservati ($X_m$, $y_m$) e la la posizione sul grafico di $\ell(\lambda)$

Modificando $\beta_{\mathrm{hat}}$ si può vedere come la curva si "adatti" ai dati osservati per $\beta_{\mathrm{hat}}\rightarrow\beta$.

In questo esempio GeoGebra si possono visualizzare distribuzione di Poisson $p(y;\lambda)$, Likelihoood $\mathcal{L}$, log-likelihood $\mathscr{L}$, log-likelihood negativa $\ell$ e componenti del gradiente $\nabla\ell$ di un GLM con una variabile $X$ e intercetta $\beta_0$


Analizzandolo da un ulteriore punto di vista, nel caso in cui si contempli una variabile $X$ e un'intercetta, date $M$ osservazioni, sappiamo che la log-likelihood negativa per ciascuna delle osservazioni è

\begin{equation} \tag{L1} \ell_m = e^{\beta_0 + \beta_1 X_m} - y_m (\beta_0 + \beta_1 X_m) \end{equation}

dove $X_m$ e $y_m$ sono le $M$ osservazioni della variabile indipendente e dipendente.

Il gradiente $\nabla\ell_m$ rispetto ai due coefficienti $\beta_n$ è quindi il vettore delle derivate parziali prime rispetto a ciascun coefficiente, ovvero

\begin{align} \nabla\ell_m &= \left[ \frac{\partial\ell_m}{\partial\beta_0}, \frac{\partial\ell_m}{\partial\beta_1} \right]^T \tag{L2.1} \\&= \begin{bmatrix} e^{\beta_0 + \beta_1 X_m} - y_m \\ X_m ( e^{\beta_0 + \beta_1 X_m} - y_m) \end{bmatrix} \tag{L2.2} \end{align}

Ponendo entrambe le componenti uguali a zero e risolvendo, ad esempio, per $\beta_1$ in funzione di $\beta_0$ otteniamo

\begin{align} & \nabla\ell_m = \mathbf{0} \;\;\;\;\;\; \textrm{per} \tag{L3.1} \\& e^{\beta_0 + \beta_1 X_m} - y_m = 0 \tag{L3.2} \\& \beta_1 = \frac{\log(y_m) - \beta_0}{X_m} \tag{L3.3} \\& \beta_1 = \frac{\log(y_m)}{X_m} -\frac{1}{X_m} \beta_0 \tag{L3.4} \end{align}

ovvero un fascio di $M$ rette di pendenza $\frac{-1}{X_m}$ e intercetta $\frac{\log(y_m)}{X_m}$ che rappresentano il minimo della log-likelihood negativa $\ell_m$ per ciascuna delle $M$ osservazioni e che si intersecano nel punto $(\beta_0, \beta_1)$ ovvero all'interesezione dei minimi tutte le $\ell_m$.

In questo esempio GeoGebra

  • $\mathrm{X}$ sono le $X_m$ osservazioni
  • $\mathrm{Y}$ sono le $y_m$ osservazioni
  • sull'asse $x$ è rappresentato $\beta_0$
  • sull'asse $y$ è rappresentato $\beta_1$
  • il fascio di rette $\mathrm{f}$ è il minimo di ciascuna log-likelihood negativa $\ell_m$
  • l'intersezione del fascio di rette è dove $\nabla\ell=\mathbf{0}$

In questo esempio GeoGebra, si possono vedere le log-likelihood negative $\ell_m$ in tre dimensioni a, b, c (sempre in funzione di $\beta_0$ e $\beta_1$) e il punto dove i minimi di tutte le $\ell_m$ si intersecano ovvero dove $\nabla\ell_m=\mathbf{0} \;,\; \forall m \in (1 ... M)$. La superficie g è la somma delle $\ell_m$ ovvero la log-likelihood negativa di tutte le osservazioni.

Ovviamente, in tutti questi esempi grafici le variabili $y_m$ e $X_m$ sono definite per rispondere perfettamente alla definizione del Modello Lineare Generalizzato ovvero $\mathbf{y} = e^{\mathbf{X}\boldsymbol\beta}$.

Nella realtà, i dati raccolti non saranno così perfettamente aderenti alla definizione ma avremo per ciascuna coppia delle $M$ variabili osservate $y_m = \exp\left\{ \sum \beta_n X_{m,n} \right\} + \varepsilon_m$ dove $\varepsilon_m$ è un errore rispetto alla definizione, che rende conto del termine $\varepsilon$ definito nella (\ref{first}). Ecco perché si distribuisce Poisson la variabile osservata $\mathbf{y}$.

Dunque la ricerca dei parametri $\boldsymbol\beta$ tali per cui $\nabla\ell=\mathbf{0}$ è in realtà la ricerca di quei valori di $\boldsymbol\beta$ tali per cui $\nabla\ell\simeq\mathbf{0}$ o meglio dove $|\nabla\ell|<k$ con $0<k<1$.


Per la stima di efficacia di un vaccino, la variabile dipendente $y$ rappresenta il numero di eventi osservati (diagnosi, ospedalizzazioni o decessi a seconda del focus) e le variabili indipendenti $X$ sono lo stato vaccinale (che nel più semplice dei casi si riduce a non vaccinato e vaccinato) ed eventuali covariate come fascia d'età, sesso e qualunque altra variabile si ritenga opportuno considerare.

Ma questo non è sufficiente perché la stima di efficacia di un vaccino (così come di qualsiasi altro trattamento) è definita come Riduzione del Rischio Relativo ($\mathrm{RRR}$) a sua volta definito come il complementare del Rischio Relativo ($\mathrm{RR}$) che è il rapporto tra l'incidenza $\mathrm{EER}$ nei vaccinati (gruppo sperimentale) e l'incidenza $\mathrm{CER}$ nei non vaccinati (gruppo di controllo)

\begin{align} \mathrm{RRR} = 1 - \mathrm{RR} = 1 - \frac{\mathrm{EER}}{\mathrm{CER}} \tag{28} \end{align}

dove

\begin{align} \mathrm{EER} = \frac{\mathrm{EE}}{\mathrm{EN}} \tag{29.1} \\ \mathrm{CER} = \frac{\mathrm{CE}}{\mathrm{CN}} \tag{29.2} \end{align}

in cui chiamando $\mathbf{X_1}$ lo stato vaccinale e definendolo pari a 0 se non vaccinato e ad 1 se vaccinato

\begin{align} \mathrm{EE} = \mathbf{y}[\mathbf{X_1}=1] \tag{30.1} \\ \mathrm{CE} = \mathbf{y}[\mathbf{X_1}=0] \tag{30.2} \end{align}

ovvero sono rispettivamente gli eventi osservati nei vaccinati e nei non vaccinati.

È chiaro che ci manca da definire $\mathrm{EN}$ e $\mathrm{CN}$ ovvero i denominatori delle incidenze osservate, il numero di vaccinati e di non vaccinati osservati.

Questo può essere ottenuto in un GLM con distribuzione di Poisson aggiungendo un termine $\mathbf{z}$ denominato offset privo di coefficiente di regressione, prendendone il logaritmo (esposizione), ovvero

\begin{align}\label{zeta} &\log(\boldsymbol\Lambda) = \log(\mathbf{z}) + \mathbf{X} \boldsymbol\beta \tag{31.1} \\ &\mathbf{y} \sim \mathbf{Poisson}(\boldsymbol\Lambda) \tag{31.2} \end{align}

dove $\mathbf{z}$ è un vettore colonna di $M$ righe contenente tutti i denominatori delle incidenze (ovvero il numero di esposti) delle $M$ osservazioni.

Infatti, con semplici passaggi algebrici ed esponenziando termine a termine nella (\ref{zeta}) otteniamo

\begin{align} \log(\boldsymbol\Lambda) - \log(\mathbf{z}) = \mathbf{X} \boldsymbol\beta \tag{32.1} \\ \log\left( \boldsymbol\Lambda \oslash \mathbf{z} \right) = \mathbf{X} \boldsymbol\beta \tag{32.2} \\ \boldsymbol\Lambda \oslash \mathbf{z} = e^{ \mathbf{X} \boldsymbol\beta} \tag{32.3} \end{align}

dove $\boldsymbol\Lambda \oslash \mathbf{z}$ indica la divisione di Hadamard termine a termine (element-wise)

\begin{equation} \begin{bmatrix} \lambda_1 \\ \vdots \\ \lambda_m \end{bmatrix} \oslash \begin{bmatrix} z_1 \\ \vdots \\ z_m \end{bmatrix} = \begin{bmatrix} \frac{\lambda_1}{z_1} \\ \vdots \\ \frac{\lambda_m}{z_m} \end{bmatrix} \tag{32.3a} \end{equation}

è esattamente l'incidenza osservata e il coefficiente $\hat{\beta}_1$ stimato (se consideriamo $X_1$ come la variabile principale dello stato vaccinale) è il logaritmo del Rischio Relativo, da cui possiamo facilmente ottenere l'efficacia come Riduzione del Rischio Relativo

\begin{align} \mathrm{RRR} &= 1 - \mathrm{RR} \tag{33.1} \\&= 1 - e^{\hat{\beta}_1} \tag{33.2} \end{align}

Lo pseudo-codice dell'algoritmo di discesa del gradiente sarà in questo caso

\begin{align*} &\textrm{INIT }\hat{\boldsymbol\beta} \tag{34} \\ &\textrm{INIT }\hat{\boldsymbol\Lambda} = e^{\mathbf{X}\hat{\boldsymbol\beta}} \\ &\textrm{while }\nabla\ell \neq \mathbf{0} \\ &\;\;\;\;\;\;\nabla\ell = \mathbf{X}^T \left( \mathbf{z} \circ \hat{\boldsymbol\Lambda} - \mathbf{y}\right) \\ &\;\;\;\;\;\;\hat{\boldsymbol\beta} = \hat{\boldsymbol\beta} - \eta \nabla\ell \\ &\;\;\;\;\;\;\hat{\boldsymbol\Lambda} = e^{\mathbf{X}\hat{\boldsymbol\beta}} \end{align*}

dove $\mathbf{z} \circ \hat{\boldsymbol\Lambda}$ indica il prodotto di Hadamard termine a termine (element-wise)

\begin{equation} \begin{bmatrix} z_1 \\ \vdots \\ z_m \end{bmatrix} \circ \begin{bmatrix} \lambda_1 \\ \vdots \\ \lambda_m \end{bmatrix} = \begin{bmatrix} z_1 \lambda_1 \\ \vdots \\ z_m \lambda_m \end{bmatrix} \tag{34a} \end{equation}

Esempi

NB: i valori negli esempi sono volutamente randomizzati e saranno ogni giorno lievemente differenti

A puro scopo esemplificativo, ci limiteremo inizialmente al caso in cui $N=1$ ovvero abbiamo una sola variabile indipendente $X_{m,1}$ (corrispondente allo stato vaccinale: 0 per non vaccinato o 1 per vaccinato) e supporremo di avere 10 osservazioni, dunque $M=10$. Contempliamo l'intercetta, dunque avremo $X_{m,0}=1$ per ogni osservazione.

\begin{equation} \log \begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_9 \\ \lambda_{10} \end{bmatrix} = \log \begin{bmatrix} z_1 \\ z_2 \\ \vdots \\ z_9 \\ z_{10} \end{bmatrix} + \begin{bmatrix} 1 & X_{1,1} \\ 1 & X_{2,1} \\ \vdots & \vdots \\ 1 & X_{9,1} \\ 1 & X_{10,1} \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix} \tag{35.1} \end{equation}\begin{equation} \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_9 \\ y_{10} \end{bmatrix} \sim \mathbf{Poisson} \begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_9 \\ \lambda_{10} \end{bmatrix} \tag{35.2} \end{equation}

Il gradiente della log-likelihood negativa $\ell(\boldsymbol\Lambda)$ risulterà

\begin{equation} \nabla\ell = \begin{bmatrix} 1 & X_{1,1} \\ 1 & X_{2,1} \\ \vdots & \vdots \\ 1 & X_{9,1} \\ 1 & X_{10,1} \end{bmatrix}^T \left( \begin{bmatrix} z_1 \lambda_1 \\ z_2 \lambda_2 \\ \vdots \\ z_9 \lambda_9 \\ z_{10} \lambda_{10} \end{bmatrix} - \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_9 \\ y_{10} \end{bmatrix} \right) = \begin{bmatrix} \frac{\partial \ell}{\partial \beta_0} \\ \frac{\partial \ell}{\partial \beta_1} \end{bmatrix} \tag{36} \end{equation}

Lo pseudo-codice dell'algoritmo di discesa del gradiente sarà in questo caso

\begin{align*} &\textrm{INIT } \hat{\beta}_0, \hat{\beta}_1 \tag{37} \\ &\textrm{INIT }\hat{\boldsymbol\Lambda} = e^{\hat{\beta}_0 + \hat{\beta}_1 \mathbf{X_1}} \\ &\textrm{while }\nabla\ell \neq \mathbf{0} \\ &\;\;\;\;\;\;\nabla\ell = \mathbf{X}^T \left( \mathbf{z} \circ \hat{\boldsymbol\Lambda} - \mathbf{y}\right) \\ &\;\;\;\;\;\;\hat{\boldsymbol\beta} = \hat{\boldsymbol\beta} - \eta \nabla\ell \\ &\;\;\;\;\;\;\hat{\boldsymbol\Lambda} = e^{\hat{\beta}_0 + \hat{\beta}_1\mathbf{X_1}} \end{align*}

Scegliamo di fermare l'algoritmo dopo almeno 500 iterazioni di burn-in e quando le derivate parziali prime della log-likelihood negativa siano almeno inferiori a 0.01 in valore assoluto, con un tasso di apprendimento $\eta=10^{-8}$. I parametri $\beta_n$ saranno inizializzati a 0.001.

vaccinated exposed events incidence
0 4460972 609441 13.66%
0 4097820 547765 13.37%
0 4854459 661916 13.64%
0 4961225 670828 13.52%
0 4638226 627889 13.54%
1 16612437 502341 3.02%
1 16051720 493331 3.07%
1 19788274 604946 3.06%
1 18272568 599920 3.28%
1 16887232 532569 3.15%

Sommando semplicemente tutti gli eventi e tutti gli esposti per stato vaccinale e calcolando la Riduzione del Rischio Relativo (RRR) otteniamo

$$ \mathrm{RRR} = 1 - \frac{\mathrm{EER}}{\mathrm{CER}} = 1 - \frac{ 3.12\% }{ 13.55\% } = 76.97\% $$

Proviamo ora ad applicare l'algoritmo di discesa del gradiente.

Un'implementazione semplice in python può essere ad esempio

## DISCESA DEL GRADIENTE ##
learning_rate = 1e-8
grad_limit = 0.01
# initialize beta_hats and gradient
beta_hat = np.array([0.001,0.001])
grad = np.array([np.nan, np.nan])
# initialize lam_hats
lam_hat = np.exp(X @ beta_hat)

run = 0
while True:
    cond1 = abs(grad[0]) < grad_limit
    cond2 = abs(grad[1]) < grad_limit
    cond3 = run > 500
    if cond1 and cond2 and cond3:
        # stop iterations
        break
    # get gradient 
    c = Z*lam_hat - y
    grad = X.T @ c
    # adjust beta hats
    beta_hat -= learning_rate * grad
    # adjust lambda hats
    theta = X @ beta_hat
    #theta[theta==0] = 1e-5
    lam_hat = np.exp(theta)

    run += 1

Gradiente finale $\nabla\ell=$ [-0.00579024 0.00996857]

Coefficienti beta:

  • $\hat{\beta}_0=-1.99891$

  • $\hat{\beta}_1=-1.46858$

Verifichiamo la stima di efficacia ottenuta con il Modello Lineare Generalizzato

$$ E = 1 - e^{\beta_1} = 1 - e^{-1.46858} = 76.97\% $$

uguale all'efficacia stimata con la semplice formula del Rischio Relativo 76.97%.


Complichiamo ora il problema aggiungendo una covariata: la fascia d'età $X_2$. Supponiamo di avere solo 2 fasce d'età ovvero "0-50" e "50+" (che corrisponderanno a variabili dummies codificate come 0 e 1) e che l'efficacia sia superiore nella fascia più anziana ma l'incidenza degli eventi più elevata nelle fasce più giovani e la popolazione formata da individui più anziani.

vaccinated age exposed events incidence
0 0 1286970 261730 20.34%
0 0 1404676 312374 22.24%
0 0 1306455 267483 20.47%
0 0 1375076 281098 20.44%
0 0 1301730 307789 23.64%
0 1 4669733 638466 13.67%
0 1 4263495 589864 13.84%
0 1 4379163 580169 13.25%
0 1 4528152 623796 13.78%
0 1 4256593 576165 13.54%
1 0 16686966 2607480 15.63%
1 0 17657159 2883793 16.33%
1 0 19968294 3270225 16.38%
1 0 17256122 2844391 16.48%
1 0 15397840 2378999 15.45%
1 1 28981520 928988 3.21%
1 1 21696122 662640 3.05%
1 1 22931122 747580 3.26%
1 1 27944845 892531 3.19%
1 1 25475233 814555 3.20%

Sommando semplicemente tutti gli eventi e tutti gli esposti per stato vaccinale, senza tener conto delle fasce d'età, calcoliamo la Riduzione del Rischio Relativo (RRR)

$$ \mathrm{RRR} = 1 - \frac{\mathrm{EER}}{\mathrm{CER}} = 1 - \frac{ 8.43\% }{ 15.43\% } = 45.38\% $$

Procediamo con l'algoritmo di discesa del gradiente.

Gradiente finale $\nabla\ell=$ [-0.00997112 0.00983944 0.00404835]

Coefficienti beta stimati

  • $\hat{\beta}_0=-1.03269$

  • $\hat{\beta}_1=-0.86490$

  • $\hat{\beta}_2=-1.33890$

Verifichiamo la stima di efficacia ottenuta con il Modello Lineare Generalizzato

$$ E = 1 - e^{\beta_1} = 1 - e^{-0.86490} = 57.89\% $$

decisamente superiore all'efficacia stimata con la semplice formula del Rischio Relativo 45.38%.

Questo avviene a causa della differente distribuzione dell'incidenza nelle fasce d'età e per la composizione della popolazione. La stima fatta precedentemente con il semplice calcolo della Riduzione del Rischio Relativo non era in grado di tener conto di queste variabili. La stima con GLM invece assicura che il risultato sia a parità delle altre variabili.

Verifichiamo anche che i coefficienti $\beta_n$ ottenuti "manualmente" con il metodo della discesa del gradiente corrispondono a quelli ottenuti con statsmodels

mod = sm.GLM.from_formula(
    "events ~ vaccinated + age",
    exposure=df.exposed, 
    data=df, family=sm.families.Poisson()
)
fit = mod.fit()
fit.summary()
Generalized Linear Model Regression Results
Dep. Variable: events No. Observations: 20
Model: GLM Df Residuals: 17
Model Family: Poisson Df Model: 2
Link Function: log Scale: 1.0000
Method: IRLS Log-Likelihood: -5.4789e+05
Date: Sun, 22 May 2022 Deviance: 1.0955e+06
Time: 17:36:02 Pearson chi2: 1.07e+06
No. Iterations: 5
Covariance Type: nonrobust
coef std err z P>|z| [0.025 0.975]
Intercept -1.0327 0.001 -1983.503 0.000 -1.034 -1.032
vaccinated -0.8649 0.001 -1611.044 0.000 -0.866 -0.864
age -1.3389 0.000 -2907.087 0.000 -1.340 -1.338

Coefficienti beta stimati

  • $\hat{\beta}_0=-1.03269$

  • $\hat{\beta}_1=-0.86490$

  • $\hat{\beta}_2=-1.33890$

uguali a quelli stimati con l'algoritmo di discesa del gradiente implementato sopra.


Le covariate potrebbero essere molte, ad esempio: il sesso, l'incidenza regionale, la regione di somministrazione e in generale tutto ciò che si ritenga possa influenzare il risultato e di cui si dispongano sufficienti dati. All'aumentare delle coavariate, aumenta anche la dimensione dei vettori e in particolare del gradiente $\nabla\ell$ della log-likelihood negativa.

Il Modello Lineare Generalizzato assicura che l'efficacia complessiva della popolazione in oggetto sia stimata a parità delle altre covariate ed è quindi un potente strumento statistico il cui utilizzo è facilitato da software specifici, come statsmodels ed altri più avanzati e robusti.

Appronfondimenti


© 2020 Max Pierini. Thanks to Sandra Mazzoli & Alessio Pamovio. ipynb-website © 2017 Peter Carbonetto & Gao Wang

Exported from Italia/GLM_Poisson.ipynb committed by maxdevblock on Sun May 22 17:50:34 2022 revision 126, 97ae0d9e