Modelli lineari generalizzati

Per chi vuole andare a fondo

La struttura di un GLM

Un modello lineare generalizzato ha tre componenti:

• una distribuzione della famiglia esponenziale per \(Y_i\) (normale, binomiale, Poisson, gamma…), con media \(\mu_i = E(Y_i)\);
• un predittore lineare \(\eta_i = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip}\);
• una funzione di legame \(g\) che collega le due: \(g(\mu_i) = \eta_i\).

La regressione lineare è il caso particolare con distribuzione normale e legame identità.

Regressione logistica

$$\log\frac{p_i}{1 - p_i} = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip}, \qquad p_i = \frac{1}{1 + e^{-\eta_i}}$$

\(e^{\beta_j}\) è l'odds ratio per un aumento unitario di \(x_j\). Il legame con il rischio relativo, per un rischio di base \(p_0\), è \(RR = OR / (1 - p_0 + p_0\cdot OR)\): i due coincidono solo quando \(p_0\) è piccolo.

Poisson con offset

Per un conteggio \(Y_i\) osservato su una popolazione (o un tempo) \(N_i\):

$$\log E(Y_i) = \log N_i + \beta_0 + \beta_1 x_{i1} + \dots \quad\Longleftrightarrow\quad \log\frac{E(Y_i)}{N_i} = \beta_0 + \beta_1 x_{i1} + \dots$$

Il termine \(\log N_i\) entra con coefficiente fissato a 1 (l'offset), e \(e^{\beta_j}\) è un rapporto tra tassi (incidence rate ratio).

Sovradispersione

La Poisson impone \(\mathrm{Var}(Y) = \mu\). Un indice pratico è il rapporto tra la statistica di Pearson e i gradi di libertà residui:

$$\hat\phi = \frac{1}{n - p - 1}\sum_i \frac{(y_i - \hat\mu_i)^2}{\hat\mu_i}$$

Valori chiaramente sopra 1 indicano sovradispersione. Le soluzioni sono il modello quasi-Poisson, che moltiplica gli errori standard per \(\sqrt{\hat\phi}\), o la binomiale negativa, con \(\mathrm{Var}(Y) = \mu + \mu^2/\theta\).

Stima

I coefficienti si stimano per massima verosimiglianza, con un algoritmo iterativo (IRLS). Il confronto tra modelli annidati si fa con il test del rapporto di verosimiglianza, basato sulla differenza delle devianze.

In R

Logistica e Poisson in R

set.seed(111)
n <- 2000
comorb <- rpois(n, 2)
eta    <- round(runif(n, 50, 90))
p      <- plogis(-3 + 0.45 * comorb + 0.01 * (eta - 70))
riamm  <- rbinom(n, 1, p)

logit <- glm(riamm ~ comorb + eta, family = binomial)
round(exp(cbind(OR = coef(logit), confint.default(logit))), 3)

# Tassi di ospedalizzazione per distretto: Poisson con offset
distr <- data.frame(pop = c(42000, 58000, 31000, 77000, 25000, 64000),
                    anziani = c(0.22, 0.25, 0.28, 0.20, 0.31, 0.24))
set.seed(5)
distr$ricoveri <- rnbinom(6, size = 30,
  mu = distr$pop * 0.004 * exp(3 * (distr$anziani - 0.25)))
pois <- glm(ricoveri ~ anziani + offset(log(pop)),
            family = poisson, data = distr)
sum(residuals(pois, type = "pearson")^2) / df.residual(pois)

Come leggere l'output

Ogni comorbidità in più moltiplica l'odds di riammissione per 1,42 (IC 95%: 1,30–1,55), a parità di età. L'età, con un effetto simulato molto piccolo, ha un odds ratio di 1,007 per anno e un intervallo che include 1.

Nel modello di Poisson per i distretti l'indice di dispersione vale 2,5: i conteggi variano due volte e mezzo più di quanto la Poisson ammetta, perché sono stati generati da una binomiale negativa. In questo caso conviene passare a glm.nb() del pacchetto MASS, oppure a family = quasipoisson.