Appunti statistica per l'economia e regressione su R

Statistica per l'Economia (Marilena Barbieri)

***********************************************************************

LAB 05 - Linear & multiple linear regression

28/11/2018

Nina Deliu

nina.deliu@uniroma1.it

***********************************************************************

Agenda

------------------------------------------------------------------

# ***************************************************************

• 1. R Studio interface (LAB01)
• 2. Some starting examples: data and operators (LAB01)
• 3. Data can have names (LAB01)
• 4. First data structure: vectors (LAB01)
• 5. Matrices (LAB02)
• 6. The most general data structure: List (LAB02)
• 7. Dataframes (LAB02)
• 8. Import and work with a Dataset (LAB02, LAB03)
• 9. Introduction to R graphics (LAB03)
• 10. The workspace (LAB01, LAB02)
• 11. Serie storiche univariate e multivariate (LAB03)
• 12. La regressione lineare semplice (LAB04)

## Last time we stopped here!

Today's topic

• 13. La regressione lineare semplice: gli OUTLIERS (LAB05)
• 14. La regressione lineare semplice: variabile esplicativa DUMMY (LAB05)
• 15. Regressione lineare multipla (LAB05)

# ***************************************************************

Ricorda!

# il cancelletto indica un commento

# per eseguire il codice andare sulla riga da eseguire e premere CTRL+Invio

13. La regressione lineare semplice: gli outliers

----------------

# A tal proposito, prendiamo un dataset già esistente in R:

?data

data(anscombe) # il comando data() carica nel nostro environment il dataset "anscombe" presente in R

str(anscombe)

1

# Modello di regressione lineare semplice:

# ricordiamo che la funzione R da usare per la regressione lineare è lm()

# lm() richiede la specificazione del modello (y ~ x)

mod1 <- lm(anscombe$y3 ~ anscombe$x3)

summary(mod1)

# vediamo ora la relazione tra la variabile dipendente y e quella indipendente x

# con un grafico a dispersione

plot(anscombe$x3, anscombe$y3) # con "pch=20" abbiamo dei pallini pieni

abline(coefficients(mod1),col=2)

# notiamo un valore particolarmente anomalo rispetto agli altri

# visualizziamo quindi i residui:

plot(residuals(mod1), xlab="Osservazione", ylab = "Residui")

# il residuo per la terza osservazione è un po' anomalo!

# Escludiamo quindi la terza osservazione, creando un database ridotto:

no.out=anscombe[-3,]

# Remark: il -3 prima della virgola dentro le parentesi quadre

# mi fa escludere la riga n.3, ossia la terza osservazione

# ristimo quindi il modello con il database ridotto, senza l'outlier

mod1.noout <- lm(no.out$y3 ~ no.out$x3)

summary(mod1.noout)

# verifico di nuovo l'andamento dei residui

plot(mod1.noout$residuals, xlab="Osservazione", ylab = "Residui")

# vediamo adesso in un unico grafico i miei dati e la retta di regressione stimata

# con il database completo (compreso outlier) e database ridotto (senza outlier)

plot(anscombe$x3, anscombe$y3, main = "Stima della retta di regressione",xlab = "X", ylab="Y")

abline(coefficients(mod1),col=2)

# Remark: abline() aggiunge una retta ad un grafico già esistente

# dandogli i due coefficienti, abline() capisce i parametri della retta

abline(coefficients(mod1.noout),col=3)

# legend ci inserisce una legenda nel grafico esistente

legend("topleft", c("Con outlier", "Senza outlier"),col = c("red", "green"), lwd = 1)

# Quindi, uno o più outlier possono inficiare la nostra stima e portarci

# a conclusioni fuorvianti. Conviene perciò capire quando un dato è da

# escludere o meno.

# Nel caso in questione, eliminare l'outlier vuol dire incrementare

# sensibilmente l'R², che addirittura assume valore 1, in quanto

2

# la retta interpola perfettamente tutti i dati.

14. Modello con variabile esplicativa dummy

---------------------------------

# Il regressore è una variabile binaria.

# Riprendiamo il dataset "cascschool" che abbiamo analizzato la volta precedente.

dati <- read.csv("caschool.csv", header=T, dec=".", sep=",")

# Remark: read.csv() è la funzione in R per importare automaticamente un dataset,

# di cui specifico:

# header = T. Cioè la prima riga avrà i nomi delle variabili e non

# verrà considerata come osservazione

# dec = "." e sep = ","

# Nota bene: quando si importa il db deve stare nella stessa cartella

# dello script R che abbiamo aperto

# Creiamo ora una Funzione indicatrice:

# ossia, una funzione che vale 1 se il rapporto studenti/insegnanti < 20

# e vale 0 altrove (quindi, quando quest'ultimo è pari a 20 o è maggiore)

dummy.str = I(dati$str<20)

mod2 <- lm(dati$testscr ~ dummy.str)

summary(mod2)

# Se plottiamo il diagramma a dispersione tra la variabile binaria e

# la variabile dipendente (punteggio medio dei test) e sovrapponiamo una retta,

# capiamo che, in tal caso, tale grafico non ha molto senso, in quanto

# beta1 non può essere vista come la pendenza della retta di regressione,

# dal momento che "dummy.str" non è una variabile continua, bensì

# una variabile binaria (può assumere solo due valori)

plot(dummy.str, dati$testscr, xlim=c(-1,2), ylim=c(600,720), pch=20, ylab="Punteggio test")

abline(coefficients(mod2), lwd=3, col=2)

# Quindi, quando D[i] (che può valere sia 1, sia 0, dove il pedice "i"

# sta ad indicare l'i-esima osservazione) è il regressore, abbiamo

# il seguente modello di regressione:

# y[i] = beta0 + beta1*D[i] + e[i] dove i=1,...,n

# Scindiamo in due casi diversi, quindi si ha che:

# y[i] = beta0 + e[i] se D[i]=0

# y[i] = beta0 + beta1 + e[i] se D[i]=1

# La media condizionata dell'errore è pari a 0, a prescindere da quanto

# vale D[i], cioè: E(e|D[i])=0

# Quindi, la media condizionata di Y[i] ove D[i]=1 è pari a:

# E(Y[i]|D[i]=1)=beta0+beta1

3

# Invece, la media condizionata di Y[i] ove D[i]=0 è pari a:

# E(Y[i]|D[i]=0)=beta0

# Piuttosto che come coefficiente angolare, beta1 può essere visto, dunque,

# come la differenza tra le medie condizionate.

# Quindi, se beta1 è significativamente diverso da zero, ciò vuol dire che

# avere un rapporto studenti/insegnanti minore di 20 comporta una media di

# risultati ai test diversa da quella che si ha quando lo stesso rapporto

# è maggiore o uguale a 20.

# beta1 è significativamente diverso da 0 in quanto la statistica test

# pari a 4 è maggiore di 1.96 (ossia il valore che mette davanti a sé solo

# lo 0.025 della massa di probabilità)

# Intervalli di confidenza dei parametri del modello al 95%

confint(mod2) # l'intervallo di confidenza esclude lo 0

15. Regressione lineare multipla

--------------------------------------------------------------

# Lavoriamo di nuovo sul dataset "caschool", però

# questa volta regrediamo la variabile risposta (testscr) su

# più variabili esplicative:

# str (rapporto studenti/insegnanti)

# el_pct (percentuale di studenti che studiano l'inglese)

# per facilità estraiamo queste variabili e mettiamole in un dataset ridotto:

dati.sub=dati[,c(11,14,16)]

# Facciamo una breve analisi esplorativa, usando il comando pairs() che mi

# crea dei grafici di dispersione per tutte le coppie di variabili presenti nel db

?pairs

pairs(dati.sub)

# vediamo inoltre la correlazione tra le variabili

cor(dati.sub)

# Passiamo quindi al modello di regressione lineare multipla;

# usiamo lo stesso comando della regressione lineare semplice,

# modificando un po' la formula, ossia, aggiungendo altre esplicative

# nella seguente struttura lm( y ~ x1 + x2 ), dove x1 e x2 sono le due covariate

mod3=lm(dati$testscr ~ dati$str + dati$el_pct)

summary(mod3)

# Commenti all'output di mod3:

# Vediamo innanzitutto che i coefficienti stimati ora sono 3:

# beta0 = intercetta; beta1 = coefficiente di x1; beta2 = coefficiente di x2

# con i rispettivi errori standard, t e p-values

4

# con un livello di significatività del 5%, tutti i coefficienti sono

# significativi (p-value < 5%) per t