Lab. 3 - Modelli di regressione
Giovanni Marchetti
Normale multivariata
Generiamo uan normale tripla di media 0 e matrice di covarianza \[ \Sigma = \begin{pmatrix} 1 & 0.2 & 0.8 \\ & 1 & 0.16 \\ & & 1 \\ \end{pmatrix} \]
Metodo 1:
Trasformazione lineare. Prima definisco la matrice di covarianza
S = matrix(c(1, 0.2, 0.8, 0.2, 1, 0.16, 0.8, 0.16, 1), 3, 3)
S [,1] [,2] [,3]
[1,] 1.0 0.20 0.80
[2,] 0.2 1.00 0.16
[3,] 0.8 0.16 1.00Poi calcolo la radice \(B = U \Lambda^{1/2} U^T\) tale che \(\Sigma = B B\)
e = eigen(S)
U = e$vectors
B = U %*% diag(sqrt(e$values)) %*% t(U)Notare che la moltiplicazione di matrici è %*%, non *. Quindi trasformiamo un vettore di variabili normali standard indipendenti.
Z = matrix(rnorm(3), 3, 1)
X = B %*% Z
X [,1]
[1,] 1.9953
[2,] -0.5205
[3,] 2.0147Per generare più osservazioni basta trasformare con \(B\) una matrice \(Z\).
Z = matrix(rnorm(300), 3, 100)
X = B %*% Z
X = t(X)
pairs(data.frame(X), col = "darkgreen", pch = 19)
cov(X) [,1] [,2] [,3]
[1,] 0.9347 0.1001 0.8165
[2,] 0.1001 0.8895 0.1998
[3,] 0.8165 0.1998 1.1395Metodo 2
Usando una funzione apposita.
library(MASS)
X =mvrnorm(100, mu = c(0,0,0), Sigma = S)
pairs(data.frame(X), col = "purple", pch = 19)
cov(X) [,1] [,2] [,3]
[1,] 0.80678 0.09845 0.64850
[2,] 0.09845 0.77045 0.06271
[3,] 0.64850 0.06271 0.80041Calcolo del valore atteso condizionato
Se voglio calcolare \(E(X_1\mid x_2, x_3)\) faccio così:
a = 1; b = c(2, 3)
rownames(S) = c(1, 2, 3)
colnames(S) = c(1, 2,3)
S[a,a][1] 1S[a,b] 2 3
0.2 0.8 S[b,b] 2 3
2 1.00 0.16
3 0.16 1.00Quindi calcolo i coefficienti \(B = \Sigma_{ab}\Sigma_{bb}^{-1}\):
B = S[a,b] %*% solve(S[b,b]) # solve calcola l'inversa
B 2 3
[1,] 0.07389 0.7882Quindi \(E(X_1\mid x_2, x_3) = 0.07389 x_2 + 0.7882 x_3\).
Funzione di regressione
Supponiamo che nella realtà i dati siano generati come segue \[ Y_i = \exp[(X_i - 0.3)^2] - 1 + \epsilon_i = r(X_i) + \epsilon_i \] dove \(x_i = 0, 0.02, \dots, 1\) e \(\epsilon_i \sim N(0, 0.02)\). Generiamo i dati.
set.seed(1231)
x = seq(0, by = 0.02, 1)
e = rnorm(length(x), 0, sd = sqrt(0.02))
y = exp((x - 0.3)^2) - 1 + e
plot(x,y, pch = 20, col = "RoyalBlue")
curve(exp((x - 0.3)^2) - 1, 0, 1, add = TRUE, col = "gray", lwd = 2)
Poiché non conosciamo il meccanismo generatore specifichiamo un modello, il modello di regressione lineare semplice \[ Y_i = \beta_0 + \beta_1 x_i + \epsilon_i \]
Adattiamo il modello con il metodo dei minimi quadrati.
plot(x,y, pch = 20, col = "RoyalBlue")
mq = lm(y ~ x)
abline(mq, col = "red", lwd = 2)
summary(mq)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-0.3196 -0.1244 0.0027 0.0841 0.4568
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0675 0.0449 -1.50 0.14
x 0.4676 0.0774 6.04 2e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.163 on 49 degrees of freedom
Multiple R-squared: 0.427, Adjusted R-squared: 0.415
F-statistic: 36.5 on 1 and 49 DF, p-value: 2.02e-07Vediamo i residui e i valori adattati
out = cbind(y, fitted(mq), residuals(mq))
colnames(out) = c("osservati", "adattati", "residui")
round(out, 3) osservati adattati residui
1 0.083 -0.067 0.150
2 0.180 -0.058 0.238
3 -0.179 -0.049 -0.130
4 0.110 -0.039 0.150
5 0.344 -0.030 0.374
6 -0.009 -0.021 0.012
7 -0.031 -0.011 -0.020
8 0.215 -0.002 0.217
9 0.010 0.007 0.003
10 0.102 0.017 0.085
11 0.149 0.026 0.123
12 -0.054 0.035 -0.090
13 0.161 0.045 0.116
14 -0.129 0.054 -0.183
15 0.147 0.063 0.083
16 0.050 0.073 -0.023
17 -0.057 0.082 -0.139
18 -0.228 0.092 -0.320
19 -0.018 0.101 -0.119
20 -0.033 0.110 -0.143
21 0.068 0.120 -0.052
22 0.192 0.129 0.063
23 0.219 0.138 0.081
24 0.081 0.148 -0.067
25 -0.078 0.157 -0.235
26 -0.019 0.166 -0.185
27 0.147 0.176 -0.029
28 -0.105 0.185 -0.290
29 0.199 0.194 0.005
30 0.013 0.204 -0.191
31 -0.012 0.213 -0.225
32 0.025 0.222 -0.197
33 0.271 0.232 0.039
34 0.231 0.241 -0.010
35 0.141 0.251 -0.110
36 0.310 0.260 0.050
37 0.312 0.269 0.043
38 0.460 0.279 0.182
39 0.156 0.288 -0.132
40 0.330 0.297 0.033
41 0.266 0.307 -0.041
42 0.246 0.316 -0.070
43 0.142 0.325 -0.184
44 0.375 0.335 0.040
45 0.363 0.344 0.019
46 0.354 0.353 0.001
47 0.820 0.363 0.457
48 0.530 0.372 0.158
49 0.517 0.381 0.136
50 0.642 0.391 0.251
51 0.475 0.400 0.075Differenze tra errori e residui
plot(x, y - (exp((x - 0.3)^2) - 1))
plot(x, residuals(mq))
Adattamento di una funzione quadratica
Introduciamo una variabile \(X^2\).
x2 = x*x
mq2 = lm(y ~ x + x2)
b = coef(mq2)
b(Intercept) x x2
0.1322 -0.7548 1.2225 plot(x,y, pch = 20, col = "RoyalBlue")
curve(b[1] + b[2] * x + b[3]*x^2, 0, 1, add= TRUE, col = "red", lwd = 2)
Osservare che non è la “vera” funzione del meccanismo generatore.
Regressione non parametrica
scatter.smooth(x, y, pch = 20, col = "RoyalBlue", lpars = c(col = "blue", lwd = 2))
Dati sull’effetto del tempo di gestazione
Lettura dei dati
data("abdom", package = "gamlss.data")Scatter
with(abdom, plot(x, y))
Adattamento del modello di regressione lineare.
mq = lm(y ~ x, data = abdom)
summary(mq)
Call:
lm(formula = y ~ x, data = abdom)
Residuals:
Min 1Q Median 3Q Max
-84.58 -8.11 -0.18 8.06 54.32
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -55.1795 2.0010 -27.6 <2e-16 ***
x 10.3382 0.0701 147.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.6 on 608 degrees of freedom
Multiple R-squared: 0.973, Adjusted R-squared: 0.973
F-statistic: 2.18e+04 on 1 and 608 DF, p-value: <2e-16with(abdom, plot(x, y, col = "gray"))
abline(mq, lwd = 2, col = "purple")
Residui
plot(mq, which = 1, lwd = 2, col = "gray")
Devianza residua
RSS = sum((residuals(mq))^2)
RSS[1] 130135Si può anche calcolare con
deviance(mq)[1] 130135La correlazione tra le due variabili indica la forza della relazione
with(abdom, cor(x,y))[1] 0.9863Previsione della circonferenza dell’addome di un feto di 30 settimane
predict(mq, newdata = data.frame(x = 30)) 1
255 Trasformazioni
Notare che il modello seguente è ancora un modello di regressione lineare semplice tra la risposta \(Y\) e il predittore \(\log(X)\). \[ Y_i = \beta_0 + \beta_1 \log x_i + \epsilon_i \]
Esercizi 3
Studiate il modello di regressione lineare senza intercetta \[ Y_i = \beta X_i + \epsilon_i \] dove \(\epsilon_i\) soddifa le assunzioni usuali. Trovare la stima dei minimi quadrati di \(\beta\). Dimostrate che i valori stimati del modello di regressione lineare senza intercetta si possono scrivere \[ \hat y_i = \sum_{j=1}^n a_j y_j \] e determinate i coefficienti \(a_j\).
Dimostrare che la somma dei residui del modello di regressione lineare semplice con intercetta è sempre zero. La stessa proprietà valre per i residui del modello di regressione senza intercetta? Giustificare.
Caricate il data frame
Autodal packageISLR. Usate la funzionelm()per adattare un modello di regressione lineare semplice alle variabilimpg(risposta) ehorsepower. Studiate prima l’help del data set.
- C’è una relazione tra predittore e risposta?
- Quanto è forte?
- Interpretate i coefficienti
- Prevedere la percorrenza di un auto con 89 cavalli.
- Usando il grafico dei residui verificare le assunzioni
- Caricate il data frame
Animalsdal packageMASS. Contiene due variabilibodyebrainche sono il peso del corpo e del cervello di alcune specie di animali.
- Adattate un modello di regressione lineare semplice prendendo come risposta il peso del cervello. Commentare.
- Adattate un modello di regressione lineare semplice prendendo come risposta il
log(brain)e come variabile esplicattivalog(body). Commentare.
Verificare che se si ha una matrice a blocchi \[ \Sigma = \begin{pmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \\ \end{pmatrix} \] con una inversa \(D = \Sigma^{-1}\) anch’essa a blocchi \[ D = \begin{pmatrix} D_{11} & D_{12} \\ D_{21} & D_{22} \\ \end{pmatrix} \] allora, sotto opportune condizioni, \(D_{11}^{-1} = \Sigma_{11.2}\).
Considerate la nomale multivariata \(X \sim N(\mu. \Sigma)\) con
\[ \mu = \begin{pmatrix} 10\\ 5\\ 1 \end{pmatrix}, \qquad \Sigma = \begin{pmatrix} 1 & 0.2 & 0.8 \\ & 1 & 0.16 \\ & & 1 \\ \end{pmatrix}. \]
- Trovare \(E(X_1 \mid x_2, x_3)\) e \(\mathrm{var}(X_1 \mid x_2, x_3)\).
- Trovare \(E(X_2 \mid x_1, x_3)\) e \(\mathrm{var}(X_2 \mid x_1, x_3)\).
- Trovare \(E(X_2 \mid x_1)\) e \(\mathrm{var}(X_2 \mid x_1)\).
Soluzioni
Minimizza \[ \sum_i (y_i - \beta x_i)^2. \] Equazioni normali \[ (-2) \sum_i [(y_i - \beta x_i)x_i] = 0 \] da cui \(\hat \beta = \frac{\sum_j x_j y_j}{\sum_j x_i^2}.\) Quindi \[ \hat y_i = x_i \frac{\sum_j x_j y_j}{\sum_j x_j^2} = \sum_j a_j y_j \] dove \(a_j = x_i x_j/\sum_j x_j^2\).
Il sistema di equzioni normali soddisfa \[ X^T X\hat \beta = X^T y \] cioè \[ X^T (y - X \hat \beta) = 0 \] Dunque sviluppando \[ \sum_i (y_i - \hat y_i) = 0, \qquad \sum_i [(y_i - \hat y_i)x_i] = 0, \] e la somma dei residui è zero. Al contrario se manca l’intercetta c’è una sola equazione \[ \sum_i [(y_i - \hat y_i) x_i ] = 0 \] e la somma dei residui non è necessariamente zero.
Calcoli
library(ISLR)
data(Auto)
mo = lm(mpg ~ horsepower, data = Auto)
summary(mo)
Call:
lm(formula = mpg ~ horsepower, data = Auto)
Residuals:
Min 1Q Median 3Q Max
-13.571 -3.259 -0.344 2.763 16.924
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.93586 0.71750 55.7 <2e-16 ***
horsepower -0.15784 0.00645 -24.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.91 on 390 degrees of freedom
Multiple R-squared: 0.606, Adjusted R-squared: 0.605
F-statistic: 600 on 1 and 390 DF, p-value: <2e-16All’aumentare dei cavalli la percorrenza diminuisce.
- C’è una relazione tra predittore e risposta? Sì
Quanto è forte? in termini di \(R^2\) è \(0.61\). Il coefficiente di correlazione è \(-0.78\). Quindi abbastanza elevato.
Interpretate i coefficienti: per ogni cavallo in più la percorrenza per miglio diminuisce di 0.16 galloni.
Prevedere la percorrenza di un auto con 89 cavalli.
predict(mo, data.frame(horsepower = 89)) 1
25.89 - Usando il grafico dei residui verificare le assunzioni
plot(mo, which = 1)
Le assunzioni di linearità e omoschedasticità non sono verificate.
- Istruzioni
library(MASS)
data(Animals)
plot(Animals$body, Animals$brain, pch = 20, xlab = "Peso corpo", ylab = "Peso cervello")
abline(lm(brain~body, data = Animals))
I dati coprono un range molto ampio e l’assunzione di linearità non è soddifatta.
Trasformiamo le due variabili in modo logaritmico.
plot(Animals$body, Animals$brain, pch = '.', xlab = "Peso corpo", ylab = "Peso cervello", log = "xy")
text(Animals$body,Animals$brain, rownames(Animals), cex = 0.5)
abline(lm(log(brain, 10)~log(body, 10), data = Animals))
È la formula dell’inversa di una matrice a blocchi.
Abbiamo
R = matrix(c(1, 0.2, 0.8,
0.2, 1, 0.16,
0.8, 0.16, 1), 3, 3)
b1 = R[1,2:3, drop = FALSE] %*% solve(R[2:3, 2:3])
b1 [,1] [,2]
[1,] 0.07389 0.7882b2 = R[2,c(1,3), drop = FALSE] %*% solve(R[c(1,3), c(1,3)])
b2 [,1] [,2]
[1,] 0.2 0v11 = R[1,1] - R[1,2:3] %*% solve(R[2:3, 2:3]) %*% R[2:3, 1]; v11 [,1]
[1,] 0.3547v22 = R[2,2] - R[2,c(1,3)] %*% solve(R[c(1,3), c(1,3)]) %*% R[c(1,3), 2]; v22 [,1]
[1,] 0.96R[1:2, 1:2] [,1] [,2]
[1,] 1.0 0.2
[2,] 0.2 1.0Quindi \[ E(X_1\mid x_2, x_3) = 10 + 0.07389163 (x_2 - 5) + 0.7881773(x_3 - 1), \quad V(X_1\mid x_2, x_3) \simeq 0.355 \] \[ E(X_2\mid x_1, x_3) = 5 + 0.2 (x_1- 10), \quad V(X_2\mid x_1, x_3) = 0.96 \] \[ E(X_2 \mid x_1) = 5 + 0.2 (x_1 - 10), \quad V(X_2\mid x_1) = 1 - 0.2^2 = 0.96 \]