Lab. 5 - Regressione multipla 2
Giovanni Marchetti
Predittori binari
Le variabili binarie sono spesso codificate con \(0\), \(1\) oppure \(-1\), \(+1\). Se una variabile binaria è invece rappresentata come un fattore con due livelli (vedi dopo),il fattore nella formula del modello viene trasformato in una variabile binaria. In tal caso il valore 0 è assegnato al primo livello del fattore e il valore 1 al secondo.
Matrice del modello
Una volta definito il modello è possibile estrarre la matrice del modello con la funzione model.matrix. Per esempio, supponiamo di avere i dati dell’esperimento sull’effetto di un trattamento di isolamento sul consumo di gasolio da riscaldamento. Si conosce anche come variabile concomitante la temperatura.
library(MASS)
attach(whiteside)
head(whiteside) Insul Temp Gas
1 Before -0.8 7.2
2 Before -0.7 6.9
3 Before 0.4 6.4
4 Before 2.5 6.0
5 Before 2.9 5.8
6 Before 3.2 5.8consumo = Gas * 28.3168466 # scala in metri cubi
isolamento = Insul
temperatura = TempLa matrice del modello consumo ~ isolamento + temperatura
model.matrix( ~ isolamento + temperatura)[24:30,] (Intercept) isolamentoAfter temperatura
24 1 0 8.5
25 1 0 9.1
26 1 0 10.2
27 1 1 -0.7
28 1 1 0.8
29 1 1 1.0
30 1 1 1.4Notate che la codifica binaria scelta per il fattore isolamento assegna lo 0 al primo livello Before e l’1 al secondo livello After.
Ecco lo scatter separato per i due gruppi di trattamento:
plot(temperatura, consumo, pch = 20, col = c('red', 'blue')[isolamento])
Adattamento del modello delle rette parallele:
mac = lm(consumo ~ isolamento + temperatura)
summary(mac)
Call:
lm(formula = consumo ~ isolamento + temperatura)
Residuals:
Min 1Q Median 3Q Max
-21.02 -6.31 1.23 6.90 21.04
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 185.513 3.344 55.5 <2e-16 ***
isolamentoAfter -44.322 2.748 -16.1 <2e-16 ***
temperatura -9.534 0.503 -18.9 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.1 on 53 degrees of freedom
Multiple R-squared: 0.91, Adjusted R-squared: 0.906
F-statistic: 267 on 2 and 53 DF, p-value: <2e-16b = coef(mac)Grafico delle due rette parallele adattate:
plot(temperatura, consumo, pch = 20, col = c('red', 'blue')[isolamento])
abline(b[1], b[3], col = "red", lwd = 2)
abline(b[1]+b[2], b[3], col = "blue", lwd = 2)
Modelli con Interazioni
Una interazione è una variabile ottenuta moltiplicando due predittori \(X_j \cdot X_h\).
In un modello l’interazione si definisce con uno speciale operatore : nella formula.
Esempio del consumo di gasolio prima e dopo il trattamento di isolamento
Adattamento del modello con interazione tra temperatura e trattamento:
maci = lm(consumo ~ isolamento + temperatura + isolamento:temperatura)
b = coef(maci)
summary(maci)
Call:
lm(formula = consumo ~ isolamento + temperatura + isolamento:temperatura)
Residuals:
Min 1Q Median 3Q Max
-27.70 -5.10 1.06 5.93 18.07
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 194.079 3.850 50.41 < 2e-16 ***
isolamentoAfter -60.314 5.100 -11.83 2.3e-16 ***
temperatura -11.135 0.637 -17.49 < 2e-16 ***
isolamentoAfter:temperatura 3.265 0.909 3.59 0.00073 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.15 on 52 degrees of freedom
Multiple R-squared: 0.928, Adjusted R-squared: 0.924
F-statistic: 222 on 3 and 52 DF, p-value: <2e-16Grafico delle due rette adattate:
plot(temperatura, consumo, pch = 20, col = c('red', 'blue')[isolamento])
abline(b[1], b[3], col = "red", lwd = 2)
abline(b[1]+b[2], b[3]+b[4], col = "blue", lwd = 2)
Predittori qualitativi
I predittori qualitativi o fattori sono codificati con più variabili indicatrici delle sue modalità . Spesso le modalità di un fattore si chiamano livelli. Le variabili indicatrici si chiamano anche dummy variables.
Per esempio se \(G\) è una variabile qualitativa con 3 livelli si creano \(3\) variabili dummy \(X_j\), \(j = 1, \dots, 3\) definite da \[
X_j = \mathcal{I}(G = j)
\] Ogni \(X_j\) è un vettore \(n \times 1\) di 1 e 0 in cui gli 1 indicano le unità con la modalità \(j\) del fattore. Mettiamo che \(G\) sia il tipo di scuola di uno studente con 4 livelli - L liceo - T istituto tecnico - P istituto professionale - A altro istituto
Abbiamo 8 studenti, due per ogni tipo di scuola.
La variabile \(G\) è definita come segue in R
G = c("L", "L", "T", "T", "P", "P", "A", "A")
G = factor(G)
G[1] L L T T P P A A
Levels: A L P TNota che i livelli sono in ordine alfabetico.
Quando il fattore compare in un modello viene trasformato nelle dummy.
model.matrix(~ G)[,] (Intercept) GL GP GT
1 1 1 0 0
2 1 1 0 0
3 1 0 0 1
4 1 0 0 1
5 1 0 1 0
6 1 0 1 0
7 1 0 0 0
8 1 0 0 0Tuttavia siccome la somma delle dummy è sempre 1 e quindi è indentica alla prima colonna di 1, viene eliminata automaticamente la prima cioè quella associata al primo livello (in questo caso A). Questo è il livello di riferimento.
Volento cambiare il livello di riferimento per esempio predendo L come primo livello si usa il comando relevel
G = relevel(G, ref = 'L')
G[1] L L T T P P A A
Levels: L A P Tmodel.matrix(~ G)[,] (Intercept) GA GP GT
1 1 0 0 0
2 1 0 0 0
3 1 0 0 1
4 1 0 0 1
5 1 0 1 0
6 1 0 1 0
7 1 1 0 0
8 1 1 0 0Analisi del punteggio al test d’ingresso a Economia
Lettura dei dati
x = read.csv('~/Dropbox/Public/gmm/static/R/testingresso.csv')
summary(x) sesso prov comune
F:513 Min. : 1 Min. : 1.0
M:615 1st Qu.: 48 1st Qu.: 8.0
Median : 48 Median : 17.0
Mean : 56 Mean : 20.1
3rd Qu.: 51 3rd Qu.: 24.0
Max. :104 Max. :272.0
istituto
Liceo Scientifico :456
Istituto Tecnico Commerciale :247
Liceo Classico : 82
Liceo Linguistico : 71
Istituto Tecnico Commerciale e per Geometri: 40
Istituto Tecnico per Geometri : 32
(Other) :200
tipo annomat votomat
Liceo Scientifico :456 Min. : 0 Min. : 60.0
Istituto Tecnico Commerciale:341 1st Qu.:2013 1st Qu.: 68.0
Liceo Classico : 82 Median :2013 Median : 76.0
Liceo Linguistico : 71 Mean :2005 Mean : 77.2
Istituto Tecnico Industriale: 70 3rd Qu.:2013 3rd Qu.: 84.0
Altro istituto : 38 Max. :2013 Max. :100.0
(Other) : 70 NA's :1
score idtest
Min. :-1.0 Min. :28
1st Qu.:10.5 1st Qu.:28
Median :13.2 Median :28
Mean :13.3 Mean :28
3rd Qu.:16.0 3rd Qu.:28
Max. :24.0 Max. :28
Definizione delle variabili e dei fattori
VOTOMAT = x$votomat
PUNTEGGIO = x$score
SESSO = x$sesso#[sel]
SCUOLA = x$istituto
levels(SCUOLA) = c("A", "A", "A", "A", "A", "P", "P", "P", "P", "P", "P", "P", "P", "T","T","T","T","T","T","T","T","L","L", "L", "L", "L", "A", "A", "A", "A")
SCUOLA = relevel(SCUOLA, 'L')
test = data.frame(VOTOMAT, PUNTEGGIO, SESSO, SCUOLA)
class(SESSO)[1] "factor"levels(SESSO)[1] "F" "M"class(SCUOLA)[1] "factor"levels(SCUOLA)[1] "L" "A" "P" "T"Analisi descrittiva
Il punteggio dipende dal voto alla maturità ? e dal sesso?
plot(VOTOMAT, PUNTEGGIO, type = "n")
points(VOTOMAT[SESSO == "M"], PUNTEGGIO[SESSO=="M"], pch = 20, col = "RoyalBlue")
points(VOTOMAT[SESSO == "F"], PUNTEGGIO[SESSO=="F"], pch = 20, col = "Orange")
Box-plot affiancati
boxplot(PUNTEGGIO ~ SESSO, horizontal = TRUE, col = c("pink", "cyan"))
Modello di regressione multipla
mo = lm(PUNTEGGIO ~ VOTOMAT + SESSO)
summary(mo)
Call:
lm(formula = PUNTEGGIO ~ VOTOMAT + SESSO)
Residuals:
Min 1Q Median 3Q Max
-12.81 -2.41 0.00 2.49 11.42
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.6356 0.8240 5.63 2.3e-08 ***
VOTOMAT 0.0981 0.0102 9.65 < 2e-16 ***
SESSOM 2.0529 0.2266 9.06 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.72 on 1125 degrees of freedom
Multiple R-squared: 0.116, Adjusted R-squared: 0.114
F-statistic: 73.5 on 2 and 1125 DF, p-value: <2e-16L’effetto del punteggio alla maturità è piuttosto modesto anche se significativo.
Il punteggio dei maschi è in media più alto di 2 punti a parità di voto alla maturità .
Benché l’\(R^2\) sia solo circa il \(12\%\) il test \(F\) indica che le due variabili sono rilevanti.
Il test \(F\) confronta il modello completo con il modello che ha solo la costante.
In dettaglio
mo0 = lm(PUNTEGGIO ~ 1)
anova(mo0, mo, test = "F")Analysis of Variance Table
Model 1: PUNTEGGIO ~ 1
Model 2: PUNTEGGIO ~ VOTOMAT + SESSO
Res.Df RSS Df Sum of Sq F Pr(>F)
1 1127 17588
2 1125 15556 2 2032 73.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La statistica test è \[ F = \frac{(RSS_0 - RSS)/(df_0 - df)}{RSS/df} = \frac{(17588-15556)/2}{15556/1125} = 73.5 \] Sotto \(H_0\) la statistica ha distribuzione \(F_{2,1125}\). Al livello del \(5\%\) si rifiuta \(H_0: \beta_1 = \beta_2 = 0\) se \(F > f_{0.95}\) cioè se supera
qf(0.95, 2, 1125)[1] 3.004Ovviamente si rifiuta l’ipotesi perché \(f_{oss} = 73.5\).
Modello che include la scuola di provenienza.
Studiamo la matrice del modello.
mo2 = lm(PUNTEGGIO ~ VOTOMAT + SCUOLA)
head(model.matrix(mo2)) (Intercept) VOTOMAT SCUOLAA SCUOLAP SCUOLAT
1 1 77 1 0 0
2 1 73 1 0 0
3 1 72 1 0 0
4 1 60 1 0 0
5 1 100 1 0 0
6 1 60 1 0 0Interpretazione del modello. I valori medi sono \[ E( Y\mid X) = \beta_0 + \beta_1 X_1 + (\beta_2 X_2 + \beta_3 X_3 + \beta_4 X_4) \] quindi l’effetto della scuola è quello di far cambiare l’intercetta: \[ E(Y \mid X) = \begin{cases} \beta_0 + \beta_1 X_1 & \text{per il liceo}\\ (\beta_0 + \beta_2) + \beta_1 X_1 &\text{per le altre scuole}\\ (\beta_0 + \beta_3) + \beta_1 X_1 &\text{per il professionale}\\ ( \beta_0 + \beta_4) + \beta_1 X_1 &\text{per il tecnico}\\ \end{cases} \]
Risultati
summary(mo2)
Call:
lm(formula = PUNTEGGIO ~ VOTOMAT + SCUOLA)
Residuals:
Min 1Q Median 3Q Max
-13.61 -2.61 0.11 2.44 9.57
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.97896 0.76879 9.08 < 2e-16 ***
VOTOMAT 0.09533 0.00997 9.56 < 2e-16 ***
SCUOLAA -1.72214 0.63098 -2.73 0.0064 **
SCUOLAP -3.62494 0.48787 -7.43 2.1e-13 ***
SCUOLAT -2.05353 0.23624 -8.69 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.68 on 1123 degrees of freedom
Multiple R-squared: 0.136, Adjusted R-squared: 0.133
F-statistic: 44.2 on 4 and 1123 DF, p-value: <2e-16A parità di voto alla maturità in media il punteggio degli studenti del professionale sono 3 punti più bassi di quelli degli studenti del liceo. L’effetto del voto alla maturità è sempre positivo ma debole.
L’adattamento è \(R^2_{adj} = 0.13\) mentre quello del modello precedente con i predittori VOTOMAT e SESSO era \(R^2_{adj} = 0.11\).
Una figura aiuta a capire.
plot(VOTOMAT, PUNTEGGIO, pch = 20, col = "gray")
beta = coef(mo2)
abline(beta[1], beta[2], col = "red", lwd = 3)
abline(beta[1]+beta[3], beta[2], col = "green", lwd = 3)
abline(beta[1]+beta[4], beta[2], col = "RoyalBlue", lwd = 3)
abline(beta[1]+beta[5], beta[2], col = "orange", lwd = 3)
legend("bottomright", c("Liceo", "Altro","Prof","Tecnico"), col = c("red", "green", "RoyalBlue","orange"), lty = c(1, 1, 1, 1), cex = 0.5, lwd = c(3 , 3, 3, 3))
Esercizi 5
Vero o falso: la retta di regressione ha pendenza nulla se le variabili sono incorrelate. Giustificare.
Vero o falso: se il coefficiente di correlazione è \(−1\) la retta di regressione ha pen- denza negativa con inclinazione di 45 gradi. Gustificare.
Vero o falso: I valori adattati hanno una somma uguale a \(n\bar y\)Ì„. Giustificare.
Vero o falso: se le due varianze campionarie di X e Y sono uguali a 1,il coefficiente di regressione e il coefficiente di correlazione sono uguali.
Vero o falso: la devianza residua è sempre più piccola della devianza di regressione. Giustificare.
Vero o falso: Se l’indice di determinazione lineare è zero è certo che Y non dipende in alcun modo da X. Giustificare.
Vero o falso: Se il coefficiente di correlazione è superiore a 0.9 allora le ipotesi del modello di regressione lineare sono verificate. Giustificare.
Vero o falso: se l’indice di determinazione lineare è zero allora le due variabili sono indipendenti. Giustificare.
Vero o falso: se due variabili esplicative sono esattamente funzioni lineari una dell’altra non esistono le stime dei minimi quadrati.
Vero o falso: se la risposta è esattamente funzione lineare di una o più variabili esplicative non esistono le stime dei minimi quadrate dei coefficienti.
Quattro oggetti vengono pesati da soli e in combinazione ottenendo i risultati seguenti: \[ \begin{matrix} \textrm{Oggetti pesati} & 1 & 2 & 3 & 4 & 1\&2 & 2\&3 & 3\&4 & 1\&4 \\ \textrm{Peso}& 21.08 &29.78 & 20.33 &111.49 & 41.43 & 39.71 & 118.77 & 120.34 \\ \end{matrix} \] Definire un modello di regressione opportuno e stimare con i minimi quadrati i pesi dei quattro oggetti.
Analizzare i dati sul test di ingresso a Economia includendo tutti i predittori ritenuti utili. Commentare con l’analisi svolta a lezione.
Soluzioni
Vero perché \(\hat \beta_1 = r_{xy} \frac{s_y}{s_x}\)
Falso perché se \(r = -1\) i punti \((x_i, y_i)\) sono allineati su una retta di pendenza qualsiasi ma negativa.
Vero, nel modello di regressione lineare semplice la somma dei \(y_i\) è uguale alla somma dei \(\hat y_i\).
Vero per la 1.
Falso perché \[ TSS = RSS + ESS \] ma questo implica solo che \(RSS \le TSS\).
Falso. Se R^2 = 0$ allora \(r_{xy} = 0\) e quindi \(X\) e \(Y\) sono solo incorrelati, ma in generale non indipendenti.
Falso. Può capitare che \(r_{xy} > 0.9\) ma ciò nonostante il modello sia non lineare.
Falso come in 6.
Vero. In tal caso il rango di \(X\) non è pieno e non esiste una sola soluzione del sistema di equazioni normali.
Falso. Per esempio
y = 1:5
x2 = c(2, 7, 8, 10, 20)
y = 5 + 3* x2 + rnorm(5)
x1 = 3*y
plot(x1, y)
lm(y ~ x1 + x2)
Call:
lm(formula = y ~ x1 + x2)
Coefficients:
(Intercept) x1 x2
1.91e-14 3.33e-01 4.65e-15 Quindi le stime esistono ma gli altri parametri sono zero.
- Calcoli.
pesi = c(21.08, 29.78, 20.33, 111.49, 41.43, 39.71, 118.77, 120.34)
X = matrix(
c(1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1,
1, 1, 0, 0,
0, 1, 1, 0,
0, 0, 1, 1,
1, 0, 0, 1), 8, 4, byrow= TRUE)
b = solve(t(X) %*% X) %*% t(X) %*% pesi
b [,1]
[1,] 16.87
[2,] 26.18
[3,] 15.52
[4,] 106.07- Istruzioni
x = read.csv('~/Dropbox/Public/gmm/static/R/testingresso.csv')
VOTOMAT = x$votomat
PUNT = x$score
SESSO = x$sesso#[sel]
SCUOLA = x$istituto
LSC = 1 + (x$tipo == "Liceo Scientifico")
levels(SCUOLA) = c("A", "A", "A", "A", "A", "P", "P", "P", "P", "P", "P", "P", "P", "T","T","T","T","T","T","T","T","L","L", "L", "L", "L", "A", "A", "A", "A")
SCUOLA = relevel(SCUOLA, 'L')
ANNO = x$annomat
dati = data.frame(VOTOMAT, PUNT, SESSO, SCUOLA, LSC, ANNO)
mq = lm(PUNT ~ VOTOMAT + LSC + SESSO , data = dati )
summary(mq)
Call:
lm(formula = PUNT ~ VOTOMAT + LSC + SESSO, data = dati)
Residuals:
Min 1Q Median 3Q Max
-13.098 -2.270 0.154 2.367 10.831
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.15403 0.88914 -0.17 0.86
VOTOMAT 0.11658 0.00977 11.93 <2e-16 ***
LSC 2.46286 0.21824 11.28 <2e-16 ***
SESSOM 1.87404 0.21546 8.70 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.53 on 1124 degrees of freedom
Multiple R-squared: 0.206, Adjusted R-squared: 0.203
F-statistic: 96.9 on 3 and 1124 DF, p-value: <2e-16Effetti di associazione: Voto alla maturità (1.1 punto in più ogni 10 voti in più), liceo (2.5 punti in più), genere (Maschi 1.9 punti in più).
\(R^2\) 20%.
Residui
par(mfrow = c(1,2))
plot(mq, which = 1, pch = '.', col = "RoyalBlue")
plot(mq, which = 2, pch = '.')