Lab. 4 - Modelli di regressione multipla
Giovanni Marchetti
null
Scatter in 3D
Consideriamo il data frame derived. Notare il metodo per leggere un data set su un file testo.
setwd("~/Dropbox/Public/gmm/static/R/")
Z = read.table("derived.txt", header = TRUE)
head(Z) Sys Dia Age Hei Wei
1 115 80 35 163 58
2 120 80 31 162 58
3 120 80 22 170 59
4 110 70 23 164 50
5 110 70 18 179 58
6 105 60 28 167 62Costruiamo le variabili DIA (pressione minima), AGE (età) e BMI (massa corporea).
DIA = Z$Dia
AGE = Z$Age
BMI = Z$Wei/(Z$Hei/100)^2Si può ottenere uno scatter in 3D delle variabili ‘BMI’, ‘AGE’ e ‘DIA’ usiamo i comandi
library(scatterplot3d)
s3d = scatterplot3d(BMI, AGE, DIA, pch = 16, highlight.3d = FALSE,
angle = 40, xlab = "BMI", ylab = "Eta'",
zlab = "Pressione")
Tuttavia non è consigliabile perché non molto informativo.
Adattamento di un modello
La funzione è sempre lm. Per esempio per adattare un modello di regressione per la pressione dipendente dalla massa corporea e dall’età si usa l’istruzione
mq = lm(DIA ~ BMI + AGE)
summary(mq)
Call:
lm(formula = DIA ~ BMI + AGE)
Residuals:
Min 1Q Median 3Q Max
-22.86 -6.51 1.18 6.05 13.95
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 56.411 8.869 6.36 1.3e-07 ***
BMI 0.142 0.492 0.29 0.7748
AGE 0.467 0.160 2.92 0.0056 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.29 on 41 degrees of freedom
Multiple R-squared: 0.294, Adjusted R-squared: 0.26
F-statistic: 8.54 on 2 and 41 DF, p-value: 0.000793L’output è analogo a quello del modello di regressione semplice. Il Residual standard error è la stima della deviazione standard degli errori, il Multiple R-squared è l’indice di determinazione lineare \(R^2\), l’adjusted R-squared è l’indice di determinazione aggiustato. L’F-statistic è il test F per l’ipotesi di nullità di tutti i parametri tranne la costante.
Rappresentare il piano di regressione sullo scatter 3D non è semplice:
s3d = scatterplot3d(BMI, AGE, DIA, pch = 16, highlight.3d = FALSE,
angle = 40, xlab = "BMI", ylab = "Eta'",
zlab = "Pressione")
mq = lm(DIA ~ BMI + AGE)
s3d$plane3d(mq, col = "red")
La regressione lineare semplice di DIA da BMI produce in questo caso un coefficiente marginale di regressione che ignora l’età e la stima è molto diversa.
lm(DIA ~ BMI)
Call:
lm(formula = DIA ~ BMI)
Coefficients:
(Intercept) BMI
48.87 1.09 La retta di regressione marginale si può visualizzare sullo scatter 3D:
s3d = scatterplot3d(BMI, AGE, DIA, pch = 16, highlight.3d = FALSE,
angle = 40, xlab = "BMI", ylab = "Eta'",
zlab = "Pressione", type = "n")
s3d$points3d(c(15, 35), c(60, 60), 48.9 + 1.08*c(15, 35), type = "l", lwd = 2, col = "blue")
s3d$points3d(BMI, rep(60, length(AGE)), DIA, pch = 20, col = "red")
Dati Advertising
Regressioni semplici
adv = read.csv("~/Dropbox/Public/gmm/static/R/Advertising.csv")
lm(Sales ~ TV, data = adv)
Call:
lm(formula = Sales ~ TV, data = adv)
Coefficients:
(Intercept) TV
7.0326 0.0475 lm(Sales ~ Radio, data = adv)
Call:
lm(formula = Sales ~ Radio, data = adv)
Coefficients:
(Intercept) Radio
9.312 0.202 lm(Sales ~ Newspaper, data = adv)
Call:
lm(formula = Sales ~ Newspaper, data = adv)
Coefficients:
(Intercept) Newspaper
12.3514 0.0547 Regressione multipla
mod3 = lm(Sales ~ TV + Radio + Newspaper, data = adv)
summary(mod3)
Call:
lm(formula = Sales ~ TV + Radio + Newspaper, data = adv)
Residuals:
Min 1Q Median 3Q Max
-8.828 -0.891 0.242 1.189 2.829
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.93889 0.31191 9.42 <2e-16 ***
TV 0.04576 0.00139 32.81 <2e-16 ***
Radio 0.18853 0.00861 21.89 <2e-16 ***
Newspaper -0.00104 0.00587 -0.18 0.86
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.69 on 196 degrees of freedom
Multiple R-squared: 0.897, Adjusted R-squared: 0.896
F-statistic: 570 on 3 and 196 DF, p-value: <2e-16Correlazioni tra i predittori e con la risposta
round(cor(adv),2) X TV Radio Newspaper Sales
X 1.00 0.02 -0.11 -0.15 -0.05
TV 0.02 1.00 0.05 0.06 0.78
Radio -0.11 0.05 1.00 0.35 0.58
Newspaper -0.15 0.06 0.35 1.00 0.23
Sales -0.05 0.78 0.58 0.23 1.00Matrice del modello
X = model.matrix(mod3)
head(X) (Intercept) TV Radio Newspaper
1 1 230.1 37.8 69.2
2 1 44.5 39.3 45.1
3 1 17.2 45.9 69.3
4 1 151.5 41.3 58.5
5 1 180.8 10.8 58.4
6 1 8.7 48.9 75.0Previsione di vendita
new = data.frame(TV =c(100, 100) , Radio =c(40, 40) , Newspaper = c(50,0))
cbind(new, y = predict(mod3, newdata = new)) TV Radio Newspaper y
1 100 40 50 15.00
2 100 40 0 15.06Errore standard dei residui
df = nrow(X) - 4
s2 = deviance(mod3) / df
RSE = sqrt(s2)Errori standard
S = solve(t(X) %*% X) * s2
sqrt(diag(S))(Intercept) TV Radio Newspaper
0.311908 0.001395 0.008611 0.005871 summary(mod3)$coef Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.938889 0.311908 9.4223 1.267e-17
TV 0.045765 0.001395 32.8086 1.510e-81
Radio 0.188530 0.008611 21.8935 1.505e-54
Newspaper -0.001037 0.005871 -0.1767 8.599e-01Intervalli di confidenza e test
Gli intervalli di confidenza separati per ciascun coefficiente di regressione si ottengono con
confint(mod3) 2.5 % 97.5 %
(Intercept) 2.32376 3.55402
TV 0.04301 0.04852
Radio 0.17155 0.20551
Newspaper -0.01262 0.01054Proprietà dei residui
I residui sono \(\boldsymbol{e} = \boldsymbol y - \boldsymbol X \hat{ \boldsymbol \beta}\).
e = residuals(mod3)
yh = fitted(mod3)
t(X) %*% e [,1]
(Intercept) 7.327e-15
TV -1.023e-12
Radio 1.251e-12
Newspaper 4.441e-13yh %*% e [,1]
[1,] -9.948e-14Esercizi 4
Nel modello di regressione lineare multipla una delle due affermazioni seguenti è sempre vera e l’altra è tipicamente falsa. Quale e perché? \[ \boldsymbol \epsilon \perp \boldsymbol X, \qquad \boldsymbol e \perp \boldsymbol X \]
Nel modello di regressione lineare multipla una delle due affermazioni seguenti è sempre vera e l’altra è tipicamente falsa. Quale e perché? \[ \boldsymbol e \perp \boldsymbol X, \qquad \boldsymbol e \perp\hspace{-4pt}\perp \boldsymbol X \] Nota: il simbolo \(\perp\hspace{-4pt}\perp\) significa indipendenza stocastica.
Guardate le due equazioni \[ \boldsymbol Y = \boldsymbol X\boldsymbol \beta + \boldsymbol \epsilon \qquad \boldsymbol Y = \boldsymbol X \hat{\boldsymbol \beta} + \boldsymbol e \] Qual è il vero modello di regressione? Quali sono gli errori e quali sono i residui? Qual è il parametro e quale lo stimatore?
- Vero o falso e giustifica.
- Nel modello di regressione lineare \(E(\hat{\boldsymbol Y} \mid \boldsymbol X) = \boldsymbol X \hat{\boldsymbol \beta}.\)
- Nel modello di regressione lineare \(E(\hat{\boldsymbol Y} \mid \boldsymbol X) = \boldsymbol X \boldsymbol \beta.\)
- Nel modello di regressione lineare \(E(\boldsymbol Y \mid \boldsymbol X) = \boldsymbol X \boldsymbol \beta.\)
Dimostrare che nel modello di regressione lineare semplice con un predittore espresso in scarti dalla media \[ Y_i = \beta_0 + \beta_1 (x_i -\bar x) + \epsilon_i \] le stime dei minimi quadrati sono \[ \hat{\beta_0} = \bar Y, \quad \hat{\beta_1} = s_{XY} / s_{XX}. \]
Nel modello di regressione multipla \[ Y_i = \beta_0 + \beta_1 (X_{i1} - \bar x_1) + \beta_2 (X_{i2} - \bar x_2) + \epsilon_i \] è vero che se \(X_{i1}\) e \(X_{i2}\) sono incorrelate allora le stime dei minimi quadrati \(\hat\beta_1\) e \(\hat \beta_2\) sono uguali alle stime separate delle regressioni di \(Y\) da \(X_1\) e di \(Y\) da \(X_2\)? Fornire in ogni caso una dimostrazione o un controesempio.
I rendimenti \(Y_i\) di un processo chimico misurato a nove diverse temperature \(t_i\) sono i seguenti.
temp = c(15, 16, 17, 18, 19, 20, 21, 22, 23)
y = c(90.0, 90.9, 90.7, 87.9, 86.4, 82.5, 80.0, 76.0, 70.0)Si consideri il modello polinomiale \[ Y_i =\beta_0 + (t_i − 19)\beta_1 +(y_i −19)^2 β_2 + \epsilon_i,\quad \epsilon_i \sim N(0,σ^2). \]
- Stimare i parametri.
- Determinare l’indice di determinazione lineare del modello di regressione polinomiale e confrontarlo con quello del modello di regressione lineare semplice
- Fare un grafico dei residui per entrambi i modelli.
- Stimare il rendimento atteso quando la temperatura è di 17.5.
- Stimare la temperatura alla quale il rendimento atteso è maggiore. Qual è l’errore standard della stima?
Soluzioni
Falso (perché \(\epsilon\) indipendente da \(X\)) e Vero (i residui dei MQ sono ortogonali alle colonne di \(X\)).
La seconda perché i residui non sono indipendenti da \(X\) ma incorrelati.
Il modello è il primo. Gli errori sono \(\epsilon\) e i residui sono \(e\). Il parametro è \(\beta\) e lo stimatore \(\hat \beta\).
- \(\hat Y = X\hat \beta\). Quindi \(E(\hat Y \mid X) = XE(\hat \beta) = X\beta\). Falso
- Vero
- Vero
- La matrice del modello è \(X = [1, x_i-\bar x]\). Quindi \[
(X^T X)^{-1} = \begin{pmatrix}
n & 0 \\
0 & \sum_i (x_i - \bar x)^2\\
\end{pmatrix}^{-1},
\qquad
\begin{pmatrix}
\sum y_i \\
\sum_i [y_i (x_i - \bar x)]\\
\end{pmatrix}
\] Quindi
\begin{align*}
\hat \beta_0 &= \sum_i y_i/n\\
\hat \beta_1 &= \frac{\sum_i [y_i (x_i - \bar x)]}{\sum_i (x_i - \bar x)^2}
\end{align*}
l risultato segue dall’identità \[ \sum_i (y_i - \bar y)(x_i - \bar x) = \sum_i y_i (x_i - \bar x). \]
- La matrice del modello è \([1, x_{i1} - \bar x_1, x_{12} - \bar x_2]\). Quindi \[ (X^T X) = \begin{pmatrix} n & 0 & 0\\ . & \sum_i (x_{i1} - \bar x_1)^2 & \sum_i (x_{i1} - \bar x_1)(x_{i2} - \bar x_2) \\ . & . & \sum_i (x_{i2} - \bar x_2)^2\\ \end{pmatrix}, \] Perciò se \(X_1\) e \(X_2\) sono incorrelate la matrice è diagonale e si inverte facilmente. Inoltre \[ \begin{pmatrix} \sum_i y_i \\ \sum_i [y_i(x_{i1} - \bar x_1)]\\ \sum_i [y_i (x_{i2} - \bar x_2)]\\ \end{pmatrix} \] e quindi risulta che \(\hat \beta_0 = \bar y\) e \begin{align*} \hat \beta_1 &= \frac{\sum_i [y_i (x_{i1} - \bar x_1)]}{\sum_i (x_{i1} - \bar x_1)^2} \\ \hat \beta_2 &= \frac{\sum_i [y_i (x_{i2} - \bar x_2)]}{\sum_i (x_{i2} - \bar x_2)^2} \\ \end{align*}
Ecco l’analisi col modello della retta
plot(temp, y, pch = 20 , col = "purple", xlab = "Temperatura", ylab = "Rendimento")
ml = lm(y ~ temp)
summary(ml)
Call:
lm(formula = y ~ temp)
Residuals:
Min 1Q Median 3Q Max
-3.922 -0.497 1.203 1.553 2.578
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 131.797 6.269 21.02 1.4e-07 ***
temp -2.525 0.327 -7.72 0.00011 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.53 on 7 degrees of freedom
Multiple R-squared: 0.895, Adjusted R-squared: 0.88
F-statistic: 59.6 on 1 and 7 DF, p-value: 0.000114abline(ml, col = "Green", lwd = 2)
Con il grafico dei residui
plot(ml, which = 1)
Questo è il modello quadratico
plot(temp, y, pch = 20 , col = "purple", xlab = "Temperatura", ylab = "Rendimento")
mq = lm(y ~ I(temp - 19) + I((temp-19)^2))
summary(mq)
Call:
lm(formula = y ~ I(temp - 19) + I((temp - 19)^2))
Residuals:
Min 1Q Median 3Q Max
-0.895 -0.467 0.110 0.367 0.841
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 86.2905 0.3377 255.54 2.4e-13 ***
I(temp - 19) -2.5250 0.0863 -29.27 1.1e-07 ***
I((temp - 19)^2) -0.3702 0.0381 -9.72 6.8e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.668 on 6 degrees of freedom
Multiple R-squared: 0.994, Adjusted R-squared: 0.992
F-statistic: 476 on 2 and 6 DF, p-value: 2.46e-07b = coef(mq)
curve(b[1] + b[2] * (x-19) + b[3]*(x-19)^2, add = TRUE, col = "RoyalBlue", lwd = 2)
abline(v = 15.59, lty = 2)
Con il grafico dei residui (molto migliore)
plot(mq, which = 1)
Il rendimento atteso quando la temperatura è 17.5 è
predict(mq, data.frame(temp = 17.5)) 1
89.24 La temperatura \(\theta\) alla quale il rendimento è massimo è \[ \theta - 19 = - \beta_1/(2 \beta_1) \]
cioè si stima con
19 - b[2]/(2*b[3])I(temp - 19)
15.59 Vedi la figura.
L’errore standard della stima non è banale da ottenere perché bisogna calcolare \[ \mathrm{var}(\hat \theta) = \mathrm{var}(19 - \hat \beta_1/(2 \hat \beta_1)) = \mathrm{var}(\hat \beta_1/(2 \hat \beta_1)) \] e si può fare in modo approssimato col metodo delta.