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
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
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)
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.
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.
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
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.
Procediamo con l'algoritmo di discesa del gradiente.
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()
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¶
- Discesa del Gradiente
- Gradient Descent
- Offset nella Regressione di Poisson
- Fitting GLM by hand
- The Exposure Variable in Poisson Regression Models
- GLM in statsmodels (python)
- I modelli GLM in python
- Modello lineare generalizzato
- INTRODUCTION TO GENERALIZED LINEAR MIXED MODELS
- HOW CAN I ESTIMATE RELATIVE RISK USING GLM FOR COMMON OUTCOMES IN COHORT STUDIES? | STATA FAQ
- Poisson Regression
© 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