I diversi tipi di modelli di regressione possono essere ricondotti ad una forma generale, chiamata modello lineare generalizzato¹⁾.

Poiché la variabile dipendente può anche non avere una distribuzione normale, i modelli lineari generalizzati si prestano all’analisi esplorativa dei dati (data mining). La funzione di R glm() può essere utilizzata in diversi contesti di analisi:

• stime sulla massima verosimiglianza (MLE);
• stime sulla massima probabilità a posteriori (MAP, ML con a priori) e stime bayesiane;
• confronti basati su information criteria.

La prima equazione è relativa alla componente sistematica:

\begin{equation} \label{eq:theta} \theta = g(\eta) \end{equation}

Dove $\theta$ è il valore atteso (ad es.: $\hat Y$), e $\eta$ è il predittore lineare:

\begin{equation} \label{eq:eta} \eta = \beta_0 + \beta_j X_j \end{equation}

Nella prima equazione, $g$ è la funzione inversa del link ($L$) che linearizza la relazione, in modo che $L(\theta) = \eta$.

La terza equazione è relativa alla componente casuale o stocastica, ed è indicata nella funzione di R dalla distribuzione di probabilità (con l'argomento family).

In pratica: Per generalizzare il modello lineare, si applica la funzione di link $L$. La funzione inversa di $L$ ($g$) "riporta indietro" il predittore alla scala del valore atteso: ad esempio exp() per il link log().

Nel caso della regressione lineare, avremo:

$$\theta = \eta$$

Il risultato del modello è già nella scala dei valori attesi:

$$\eta = \mu = \beta_0 + \beta_j X_j$$

Nel caso della regressione logistica avremo invece:

$$\eta = ln(\omega) = \beta_0 + \beta_j X_j$$

dove $\omega$ è l'odd e $ln(\omega)$ è il logit.

I valori attesi (gli odds) possono essere calcolati con l'antilogaritmo, ovvero la funzione exp(), e poi riportati alle probabilità $\pi$ con la funzione logistica.

Nel caso dei modelli log-lineari, infine, avremo:

$$\theta = \pi = e^{\eta}$$

dove il logaritmo naturale (funzione log()) linearizza il modello, e la funzione exp() riporta il risultato alla scala dei valori attesi ($\pi$).

La notazione dei modelli lineari generalizzati ci consente di leggere l'aiuto di R sulla funzione, e di applicarla a seconda delle necessità:

glm(formula, 
    data, subset, na.action,          # dati
    family = gaussian,                # distribuzione
    weights,                                    
    start = NULL, 
    etastart,                                   
    mustart,                                    
    ...)

Per confrontare l'output della funzione glm() con l'output di lm(), riprendiamo l'esempio presentato in Regressione lineare bivariata:

glm.res <- glm(dist ~ speed, data = cars)
 
summary(glm.res)

## Call:
## glm(formula = dist ~ speed, data = cars)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -29.069   -9.525   -2.272    9.215   43.201  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 236.5317)
## 
##     Null deviance: 32539  on 49  degrees of freedom
## Residual deviance: 11354  on 48  degrees of freedom
## AIC: 419.16
## 
## Number of Fisher Scoring iterations: 2

Il risultato è lo stesso per quanto riguarda i coefficienti ($\eta$), ma non viene fornito il valore di $R^2$.

Ciò che andiamo a vedere come prima cosa è la differenza — o il rapporto — fra devianza totale (Null deviance) e devianza residua (Residual deviance), che fornisce la misura della capacità esplicativa del modello.

Il rapporto fra queste due grandezze corrisponde infatti al coefficiente di determinazione:

# R^2 = 1 - (Residual deviance / Null deviance)
1 - (deviance(glm.res)/glm.res$null.deviance)

## [1] 0.6510794

Nel summary, inoltre, come misura di errore viene riportata la dispersione, anziché l'errore standard dei residui (Dispersion parameter for gaussian family taken to be 236.5317).

In questo caso, le devianze sono calcolate rispetto alla distribuzione normale, e quindi la dispersione è pari alla varianza dei residui, ovvero la devianza residua divisa per il numero di gradi di libertà residui, ed equivale al quadrato dell'errore standard dei residui (vedi I residui):

# errore standard dei residui al quadrato
sigma(glm.res)^2

## [1] 236.5317

In generale, però, si tratta di una misura della variabilità dei dati rispetto al modello, il cui calcolo dipende dalla distribuzione (family) e dal link adottati.

Come detto , nei GLM, l'analisi della devianza assume un ruolo centrale per la valutazione della bontà di adattamento. La devianza nulla e la devianza residua del modello permettono infatti di quantificare e valutare la variabilità spiegata.

La significatività di tale quota viene valutata attraverso test statistici specifici, che, anche in questo caso, dipendono dal modello applicato (o, per meglio dire, dal link e dalla family).

In questo caso, per controllare la significatività statistica della devianza spiegata, si ricorre al test “F”, costruendo la tabella dell’analisi della devianza, con la funzione anova():

anova(glm.res, test = "F")

## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: dist
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev      F   Pr(>F)    
## NULL                     49      32539                    
## speed  1    21186        48      11354 89.567 1.49e-12 ***
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

L'ipotesi nulla è che la devianza residua del modello zero (NULL) e quella del modello con un regressore non siano significativamente diverse.

Non sempre il test F è appropriato. Quando non lo è (ad esempio, per famiglie di distribuzioni come la binomiale o la poissoniana, o quando si confrontano modelli basati sulla devianza), si utilizzerà un test del chi-quadrato (o un test del rapporto di verosimiglianza, LRT, che spesso segue una distribuzione chi-quadrato).

# ad esempio
anova(glm.res, test = "Chisq")

Vedi Information criteria

Stesse differenze possono essere osservate applicando la funzione glm() ai modelli utilizzati come esempio alla pagina Regressione lineare multipla (multivariata):

fit1 <- glm(mpg ~ wt + hp, data = mtcars)
summary(fit1)

## Call:
## glm(formula = mpg ~ wt + hp, data = mtcars)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -3.941  -1.600  -0.182   1.050   5.854  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.22727    1.59879  23.285  < 2e-16 ***
## wt          -3.87783    0.63273  -6.129 1.12e-06 ***
## hp          -0.03177    0.00903  -3.519  0.00145 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 6.725785)
## 
##     Null deviance: 1126.05  on 31  degrees of freedom
## Residual deviance:  195.05  on 29  degrees of freedom
## AIC: 156.65
## 
## Number of Fisher Scoring iterations: 2

La tabella per il confronto delle statistiche di due modelli presenta i valori della devianza e della divergenza, insieme agli Information Criteria:

# secondo modello (moltiplicativo)
fit2 <- glm(mpg ~ wt * hp, data = mtcars) 
 
rbind(glance(fit1), glance(fit2))

## # A tibble: 2 x 7
##   null.deviance df.null logLik   AIC   BIC deviance df.residual
## *         <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int>
## 1         1126.      31  -74.3  157.  163.     195.          29
## 2         1126.      31  -67.8  146.  153.     130.          28

Vediamo che:

• il secondo modello spiega una quota maggiore della variabilità totale, avendo una devianza residua inferiore. La differenza tra queste devianze residue ($195,1 - 129,8 = 65,3$) è attribuibile al termine di interazione; di conseguenza:
• la divergenza (dalla distribuzione di massima verosimiglianza) è minore, così come i valori degli information criteria.

Per confrontare due modelli annidati (nested: stimati sugli stessi dati, ma con parametri diversi), possiamo tornare ad usare la funzione anova(), che calcola la devianza residua e il test F per il confronto fra i modelli.

A questo proposito, è necessario ricordare che:

• la validità del p-value di questo test F dipende dal rispetto delle assunzioni del modello più generale (inclusa l'omoschedasticità dei residui);
• il test F è applicabile solo a modelli lineari generalizzati con distribuzione gaussiana e link identità. Per altri modelli, si utilizzerà un test del chi-quadrato.

anova(fit1, fit2, test = "F")
# in questo caso, dobbiamo esplicitare il test

## Analysis of Deviance Table
## 
## Model 1: mpg ~ wt + hp
## Model 2: mpg ~ wt * hp
##   Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
## 1        29     195.05                                 
## 2        28     129.76  1   65.286 14.088 0.0008108 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Applicando il test al secondo modello (quello con più regressori) possiamo valutare la rilevanza dei regressori singolarmente presi:

anova(fit2, test = "F")

## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: mpg
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev       F    Pr(>F)    
## NULL                     31    1126.05                      
## wt     1   847.73        30     278.32 182.923 8.428e-14 ***
## hp     1    83.27        29     195.05  17.969 0.0002207 ***
## wt:hp  1    65.29        28     129.76  14.088 0.0008108 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Il test F viene applicato sequenzialmente alla devianza residua dopo l'introduzione di ciascun regressore. Di conseguenza, il risultato dipende anche dall'ordine delle variabili nel modello (vedi: Modelli lineari e scomposizione della devianza).

E' possibile scaricare ed eseguire lo script dell'esempio:

es_glm.R

 
# REGRESSIONE LINEARE SEMPLICE ------------------------------------------------
glm.res <- glm(dist ~ speed, data = cars) 
summary(glm.res)
 
# test F
anova(glm.res, test = "F")
 
# tabella delle statistiche
library(broom)
glance(glm.res)
 
# REGRESSIONE LINEARE MULTIPLA ------------------------------------------------
fit1 <- glm(mpg ~ wt + hp, data = mtcars) 
summary(fit1)
 
fit2 <- glm(mpg ~ wt * hp, data = mtcars) 
 
# tabella per il confronto delle statistiche
rbind(glance(fit1), glance(fit2))
 
# Test F
anova(fit1, fit2, test = "F")
 
# Test F sui regressori
anova(fit2, test = "F")