Indice

Modelli lineari e scomposizione della devianza

Modelli lineari e scomposizione della devianza

È spesso utile, per valutare i modelli lineari nel dettaglio delle variabili e delle loro interazioni, effettuare l'analisi della varianza (o devianza) spiegata dai singoli regressori.

Caso della regressione lineare semplice

Vedi: Regressione lineare bivariata

Consideriamo il seguente esempio:

fit1 <- lm(data=mtcars, mpg ~ wt)
summary(fit1)

## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528,	Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Tabella dell'Anova

Per ottenere una tabella dell'Anova, usiamo la funzione anova()

anova(fit1)

## Analysis of Variance Table
## 
## Response: mpg
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## wt         1 847.73  847.73  91.375 1.294e-10 ***
## Residuals 30 278.32    9.28                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Per facilitare i calcoli e la scrittura del codice, utilizziamo il pacchetto [[r:tidyverse:pacchetto_broom|]]:

library(broom)
 
tab1 <- tidy(anova(fit1))
tab1

Tab. 1: Caso della regressione semplice

term         df     sumsq    meansq   statistic   p.value
----------  ---  --------  --------  ----------  --------
wt            1   847,725   847,725      91,375         0
Residuals    30   278,322     9,277          --        --

L'analisi della devianza

Tab. 2: I quadrati

term	df	sumsq	meansq
TOTALE	31	1.126,05	36,32	(di mpg)
wt	1	847,73	847,73	spiegata da wt
Residuals	30	278,32	9,28	non spiegata, del modello

Somma dei quadrati (devianza) totale: sum(tab1$sumsq) = 1126.05
Somma dei quadrati spiegata: tab1[[1,3]] = 847.73
Somma dei quadrati residua: tab1[[2,3]] = 278.32

Considerando che

$$F = \frac{\text{varianza spiegata}}{\text{varianza non spiegata}} = \frac{R^2/\text{df.modello}}{(1-R^2)/\text{df.residui}}$$

e che (vedi Regressione lineare bivariata):

$$R^2 = \frac{\text{devianza spiegata}}{\text{devianza totale}}$$

in questo caso, il valore F che troviamo nella tabella è lo stesso del test di significatività del coefficiente di determinazione (in summary(), sopra), in quanto la variabile indipendente è solo una. Infatti:

Il coefficiente di determinazione R² è semplicemente pari alla devianza spiegata da wt (847.73), divisa per la devianza totale (1126.05).

# R^2
R.q <- tab1[[1,3]] / sum(tab1$sumsq)
R.q

## [1] 0.7528328

F è quindi uguale a:

# F = varianza spiegata / varianza non spiegata
(R.q / 1) / ((1-R.q) / 30)

## [1] 91.37533

Caso della regressione multipla

Vedi Regressione lineare multipla (multivariata)

Aggiungiamo una variabile indipendente al modello precedente:

fit2 <- lm(mpg ~ wt + hp, data = mtcars)
summary(fit2)

## 
## Call:
## lm(formula = mpg ~ wt + hp, data = mtcars)
## 
## 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
## 
## Residual standard error: 2.593 on 29 degrees of freedom
## Multiple R-squared:  0.8268,	Adjusted R-squared:  0.8148 
## F-statistic: 69.21 on 2 and 29 DF,  p-value: 9.109e-12

anova(fit2)

## Analysis of Variance Table
## 
## Response: mpg
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## wt         1 847.73  847.73 126.041 4.488e-12 ***
## hp         1  83.27   83.27  12.381  0.001451 ** 
## Residuals 29 195.05    6.73                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La varianza è scomposta nelle quote spiegate da ciascuna variabile. La quota di wt è la stessa di prima (847,73). hp spiega una quota di varianza aggiuntiva pari al 8,9% del totale.

Aggiungiamo come terzo regressore l'interazione fra le due variabili indipendenti:

fit3 <- lm(mpg ~ wt * hp, data = mtcars) 
summary(fit3)

## 
## Call:
## lm(formula = mpg ~ wt * hp, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0632 -1.6491 -0.7362  1.4211  4.5513 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 49.80842    3.60516  13.816 5.01e-14 ***
## wt          -8.21662    1.26971  -6.471 5.20e-07 ***
## hp          -0.12010    0.02470  -4.863 4.04e-05 ***
## wt:hp        0.02785    0.00742   3.753 0.000811 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.153 on 28 degrees of freedom
## Multiple R-squared:  0.8848,	Adjusted R-squared:  0.8724 
## F-statistic: 71.66 on 3 and 28 DF,  p-value: 2.981e-13

Costruiamo la tabella dell'Anova:

tab2 <- tidy(anova(fit3))
tab2

Tab. 3: Caso della regressione multipla

term         df     sumsq    meansq   statistic   p.value
----------  ---  --------  --------  ----------  --------
wt            1   847,725   847,725     182,923     0,000
hp            1    83,274    83,274      17,969     0,000
wt:hp         1    65,286    65,286      14,088     0,001
Residuals    28   129,761     4,634          --        --

Il terzo regressore spiega una ulteriore — e statisticamente significativa — quota di varianza.

Analisi della devianza

La tabella che segue sintetizza il calcolo della devianza e della varianza spiegata e totale del modello:

Tab. 4: I quadrati

termine	df	sumsq	meansq = sumsq/df
TOTALE	31	1126,05	36,32	= var(mpg)
wt	1	847,73	847,73	spiegata da wt
hp	1	83,27	83,27	spiegata da hp
wt:hp	1	65,29	65,29	spiegata da wt:hp
SPIEGATA	3	996,29	332,10	spiegata dal modello
Residua	28	129,76	4,63	non spiegata, del modello

somma dei quadrati e quadrati medi per i regressori sono uguali (meansq = sumsq/df);
la devianza spiegata dal modello è pari alla somma delle devianze: 847,73 + 83,27 + 65,29 = 996,29;
la varianza spiegata dal modello è pari alla devianza spiegata, divisa per i gradi di libertà del modello (3): 996,29 / 3 = 332,095;
i gradi di libertà del modello (devianza spiegata) sono 3 (i regressori);
i gradi di libertà dei residui sono il numero dei casi meno il numero delle variabili (32 - 4 = 28).

Il coefficiente di determinazione R² è dunque uguale a 996,29 (devianza spiegata) / 1126,05 (devianza totale). Con riferimento agli indici della tab2:

(tab2[[1,3]] + tab2[[2,3]] + tab2[[3,3]]) / sum(tab2$sumsq)

o, più semplicemente:

R.q <- 1 - (tab2[[4,3]] / sum(tab2$sumsq))
R.q

## [1] 0.8847637

In tabella 3 non troviamo il valore F relativo al coefficiente di determinazione, pari a:

(R.q / 3) / ((1-R.q) / 28)

## [1] 71.65967

Analisi della devianza nel confronto di modelli

Avendo chiari i significati di queste misure, e i rapporti fra queste grandezze, è possibile utilizzare in maniera corretta l'analisi della devianza (la tabella prodotta da anova(), che in questo caso non presenta i quadrati medi) per confrontare due o più modelli alternativi, purché riferiti allo stesso dataset (nested models)¹⁾:

anova(fit1, fit2, fit3)

## Analysis of Variance Table
## 
## Model 1: mpg ~ wt
## Model 2: mpg ~ wt + hp
## Model 3: mpg ~ wt * hp
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     30 278.32                                  
## 2     29 195.05  1    83.274 17.969 0.0002207 ***
## 3     28 129.76  1    65.286 14.088 0.0008108 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vengono qui confrontate le differenze fra le devianze dei tre modelli lineari:

278,32 - 195,05 = 83,274 (devianza spiegata in più dal secondo modello rispetto al primo, ovvero la devianza spiegata da hp);
195,05 - 129,76 = 65,286 (devianza spiegata in più dal terzo modello rispetto al secondo, ovvero la devianza spiegata da wt:hp).

Come si vede, i modelli sono concatenati, ovvero vengono confrontati progressivamente.

Il valore F dei modelli è riferito alla devianza più piccola (e non corrisponde ai valori F dei regressori):

83,274 / (129,76 / 28);
65,286 / (129,76 / 28).

L'ipotesi nulla è che ciascuna di queste differenze sia statisticamente non significativa. In questo caso, l'ipotesi nulla può essere respinta (le differenze sono statisticamente significative).

Cambiando l'ordine dei modelli:

anova(fit3, fit1, fit2)

## Analysis of Variance Table
## 
## Model 1: mpg ~ wt * hp
## Model 2: mpg ~ wt
## Model 3: mpg ~ wt + hp
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     28 129.76                                  
## 2     30 278.32 -2  -148.560 16.028 2.293e-05 ***
## 3     29 195.05  1    83.274 17.969 0.0002207 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

si conclude che il secondo modello spiega significativamente meno del primo (-148.560), e il terzo significativamente più del primo. I valori di F sono sempre riferiti a 129,76, e sono uguali a:

(-148,560 / -2) / (129,76 / 28);
83,274 / (129,76 / 28).

Con i Modelli lineari generalizzati

Con i risultati dei modelli lineari generalizzati, la funzione anova() produce direttamente la tabella dell'analisi della devianza:

anova(glm(data=mtcars, mpg ~ wt))

## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: mpg
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev
## NULL                    31    1126.05
## wt    1   847.73        30     278.32

Per avere il test F, scriveremo:

anova(glm(data=mtcars, mpg ~ wt), 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 91.375 1.294e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tab. 5: I quadrati (glm)

	Devianza
di mpg (modello 0)	1126,05
spiegata da wt	847,73	sistematica
residua (del modello)	278,32	casuale

Anche per il confronto fra più modelli, utilizzeremo lo stesso comando:

gfit1 <- glm(mpg ~ wt + hp, data = mtcars) 
gfit2 <- glm(mpg ~ wt * hp, data = mtcars) 
anova(gfit1, gfit2, test = "F")

## 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

Script di esempio

Es-anova-lm-glm.R

 
# ANOVA DEI MODELLI LINEARI
 
# primo modello
fit1 <- lm(data=mtcars, mpg ~ wt)
summary(fit1)
anova(fit1)
 
# tabella 1
library(broom)
tab1 <- tidy(anova(fit1))
tab1
 
# secondo modello
fit2 <- lm(mpg ~ wt + hp, data = mtcars)
summary(fit2)
anova(fit2)
 
# terzo modello
fit3 <- lm(mpg ~ wt * hp, data = mtcars) 
summary(fit3)
 
# tabella 2
tab2 <- tidy(anova(fit3))
tab2
 
 
# ANOVA per confrontare modelli
 
anova(fit1, fit2, fit3)
 
anova(fit3, fit1, fit2)
 
# nei modelli lineari generalizzati
anova(glm(data=mtcars, mpg ~ wt))
 
anova(glm(data=mtcars, mpg ~ wt), test = "F")
 
gfit1 <- glm(mpg ~ wt + hp, data = mtcars) 
gfit2 <- glm(mpg ~ wt * hp, data = mtcars) 
anova(gfit1, gfit2, test = "F")
 
rm(fit1, fit2, fit3, gfit1, gfit2, tab1, tab2, R.q)

Analisi multivariata, Regressione lineare, Modelli lineari

¹⁾

Si tratta del requisito di omoschedasticità; bisogna anche fare attenzione ai casi mancanti, che potrebbero incidere sul numero dei casi e i gradi di libertà dei modelli