Modello lineare generalizzato

I modelli lineari generalizzati (GLM) sono una generalizzazione del più classico modello lineare nell'ambito della regressione lineare. Mentre nel modello lineare classico si ipotizza che la variabile endogena sia distribuita in modo normale, nell'ambito dei modelli lineari generalizzati la variabile endogena può essere distribuita come una qualsiasi variabile casuale della famiglia esponenziale e dunque, oltre alla variabile casuale normale, entrano in gioco anche altre variabili casuali, quali la binomiale, la poissoniana, la gamma, la normale inversa e altre.

I modelli lineari generalizzati vennero formulati da John Nelder e Robert Wedderburn come un modo per uniformare all'interno di un unico modello diversi altri modelli statistici, compreso il modello lineare, le regressione logistica e la regressione poissoniana. Si riesce in questo modo a incorporare in un unico modello anche altri modelli oltre a quello lineare.^[1]

Intuizione

La regressione lineare ordinaria prevede il valore atteso di una data quantità sconosciuta (la variabile di risposta, una variabile casuale) come combinazione lineare di un insieme di valori osservati (predittori). Ciò implica che un cambiamento costante in un predittore porta a un cambiamento costante nella variabile di risposta (cioè un modello a risposta lineare). Ciò è appropriato quando la variabile di risposta può variare, con buona approssimazione, indefinitamente in entrambe le direzioni, o più generalmente per qualsiasi quantità che varia solo di una quantità relativamente piccola rispetto alla variazione delle variabili predittive, ad es. altezze umane.

Tuttavia, queste ipotesi non sono appropriate per alcuni tipi di variabili di risposta. Ad esempio, nei casi in cui si prevede che la variabile di risposta sia sempre positiva e che vari in un ampio intervallo, i cambiamenti costanti dell'input portano a variazioni di output che variano geometricamente (cioè in modo esponenziale), piuttosto che costantemente. Ad esempio, supponiamo che un modello di previsione lineare apprenda da alcuni dati (forse tratti principalmente da grandi spiagge) che una diminuzione della temperatura di 10 gradi porterebbe a 1.000 persone in meno a visitare la spiaggia. È improbabile che questo modello si generalizzi bene su spiagge di dimensioni diverse. Più specificamente, il problema è che se si utilizza il modello per prevedere la nuova presenza con un calo di temperatura di 10 gradi per una spiaggia che riceve regolarmente 50 bagnanti, si prevede un valore di presenza impossibile di -950 bagnanti. Logicamente, un modello più realistico prevedrebbe invece un tasso costante di maggiore frequentazione della spiaggia (ad es. un aumento di 10 gradi porta a un raddoppio della frequentazione della spiaggia e un calo di 10 gradi porta a un dimezzamento delle presenze). Tale modello è definito modello a risposta esponenziale (o modello log-lineare, poiché si prevede che il logaritmo della risposta vari linearmente).

Allo stesso modo, un modello che predice una probabilità di fare una scelta sì/no (una variabile di Bernoulli) è ancora meno adatto come modello a risposta lineare, poiché le probabilità sono limitate su entrambe le estremità (devono essere comprese tra 0 e 1). Si immagini, ad esempio, un modello che prevede la probabilità che una determinata persona vada in spiaggia in funzione della temperatura. Un modello ragionevole potrebbe prevedere, ad esempio, che una variazione di 10 gradi renda una persona due volte più o meno propensa ad andare in spiaggia. Ma cosa significa "due volte più probabile" in termini di probabilità? Non può letteralmente significare raddoppiare il valore di probabilità (ad es. 50% diventa 100%, 75% diventa 150%, ecc.). Piuttosto, sono le probabilità che raddoppiano: da 2:1 a 4:1, a 8:1, ecc. Tale modello è un modello logistico.

I modelli lineari generalizzati coprono tutte queste situazioni consentendo variabili di risposta che hanno distribuzioni arbitrarie (piuttosto che distribuzioni semplicemente normali) e che una funzione arbitraria della variabile di risposta (la funzione di collegamento) vari linearmente con i valori previsti (piuttosto che assumere che la risposta stessa debba variare linearmente). Ad esempio, il caso precedente del numero previsto di partecipanti alla spiaggia sarebbe tipicamente modellato con una distribuzione di Poisson e un collegamento logaritmico, mentre il caso della probabilità prevista di frequentazione della spiaggia sarebbe tipicamente modellato con una distribuzione di Bernoulli (o distribuzione binomiale, a seconda di come viene esattamente formulato il problema) e una funzione di collegamento log-odds (o logit).^[2]

Panoramica

In un GLM, ciascun valore dalla variabile dipendente $\mathbf {Y}$ si assume venga generato da una particolare variabile casuale della famiglia esponenziale, la quale comprende parecchie variabili casuali quali binomiale, poissoniana, gamma, normale inversa e altre. La media ${\boldsymbol {\mu }}$ della distribuzione dipende dalla variabile indipendente $\mathbf {X}$ :

\operatorname {E} (\mathbf {Y} )={\boldsymbol {\mu }}=g^{-1}(\mathbf {X} {\boldsymbol {\beta }}),

dove $\operatorname {E} (\mathbf {Y} )$ è il valore atteso di $\mathbf {Y}$ ; $\mathbf {X} {\boldsymbol {\beta }}$ è il predittore lineare, ovvero una combinazione lineare di $\mathbf {X}$ e parametri ignoti ${\boldsymbol {\beta }}$ ; $g$ è la cosiddetta funzione di collegamento.

In questo ambito, la varianza è tipicamente una funzione $\operatorname {V}$ della media:

\operatorname {Var} (\mathbf {Y} )=\operatorname {V} ({\boldsymbol {\mu }})=\operatorname {V} (g^{-1}(\mathbf {X} {\boldsymbol {\beta }})).

Ciò risulta conveniente se $\operatorname {V}$ è distribuita come una variabile aleatoria della famiglia esponenziale, ma la varianza può essere semplicemente una funzione del valore stimato.

I parametri ignoti ${\boldsymbol {\beta }}$ vengono stimati solitamente con il metodo della massima verosimiglianza, quello della massima quasi-verosimiglianza o con tecniche bayesiane.

Le componenti del modello

Il GLM è composto da tre elementi^[3]:

la funzione di distribuzione $f$ , facente parte della famiglia esponenziale;
il predittore lineare $\eta =\mathbf {X} {\boldsymbol {\beta }}$ ;
una funzione $g$ , detta di collegamento, tale che $\operatorname {E} (\mathbf {Y} )={\boldsymbol {\mu }}=g^{-1}(\eta )$ .

Distribuzione della probabilità

Una famiglia esponenziale iperdispersa di distribuzioni è una generalizzazione di una famiglia esponenziale e il modello di dispersione esponenziale di distribuzioni e include quelle famiglie di distribuzioni di probabilità, con parametri ${\boldsymbol {\theta }}$ e $\tau$ , mentre la funzione di densità $f$ , per il caso di una distribuzione discreta può essere espressa nella forma:

f_{Y}(\mathbf {y} \mid {\boldsymbol {\theta }},\tau )=h(\mathbf {y} ,\tau )\exp \left({\frac {\mathbf {b} ({\boldsymbol {\theta }})^{\rm {T}}\mathbf {T} (\mathbf {y} )-A({\boldsymbol {\theta }})}{d(\tau )}}\right).

Il parametro di dispersione, $\tau$ , tipicamente è noto ed è solitamente correlato alla varianza della distribuzione. Le funzioni $h(\mathbf {y} ,\tau )$ , $\mathbf {b} ({\boldsymbol {\theta }})$ , $\mathbf {T} (\mathbf {y} )$ , $A({\boldsymbol {\theta }})$ , e $d(\tau )$ sono conosciute. Molte distribuzioni comuni appartengono a questa famiglia, tra cui la normale, l'esponenziale, la gamma, la Poisson, la Bernoulli e (per un numero fisso di prove) la binomiale, la multinomiale e la binomiale negativa.

Nel caso in cui $\mathbf {y}$ e ${\boldsymbol {\theta }}$ siano dei semplici scalari $y$ e $\theta$ , il modello diventa:

f_{Y}(y\mid \theta ,\tau )=h(y,\tau )\exp \left({\frac {b(\theta )T(y)-A(\theta )}{d(\tau )}}\right).

Il parametro ${\boldsymbol {\theta }}$ è correlato alla media della distribuzione. Se $\mathbf {b} ({\boldsymbol {\theta }})$ è la funzione identità, si suol dire che la distribuzione è nella forma canonica (o forma naturale). Si noti che qualsiasi distribuzione può essere convertita in forma canonica mediante la sostituzione di ${\boldsymbol {\theta }}$ con ${\boldsymbol {\theta }}'$ per mezzo della trasformazione ${\boldsymbol {\theta }}=\mathbf {b} ({\boldsymbol {\theta }}')$ . È sempre possibile convertire $A({\boldsymbol {\theta }})$ in termini del nuovo parametro ${\boldsymbol {\theta }}'$ , anche se $\mathbf {b} ({\boldsymbol {\theta }}')$ non è una funzione invertibile. Se inoltre, $\mathbf {T} (\mathbf {y} )$ è l'identità e $\tau$ è conosciuto, allora ${\boldsymbol {\theta }}$ è detto parametro canonico (o parametro naturale) ed è correlato alla media dalla relazione

{\boldsymbol {\mu }}=\operatorname {E} (\mathbf {Y} )=\nabla A({\boldsymbol {\theta }}).

Specificazione del modello

Come sopra citato, i modelli lineari generalizzati comprendono una vasta gamma di modelli. Dunque per individuare un particolare modello é necessario specificare:

la funzione di collegamento $g$ (nota anche come link function) da applicare al valore atteso della variabile risposta $g(E[Y])=\eta$ ; tale funzione deve essere nota, monotona e derivabile
la distribuzione della variabile risposta, che deve appartenere alla famiglia di dispersione esponenziale; può essere Bernoulli, Binomiale (risposta dicotomica), Poisson (risposta conteggio), Gamma (fenomeni di durata), Normale, Normale Inversa

Ad esempio: la regressione lineare è un glm con funzione di collegamento la funzione identità e distribuzione della variabile risposta Normale; la regressione logistica è un glm con funzione di collegamento la funzione logit e distribuzione della variabile risposta Bernoulli; la regressione poissoniana è un glm con funzione di collegamento la funzione logaritmica e distribuzione della variabile risposta Poisson^[4].

Stima dei parametri

Stima dei coefficienti di regressione

La stima dei coefficienti di regressione $\beta$ avviene tramite il metodo della massima verosimiglianza, che consiste nel massimizzare la funzione di verosimiglianza. Nel caso dei glm, che conservano l'indipendenza della variabile risposta, tale funzione è data dalla produttoria della funzione di denistà o di probabilità. Si procede facendo la derivata prima della funzione di verosimiglianza e uguagliandola a zero. Da tale equazione si dovrebbe ottenere la stima ${\widehat {\beta }}$ dei coefficienti $\beta$ . Tuttavia tale equazione non ammette soluzione analitica. Pertanto, è necessario ricorrere ad algoritmi numerici per ottenere tale stima.

Il metodo più usato è l'algoritmo IRLS (iterative reweighted leasts squares). Dalla denominazione dell'algoritmo si può desumere il suo funzionamento:

iterative: l'algoritmo itera fino a convergenza l'algoritmo Fisher-Scoring (algoritmo numerico di massimizzazione che deriva dallo sviluppo in serie di Taylor arrestato al primo ordine, che si sintetizza nella seguente formula: $\beta ^{(r+1)}=\beta ^{(r)}+I(\beta ^{(r)})^{-1}*s(\beta ^{(r)})$ , dove $I(\beta ^{(r)})$ è l'informazione attesa di Fisher e $s(\beta ^{(r)})$ è la funzione punteggio o score function.
reweighted: le singole unità statistiche sono riponderate tramite la matrice $W=diag(w_{i})$ dove $w_{i}=\left({\frac {d\mu _{i}}{d\eta _{i}}}\right)^{2}*\left({\frac {1}{Var(Y_{i})}}\right)$
leasts squares: si dimostra che la formula dell'algoritmo di Fisher-Scoring nel caso dei glm si riconduce alla seguente $\beta ^{(r+1)}=(X'W^{(r)}X)^{-1}X'W^{(r)}z^{(r)}$ , dove $z$ è il vettore delle pseudo-risposte $z_{i}=\eta _{i}+(y_{i}-\mu _{i})*{d\eta _{i} \over d\mu _{i}}$ ; che, ponendo $W=I$ e $y_{i}=z_{i}$ , si riconduce alla formula dei minimi quadrati per il modello di regressione lineare normale.^[4]

Stima del parametro di dispersione

Per la stima del parametro di dispersione $\tau$ si ricorre al metodo dei momenti e si ottengono i seguenti risultati:

${\tilde {d(\tau )}}={1 \over n-k-1}*\sum _{i=1}^{N}{(y_{i}-{\widehat {\mu _{i}}})^{2} \over V({\widehat {\mu _{i}}})}$
${\tilde {d(\tau )}}={D \over n-k-1}$ dove $D$ è la devianza del modello.^[4]

Interpretazione dei coefficienti di regressione

Il significato delle stime dei coefficienti di regressione ${\widehat {\beta _{j}}}$ varia a seconda della categoria di modello di modello lineare generalizzato preso in considerazione:

regressioni lineare normale: ${\widehat {\beta _{j}}}$ corrisponde alla variazione assoluta media della variabile risposta in corrispondenza della variazione univariata della j-esima covariata a parità della altre condizioni;
regressione logistica e binomiale: ${\widehat {OR}}=\exp {\beta _{j}}$ corrisponde all'aumento percentuale dell'Odds Ratio all'aumentare di un'unità della j-esima covariata a parità delle altre condizioni
regressione poissoniana: ${\widehat {IRR}}=\exp {\widehat {\beta _{j}}}$ corrisponde alla variazione percentuale della risposta media all'aumentare di un'unità della j-esima covariata a parità delle altre condizioni (Incident Rate Ratio)

Inoltre, è possibile stabilire tramite verifica d'ipotesi la significatività dei singoli coefficienti di regressione. In particolare, si ricorre al test di Wald, noto anche come test Z: $H_{0}:\beta _{j}=0$ vs $H_{1}:\beta _{j}\neq 0$ .

La statistica test è: $z.value={{\widehat {\beta _{j}}}-0 \over SE({\widehat {B_{j}}})}\sim N(0,1)$ , il p-value si calcola come: $2P(Z>|z.value|)$ . Se il p-value è inferiore al valore $\alpha$ fissato a priori allora la covariata è significativa al livello $\alpha$ .

Bonta del modello e confronto tra modelli

Innanzitutto, si definisce il concetto di devianza per un modello lineare generalizzato, interpretabile come la distanza dal modello saturo in termini di estremo superiore della log-verosimiglianza. Il modello saturo è quello con $k=n$ , ha un fit perfetto ma genera overfitting. In formule: $D=-2*d(\tau )*[l({\widehat {\beta }},y)-l({\widehat {\beta }}_{s},y)]=-2*\sum _{i=1}^{N}[y_{i}{\widehat {\theta }}_{i}-c({\widehat {\theta }}_{i})]-[y_{i}{\widehat {\theta }}_{iS}-c({\widehat {\theta }}_{iS})=\sum _{i=1}^{N}d_{i}$ . Inoltre, si introduce la devianza normalizzata $D^{*}={D \over d(\tau )}$

Bontà del modello

Per verificare la bontà del modello si può ricorrere a due test statistici: uno basato sulla devianza ed uno basato sulla $\mathrm {X} ^{2}$ di Pearson. Entrambi hanno come ipotesi nulla l'adeguatezza del modello.

Test basato devianza: la statistica test è $D^{*}{\dot {\sim }}\chi _{n-k-1}^{2}$ (per n grande e parametro di dispersione $\tau$ noto e piccolo^[5]. Per il modello lineare normale la distribuzione del test è esatta. Per il glm Poisson il risultato va bene per $\mu _{i}$ o ${\hat {\mu }}_{i}$ grandi, ad esempio maggiori di 5);
Test basato sulla $\mathrm {X} ^{2}$ di Pearson: la statistica test è $\mathrm {X} ^{2}=\sum _{i=1}^{N}{(y_{i}-{\widehat {\mu }}_{i})^{2} \over \mathrm {V} ({\widehat {\mu }}_{i})d(\tau )}\sim \chi _{n-k-1}^{2}$ . (per n grande).

In entrambi i casi se il p-value è maggiore del livello di significatività fissato a priori, non rifiuto l'ipotesi nulla e concludo che il modello è adeguato.

Confronto tra modelli annidati

Due modelli $M_{1}$ (con $k_{1}$ covariate) ed $M_{2}$ (con $k_{2}$ covariate), tali che $k_{1}<k_{2}$ , si dicono annidati se hanno la stessa specificazione e le prime $k_{1}$ covariate di $M_{2}$ coincidono con le covariate di $M_{1}$ a cui se ne aggiungono altre. Per confrontare due modelli annidati si può ricorrere a due test statistici: il test basato sulla devianza e quello basato sulla statistica F.

Per entrambi i test l'ipotesi nulla è: $H_{0}:\beta _{k_{1}+1}=\beta _{k_{2}+2}=...=\beta _{k_{2}}=0$ .

Le statistiche test sono:

Test basato sulla devianza: $D_{1}^{*}-D_{2}^{*}\sim \chi _{k_{2}-k_{1}}^{2}$
Test basato sulla statistica F: $F={(D_{1}-D_{2})/(k_{2}-k_{1}) \over D_{2}/(n-k_{2}-1)}\sim {\mathcal {F}}_{k_{2}-k_{1},n-k_{2}-1}$ .

In entrambi i test se il p-value è inferiore al valore soglia fissato a priori e scelgo il modello $M_{2}$ .

Confronto tra modelli generici

Per il confronto tra due generici modelli si può ricorrere all' $R^{2}$ di Naglekerke oppure ai criteri AIC (Akaike Information Criterion) e BIC (Bayesian Information Criterion).

$R^{2}={1-exp((-D_{M_{0}}+D_{M})/n) \over 1-exp(-D_{M_{0}}/n)}\in [0,1]$ , il modello migliore ha valore più alto
$AIC=-2*[l_{M}-(k+1)]$ , il modello migliore ha valore più basso
$BIC=-2l_{M}+(k+1)log(n)$ , il modello migliore ha valore più basso

Note

^ J. A. Nelder e R. W. M. Wedderburn, Generalized Linear Models, in Journal of the Royal Statistical Society. Series A (General), vol. 135, n. 3, 1972, pp. 370–384, DOI:10.2307/2344614. URL consultato il 1º novembre 2020.
^ (EN) Stephen Senn, A Conversation with John Nelder, in Statistical Science, vol. 18, n. 1, 2003-02, pp. 118–131, DOI:10.1214/ss/1056397489. URL consultato il 1º novembre 2020.
^ 6.1 - Introduction to Generalized Linear Models | STAT 504, su online.stat.psu.edu. URL consultato il 1º novembre 2020 (archiviato dall'url originale il 27 ottobre 2020).
^ ^a ^b ^c Hardin J.W., Hilbe J.M., Generalized Linear Models and Extensions, Stata Press, 2007.
^ Bent Jørgensen, Exponential Dispersion Models, in Journal of the Royal Statistical Society. Series B (Methodological), vol. 49, n. 2, 1987, pp. 127-162.

Altri progetti

Wikimedia Commons contiene immagini o altri file sul modello lineare generalizzato

Controllo di autorità	Thesaurus BNCF 58064 · BNF (FR) cb167136343 (data)

Portale Economia: accedi alle voci di Wikipedia che trattano di economia

[1] J. A. Nelder e R. W. M. Wedderburn, Generalized Linear Models, in Journal of the Royal Statistical Society. Series A (General), vol. 135, n. 3, 1972, pp. 370–384, DOI:10.2307/2344614. URL consultato il 1º novembre 2020.

[2] (EN) Stephen Senn, A Conversation with John Nelder, in Statistical Science, vol. 18, n. 1, 2003-02, pp. 118–131, DOI:10.1214/ss/1056397489. URL consultato il 1º novembre 2020.

[3] 6.1 - Introduction to Generalized Linear Models | STAT 504, su online.stat.psu.edu. URL consultato il 1º novembre 2020 (archiviato dall'url originale il 27 ottobre 2020).

[:0-4] Hardin J.W., Hilbe J.M., Generalized Linear Models and Extensions, Stata Press, 2007.

[5] Bent Jørgensen, Exponential Dispersion Models, in Journal of the Royal Statistical Society. Series B (Methodological), vol. 49, n. 2, 1987, pp. 127-162.

[1]

[2]

[3]

[4]

[5]

Modello lineare generalizzato

Indice

Intuizione

Panoramica