LABORATORIO DI STATISTICA AZIENDALE

1. 3. Modello di Regressione lineare con R LABORATORIO DI STATISTICA AZIENDALE Enrico Properzi - enrico.properzi3@unibo.it A.A. 2010/2011

2. 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:

3. 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

4. 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")

5. > summary (modello) Call: lm(formula = peso ~ altezza) Residuals: Min 1Q Median 3Q Max -3.2293 -0.8445 0.1051 1.0207 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 degreesoffreedom Multiple R-squared: 0.9701, AdjustedR-squared: 0.9664 F-statistic: 259.7 on 1 and 8 DF, p-value: 2.206e-07

6. Dall’oggetto di classe lm è possibile estrarre i singoli valori avvalendoci di funzioni specifiche: coef () coefficienti di regressione residuals () residui del modello fitted() valori risposta stimati dal modello deviance() devianza dei residui formula() formula del modello anova() 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

7. 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: • Il grafico dei residui contro i valori stimati permette di valutare se le ipotesi di omoschedasticità, media nulla e incorrelazione dei residui sono verificate. In un buon modello questo grafico dovrebbe apparire completamente casuale • plot(fitted(modello), residuals(modello)) • abline (0,0)

9. 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à

10. Il comando plot(modello) fornisce quattro grafici per la valutazione delle ipotesi sottostanti al modello; 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.

11. 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

12. 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..)

13. 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)

14. > summary(lmcase) Call: lm(formula = Price ~ SQFT + Features + Age) Residuals: Min 1Q Median 3Q Max -1008.623 -103.122 -5.155 73.177 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

15. 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.

16. 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

17. 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

18. 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

19. 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)

20. 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

21. Quale modello finale???