4. Modello di Regressione lineare con R parte - 2 Enrico Properzi -...

4. Modello di Regressione lineare con Rparte - 2

Enrico Properzi - [email protected] A.A. 2010/2011

Il file azfrutt.txt contiene alcuni dati tratti dai bilanci di un campione di aziende agricole:

codaz Codice dell’azienda ote orientamento produttivo ricavi produzione lorda vendibile capimp capitale impiegato ulut unità lavorative totali sau superficie agricola

frutta <- read.table("azfrutt.txt", header=T)

str(frutta)'data.frame': 36 obs. of 6 variables: $ cod : int 1 2 3 4 5 6 7 8 9 10 ... $ ote : int 31 31 32 31 32 32 31 32 31 31 ... $ ricavi: int 19688 48313 39915 50441 79577 84774 15499 46565 31121 ... $ capimp: int 169459 556788 206464 130666 454788 658214 123123 474778... $ ulut : num 0.62 1.44 2 0.8 2.8 0.7 0.44 1.32 2.1 0.9 ... $ sau : num 6.60 24.53 5.50 7.89 17.12 ...

La variabile ote riporta le tipologie di coltura prevalenti31 viticoltura32 frutta

Dobbiamo trasformarle in fattori:

frutta$ote <- factor(frutta$ote, labels=c("viticoltura", "frutta"))

Studiamo inizialmente un modello di regressione in cui la produzione sia funzione della superficie agricola (sau), del lavoro (ulut) e del capitale (capimp).

ricavi = β0 + β1sau + β2ulut + β3capimp + ε

> m <- lm(ricavi ~ sau + ulut + capimp)> summary (m)

Call:lm(formula = ricavi ~ sau + ulut + capimp)

Residuals: Min 1Q Median 3Q Max -29662 -7327 1290 8016 21822

Stime dei coefficienti e test per la verifica di ipotesi che siano significativamente diversi da 0:

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.393e+04 5.920e+03 2.352 0.0250 * sau 1.907e+03 3.171e+02 6.014 1.04e-06 ***ulut 2.500e+03 3.078e+03 0.812 0.4227 capimp 2.456e-02 9.175e-03 2.677 0.0116 * ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Indicatori per la valutazione della bontà di adattamento del modello:

Residual standard error: 11410 on 32 degrees of freedomMultiple R-squared: 0.7529, Adjusted R-squared: 0.7298 F-statistic: 32.51 on 3 and 32 DF, p-value: 7.803e-10

Per poter confrontare i coefficienti di regressione parziale tra loro , è necessario calcolare i coefficienti di regressione standardizzati:

coef(m)[2]*(sd(sau)/sd(ricavi)) sau 0.6595706

coef(m)[3]*(sd(ulut)/sd(ricavi)) ulut 0.07225166

coef(m)[4]*(sd(capimp)/sd(ricavi)) capimp 0.2963803

Verifichiamo le assunzioni del modello OLS:

par(mfrow=c(2,2)) plot(m) par(mfrow=c(1,1))

Tra le variabili indipendenti della formula del comando lm è possibile inserire anche variabili di tipo factor

Utilizzando i dati dell’esempio precedente(aziende agricole) possiamo introdurre un nuovo regressore di tipo qualitativo: ote

Consideriamo il modello:

ricavi = β0 + β1sau + β2ote + β3capimp + ε

in cui la variabile fattore ote può assumere due modalità: “viticoltura” o “frutta”.

Nell’utilizzare questo fattore, R considera la variabile dummy con valore0 quando ote = “viticoltura”1 quando ote = “frutta”

m1 <- lm(ricavi ~ sau + ote + capimp)> summary(m1)

Call:lm(formula = ricavi ~ sau + ote + capimp)

Residuals: Min 1Q Median 3Q Max -22025 -7541 1021 5704 20822

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.660e+04 4.113e+03 4.036 0.000317 ***sau 1.680e+03 3.207e+02 5.239 9.91e-06 ***otefrutta 8.663e+03 4.084e+03 2.121 0.041737 * capimp 2.249e-02 8.585e-03 2.620 0.013349 * ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10790 on 32 degrees of freedomMultiple R-squared: 0.7789, Adjusted R-squared: 0.7582 F-statistic: 37.58 on 3 and 32 DF, p-value: 1.338e-10

Volendo interpretare il significato dei coefficienti possiamo dire che:

Ogni ettaro di superificie agricola in più determina un incremento di 1.680 euro di ricavi. Questo è valido sia per le aziende con vite sia per quelle con frutta.

Ogni euro di capitale investito determina un incremento di 0.0225 euro di ricavi. Questo è valido sia per le aziende con vite sia per quelle con frutta.

A parità degli altri fattori, le aziende frutticole tendono ad avere ricavi di 8.663 euro più alti rispetto a quelle viticole.

ricavi = β0 + β1sau + β3capimp se vinicola

ricavi = (β0 + β2) + β1sau + β3capimp se frutticola

>plot(sau, ricavi, col=c("blue", "red")[ote]) > plot(capimp, ricavi, col=c("blue", "red")[ote])

Il problema che spesso si deve affrontare nella costruzione di un modello di regressione riguarda la scelta delle variabili esplicative.

Nel modello bisognerebbe includere quelle variabili esplicative la cui varaizione apporta un contributo reale alla variazione della variabile risposta. In genere aumentando il numero di regressori inseriti nel modello la devianza dei residui tende a diminuire.

Tuttavia alcune variabili esplicative potrebbero risultare statisticamente significative, e quindi venire incluse nel modello, unicamente per fattori dovuti al caso. AL contrario variabili esplicative logicamente fondamentali potrebbero risultare statisticamente non significative ed essere così escluse dal modello.

La strategia complessiva della scelta di varaibili si può articolare nelle seguenti fasi:

Decidere quali sono le variabili che costituiscono l’insieme più ampio dei k regressori

Trovare uno o più sottoinsiemi di variabili (p) che spiegano bene la variabile risposta

Applicare una regola di arresto per decidere quante variabili esplicative usare

Stimare i coefficienti di regressione Saggiare la bontà del modello ottenuto

SELEZIONE DI VARIABILI E AGGIORNAMENTO DEL MODELLO

Il problema della scelta della regola di arresto viene risolto con l’Akaike Information Criterion (AIC) o il Bayes Information Criterion (BIC) entrambi basati sui logaritmi della verosimiglianza

AIC = -2*logL + k*df

dove: L è la verosimiglianza df sono i gradi di libertàk i parametri

Il comando AIC(modello) fornisce il valore di AIC per il modello stimato

Il comando AIC(modello, k=log(n)), dove n = numero di osservazioni, fornisce il valore di BIC del modello.

Consideriamo il dataset risparmi.txt relativo ai dati economici di 50 paesi.

dpi è il reddito disponibile pro-capite in dollari USA, ddpi è il tasso percentuale di variazione del reddito disponibile pro

capite, RS è il risparmio personale aggregato diviso per reddito disponibile. pop15 è la percentuale di popolazione sotto i 15 anni pop75 è la percentuale di popolazione oltre i 75 anni.

risparmi <- read.table("risparmi.txt", header=T) str(risparmi)'data.frame': 50 obs. of 6 variables: $ Nazione: Factor w/ 50 levels "Australia","Austria",..: 1 2 3 4 5 6 7 ... $ sr : num 11.43 12.07 13.17 5.75 12.88 ... $ pop15 : num 29.4 23.3 23.8 41.9 42.2 ... $ pop75 : num 2.87 4.41 4.43 1.67 0.83 2.85 1.34 0.67 1.06 1.14 ... $ dpi : num 2330 1508 2108 189 728 ... $ ddpi : num 2.87 3.93 3.82 0.22 4.56 2.43 2.67 6.51 3.08 2.8 ...

Innanzitutto consideriamo un modello con tutti i predittori:

g <- lm(sr ~ pop15 + pop75 + dpi + ddpi)

summary(g)

Call:lm(formula = sr ~ pop15 + pop75 + dpi + ddpi)

Residuals:14 Min 1Q Median 3Q Max -8.2422 -2.6857 -0.2488 2.4280 9.7509

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 28.5660865 7.3545161 3.884 0.000334 ***pop15 -0.4611931 0.1446422 -3.189 0.002603 ** pop75 -1.6914977 1.0835989 -1.561 0.125530 dpi -0.0003369 0.0009311 -0.362 0.719173 ddpi 0.4096949 0.1961971 2.088 0.042471 * ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.803 on 45 degrees of freedomMultiple R-squared: 0.3385, Adjusted R-squared: 0.2797 F-statistic: 5.756 on 4 and 45 DF, p-value: 0.0007904

Notiamo che alcuni regressori risultano poco significativi, in particolare dpi.

Potremmo stimare un secondo moello escludendo questa variabile dai regressori:

g2 <- lm(sr ~ pop15 + pop75 + ddpi)summary(g2)

Call:lm(formula = sr ~ pop15 + pop75 + ddpi)

Residuals: Min 1Q Median 3Q Max -8.2539 -2.6159 -0.3913 2.3344 9.7070

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 28.1247 7.1838 3.915 0.000297 ***pop15 -0.4518 0.1409 -3.206 0.002452 ** pop75 -1.8354 0.9984 -1.838 0.072473 . ddpi 0.4278 0.1879 2.277 0.027478 * ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Potrei eliminare anche il regressore pop75, che non risulta significativo

g3 <- lm(sr ~ pop15 + ddpi)

summary(g3)

Call:lm(formula = sr ~ pop15 + ddpi)

Residuals: Min 1Q Median 3Q Max -7.5831 -2.8632 0.0453 2.2273 10.4753

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 15.59958 2.33439 6.682 2.48e-08 ***pop15 -0.21638 0.06033 -3.586 0.000796 ***ddpi 0.44283 0.19240 2.302 0.025837 * ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Per la scelta del modello si può far ricorso all’AIC

AIC(g) AIC(g2) AIC(g3)[1] 282.1961 [1] 280.3414 [1] 281.8861

Il comando steo consente di effettuare un’analisi di regressione stepwise, basata sul criterio AIC.

Nell’argomento direction si può indicare se usare una procedura:

Backward Forward Oppure entrambe (both)

gstep <- step(g, direction="backward")

Start: AIC=138.3sr ~ pop15 + pop75 + dpi + ddpi

Df Sum of Sq RSS AIC- dpi 1 1.893 652.61 136.45<none> 650.71 138.30- pop75 1 35.236 685.95 138.94- ddpi 1 63.054 713.77 140.93- pop15 1 147.012 797.72 146.49

Step: AIC=136.45sr ~ pop15 + pop75 + ddpi

Df Sum of Sq RSS AIC<none> 652.61 136.45- pop75 1 47.946 700.55 137.99- ddpi 1 73.562 726.17 139.79- pop15 1 145.789 798.40 144.53

4. Modello di Regressione lineare con R parte - 2 Enrico Properzi -...

Documents

Transcript of 4. Modello di Regressione lineare con R parte - 2 Enrico Properzi -...