Vedi:

• Regressione lineare bivariata
• Funzioni per estrarre informazioni dai modelli

Vedi Le formule dei modelli

I regressori (o predittori) del modello sono le variabili indipendenti (i regressori principali $X_j$) e le loro eventuali interazioni (regressori subordinati).

Il numero di coefficienti ($b_j$ o $\beta_j$) è pari al numero di regressori, più l’eventuale intercetta ($b_0$ o $\beta_0$).

modello additivo

$$\hat Y = a + b_1X1 + b_2X2$$

mpg ~ wt + hp

In questo caso, i regressori sono due ($b_1$, $b_2$), i coefficienti sono tre (con $a$).

modello con interazione

$$\hat Y = a + b_1X1 + b_2X2 + b_3X1X2$$

mpg ~ wt * hp

In questo caso, i regressori stimati saranno invece tre: oltre a $b_1$, $b_2$, avremo infatti anche un $b_3$, per l'interazione fra i primi due.

summary(mtcars[c("mpg","hp","wt")])

##       mpg              hp              wt       
##  Min.   :10.40   Min.   : 52.0   Min.   :1.513  
##  1st Qu.:15.43   1st Qu.: 96.5   1st Qu.:2.581  
##  Median :19.20   Median :123.0   Median :3.325  
##  Mean   :20.09   Mean   :146.7   Mean   :3.217  
##  3rd Qu.:22.80   3rd Qu.:180.0   3rd Qu.:3.610  
##  Max.   :33.90   Max.   :335.0   Max.   :5.424  

Gli assunti del modello sono:

• la linearità della relazione;
• indipendenza dei residui;
• omoschedasticità (varianza costante) dei residui;
• normalità della distribuzione dei residui, con media pari a zero;
• Assenza di multicollinearità perfetta: nessuna variabile indipendente nel modello deve essere una combinazione lineare perfetta di una o più altre variabili indipendenti.

plot(mtcars[c("mpg","hp","wt")])

oppure, meglio:

library(psych)
pairs.panels(mtcars[c("mpg","hp","wt")])

che produce il seguente grafico:

Anche se gli assunti di normalità riguardano i residui, e non le variabili del modello, il controllo della normalità è un'analisi utile per anticipare potenziali problemi di eteroschedasticità o di influenza sproporzionata di alcuni dati (asimmetria e outliers).

Utilizziamo la funzione shapiro.test() (vedi test di Shapiro).

Il test conferma che gli outliers della variabile hp determinano uno scostamento dalla normalità.

boxplot(scale(mtcars[c("mpg","hp","wt")]))

• Matrice di correlazione: Si calcola la correlazione (es. Pearson) tra ogni coppia di predittori. Valori molto alti ($|r| > 0.8$ o $0.9$) indicano un problema di collinearità semplice.
• VIF - Variance Inflation Factor: calcolato per ciascun predittore e indica quanto la varianza del suo coefficiente stimato è “gonfiata” dalla presenza degli altri predittori.Un VIF di 1 indica assenza di correlazione.Valori VIF superiori a 5 o, più conservativamente, a 10, sono generalmente considerati problematici e indicano multicollinearità.

cor.test(mtcars$hp, mtcars$wt)

## 	Pearson's product-moment correlation
## 
## data:  mtcars$hp and mtcars$wt
## t = 4.7957, df = 30, p-value = 4.146e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4025113 0.8192573
## sample estimates:
##       cor 
## 0.6587479 

Se si tratta di più variabili, si può creare una matrice di correlazione.

Il test VIF è disponibile in diversi pacchetti di R, dei quali i più diffusi sono DescTools e car. Qui utilizzeremo quest'ultimo. La sintassi è vif(modello):

library(car)
vif(lm(mpg ~ wt + hp, data = mtcars))

##       wt       hp 
## 1.766625 1.766625 

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

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

`

Ovvero: $\hat{mpg} = 32,23 -3,88wt -0,03hp$

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

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

Ovvero: $\hat {mpg} = 49,81 -8,21 wt -0,12 hp +0,03 wt \cdot hp$

Aumentando il numero delle variabili indipendenti, vi è una certa quota di varianza aggiuntiva spiegata, indipendentemente dalla significatività del contributo di ciascuna variabile. Per questa ragione, si utilizza un coefficiente di determinazione “corretto” (Adjusted R-squared), rispetto al numero dei regressori ($p$, parametri):

$$\text{adj.R2} = 1 - (1- R^2)\frac{n-1} {n - p - 1}$$

# per il primo modello
1 - (1 - 0.8268) * (31 / 29)
 
# oppure
1 - (1 - summary(fit1)$r.squared) * (31 / fit1$df.residual)

## [1] 0.8148396

# per il secondo modello
1 - (1 - 0.8848) * (31 / 28)
 
# oppure
1 - (1 - summary(fit2)$r.squared) * (31 / fit2$df.residual)

## [1] 0.872417

Vedi:

• Analisi della varianza e test F nel caso della regressione multipla
• Analisi della devianza nel confronto di modelli

Con la funzione glance() del pacchetto broom (vedi la pagina Broom, e rbind(), è semplice creare una tabella in cui vengono messe a confronto le statistiche dei modelli.

library(broom)
 
rbind("fit1" = glance(fit1), 
      "fit2" = glance(fit2))

## # A tibble: 2 x 11
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC deviance df.residual
## *     <dbl>         <dbl> <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>    <dbl>       <int>
## 1     0.827         0.815  2.59      69.2 9.11e-12     3  -74.3  157.  163.     195.          29
## 2     0.885         0.872  2.15      71.7 2.98e-13     4  -67.8  146.  153.     130.          28

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

Es-regr-mult.R

# analisi esplorativa ------------------
library(psych)
 
summary(mtcars[c("mpg","hp","wp")])
 
# grafici
plot(mtcars[c("mpg","hp","wp")])
pairs.panels(mtcars[c("mpg","hp","wt")])
boxplot(scale(mtcars[c("mpg","hp","wt")]))
 
# test VIF
library(car)
vif(lm(mpg ~ wt + hp, data = mtcars))
 
# analisi di regressione ----------------
 
# modello additivo
fit1 <- lm(mpg ~ wt + hp, data = mtcars) 
summary(fit1)
 
# modello con interazione
fit2 <- lm(mpg ~ wt * hp, data = mtcars) 
summary(fit2)
 
# confronto
library(broom)
rbind("fit1" = glance(fit1), 
      "fit2" = glance(fit2))