Transcript plot(modello)

LABORATORIO DI STATISTICA
AZIENDALE
3. Modello di Regressione lineare con R
Enrico Properzi - [email protected]
A.A. 2010/2011
L'analisi della regressione è utilizzata per spiegare la
relazione esistente tra una variabile Y (continua), detta
variabile risposta, e una o più variabili esplicative dette
anche
covaraite,
predittori,
regressori
o
variabili
indipendenti.
In termini di funzione abbiamo:
Alla prima componente, che rappresenta la parte della
variabile risposta spiegata dai predittori, si deve
aggiungere unna seconda componente, detta casuale, che
rappresenta quella parte di variabilità della variabile
dipendente che non può essere ricondotta a fattori
sistematici oppure facilmente individuabili, ma dovuti al
caso.
Il legame funzionale può essere di qualunque tipo, e quando
si utilizza una funzione di tipo lineare si parla di
regressione lineare multipla o modello lineare che assume
la seguente formulazione:
dove:
- β0 è l'intercetta o termine nullo
- β1, … , βk sono i coefficienti di regressione delle variabili
esplicative
Questi parametri, insieme alla varianza dell'errore, sono gli
elementi del modello da stimare sulla base delle osservazioni
campionarie.
La stima dei parametri di un modello di regressione multipla con il
software R viene eettuata mediante il comando lm():
lm(formula, data, ....)
dove:
formula è il modello che deve essere analizzato e viene espresso
nella forma:
Y = X1 + X2 + …
data è il data frame che contiene le variabili da analizzare
E' possibile creare un oggetto di classe lm contenente i risultati
e i parametri
del modello stimato
Esempio:
>
>
>
>
>
>
>
peso <- c(80, 68, 72, 75, 70, 65, 62, 60, 85, 90)
altezza <- c(175, 168, 170, 171, 169, 165, 165, 160, 180, 186)
modello <- lm(peso ~ altezza)
plot(altezza, peso)
abline (modello, col="blue")
segments(altezza, fitted(modello), altezza, peso, lty=2)
title (main="Plot della regressione peso altezza")
> summary (modello)
Call:
lm(formula = peso ~ altezza)
Residuals:
Min
1Q
-3.2293 -0.8445
Median
0.1051
3Q
1.0207
Max
2.1734
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -143.69578
13.43899 -10.69 5.14e-06 ***
altezza
1.26621
0.07857
16.12 2.21e-07 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.807 on 8 degrees of freedom
Multiple R-squared: 0.9701,
Adjusted R-squared: 0.9664
F-statistic: 259.7 on 1 and 8 DF, p-value: 2.206e-07
Dall’oggetto di classe lm è possibile estrarre i singoli valori
avvalendoci di funzioni specifiche:






coef ()
residuals ()
fitted()
deviance()
formula()
anova()
coefficienti di regressione
residui del modello
valori risposta stimati dal modello
devianza dei residui
formula del modello
produce la tavola della varianza
associata alla regressione
> anova(modello)
Analysis of Variance Table
Response: peso
Df Sum Sq Mean Sq F value
Pr(>F)
altezza
1 847.98 847.98 259.75 2.206e-07 ***
Residuals 8 26.12
3.26
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Nell’analizzare i risultati di una regressione deve essere
verificato se le ipotesi sottostanti al modello dei minimi
quadrati ordinari sono soddisfatte
Per fare questo si può far ricorso ad alcuni grafici:
1.
Il grafico dei
valutare se le
incorrelazione
questo grafico

plot(fitted(modello), residuals(modello))
abline (0,0)

residui contro i valori stimati permette di
ipotesi di omoschedasticità, media nulla e
dei residui sono verificate. In un buon modello
dovrebbe apparire completamente casuale
Il Q-Q plot permette di verificare se i residui sono distribuiti
normalmente
qqnorm(residuals(modello))
qqline(residuals(modello), col="red")
Quanto più i punti che
rappresentano i residui
ordinati giacciono in
prossimità della linea Q-Q
tanto più è plausibile
l’assunzione di normalità
Il comando plot(modello) fornisce quattro grafici per la
valutazione delle ipotesi sottostanti al modello;
1.
2.
3.
4.
Il grafico dei residuoi rispetto ai valori previsti consente di
valutare se ci sono delle tendenze nella distribuzione dei
residui stessi oppure una variabilità costante
Il Q-Q plot dei residui standardizzati permette di verificare
se gli errori sono distribuiti normalmente
Lo scale-location plot è utilizzato per determinare se la
distribuzione dei residui è costante su tutto il range dei
valori previsti ed è utile nell’individuazione di valori
outlier
Il grafico residuals vs leverage (residui standardizzati
rispetto ai punti di leva) ha lo scopo di evidenziare se ci
sono dati anomali che possono influenzare la stima del modello.
Per verificare che siano valide le ipotesi di base è possibile
fare ricorso ad alcuni test specifici:
Ad esempio, per verificare che la media degli errori non sia
significativamente diversa da zero si utilizza il test t di
Student:
> t.test(resid(modello))
One Sample t-test
data: resid(modello)
t = 0, df = 9, p-value = 1
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-1.218611 1.218611
sample estimates:
mean of x
6.102974e-17
Il file prezzicase.txt contiene alcune informazioni su un un
gruppo di immobili in vendita
> case <- read.table("prezzicase.txt", header=T)
> str(case)
'data.frame':
117 obs. of 7 variables:
$ Price
: int 2050 2080 2150 2150 1999 1900 1800 1560 1450 1449 ...
$ SQFT
: int 2650 2600 2664 2921 2580 2580 2774 1920 2150 1710 ...
$ Age
: Factor w/ 31 levels "*","1","13","14",..: 3 1 28 17 20 20 10 2
1 2 ...
$ Features: int 7 4 5 6 4 4 4 5 4 3 ...
$ NE
: int 1 1 1 1 1 1 1 1 1 1 ...
$ Corner : int 0 0 0 0 0 0 0 0 0 0 ...
$ Tax
: Factor w/ 96 levels "*","1010","1035",..: 20 5 12 19 21 18 22
9 1 2 ...
Si vuole studiare un modello di regressione in cui il prezzo
delle case è funzione della superficie (SQFT), dell’età
dell’immobile (Age) e del numero di caratteristiche presenti
tra un gruppo di 11 (Features: lavastoviglie, frigorifero,
forno a microonde, ecc..)
Un primo step consiste nell’analisi delle relazioni che
sussistono tra le variabili utilizzate mediante la valutazione
della loro matrice di correlazione:
> cor (case[,c("Price", "SQFT", "Features", "Age")])
Price
SQFT
Features
Age
Price
1.0000000 0.8447951 0.4202725 -0.1956308
SQFT
0.8447951 1.0000000 0.3949250 -0.1515718
Features 0.4202725 0.3949250 1.0000000 -0.2069562
Age
-0.1956308 -0.1515718 -0.2069562 1.0000000
Successivamente stimiamo
il modello di regressione:
> lmcase <- lm(Price~SQFT+Features+Age)
> summary(lmcase)
Call:
lm(formula = Price ~ SQFT + Features + Age)
Residuals:
Min
1Q
-1008.623 -103.122
Median
-5.155
3Q
73.177
Max
796.825
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 41.88434
78.13318
0.536
0.5930
SQFT
0.58090
0.03909 14.860
<2e-16 ***
Features
25.18350
14.71659
1.711
0.0898 .
Age
-1.73798
1.59204 -1.092
0.2773
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 202 on 113 degrees of freedom
Multiple R-squared: 0.7255,
Adjusted R-squared: 0.7182
F-statistic: 99.54 on 3 and 113 DF, p-value: < 2.2e-16
BONTA’ DI ADATTAMENTO:
Un primo aspetto utile per la valutazione del modello di
regressione è la valutazione della bontà di adattamento del
modello (goodness of fit).
Le statistiche maggiormente usate a questo scopo sono l’errore
standard della stima e l’R2
L’errore standard della stima (residual standard error)
corrisponde all’errore standard dei residui e rappresenta un
indice che esprime l’ampiezza dell’errore di misura del modello
considerato.
R2 esprime la parte di varianza spiegata attraverso il modello.
Varia sempre tra 0 e 1 e può essere interpretato come la
percentuale di varianza della variabile dipendente spiegata
dalle variabili indipendenti utilizzate nel modello.
VERIFICA CAPACITA’ PREDITTIVA DEL MODELLO
Per verificare se la previsione della variabile dipendente Y
migliora significativamente mediante il modello di regressione
si pone a confronto la varianza spiegata dal modello con la
varianza residua (non spiegata).
L’ipotesi H0 che si sottopone a verifica è che la varianza
spiegata sia uguale alla varianza residua, cioè che il modello
non migliora l’errore di rpevisione della variabile dipendente.
Per la verifica dell’ipotesi si usa il test F (rapporto tra le
varianze) che si distribuisce come una variabile casuale F di
Fisher
Si può anche effettuare l’analisi della varianza su tutte le
variabili inserite nel modello
> anova(lmcase)
Analysis of Variance Table
Response: Price
Df
Sum Sq Mean Sq F value Pr(>F)
SQFT
1 11981915 11981915 293.7570 < 2e-16 ***
Features
1
149320
149320
3.6608 0.05824 .
Age
1
48609
48609
1.1917 0.27730
Residuals 113 4609103
40789
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
CONTRIBUTO DEI SINGOLI PREDITTORI
Se abbiamo scartato H0 possiamo approfondire l’analisi indagando
il contributo di ciascun predittore considerato singolarmente.
Si testa l’ipotesi che ciascun coefficiente stimato sia uguale a
0. Si utilizza per questo la statistica t:
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 41.88434
78.13318
0.536
0.5930
SQFT
0.58090
0.03909 14.860
<2e-16 ***
Features
25.18350
14.71659
1.711
0.0898 .
Age
-1.73798
1.59204 -1.092
0.2773
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
VERIFICA ASSUNZIONI OLS
Una prima valutazione del rispetto delle ipotesi sottostanti al
modello lineare si effettua utilizzando i grafici standard
prodotti dal comando plot(modello)
> par (mfrow=c(2,2))
> plot(lmcase)
La media dei residui, inoltre, non è significativamente diversa
da zero:
> t.test(resid(lmcase))
One Sample t-test
data: resid(lmcase)
t = 0, df = 116, p-value = 1
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-36.49965 36.49965
sample estimates:
mean of x
-3.491835e-15
Quale modello finale???