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 lime 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 gia' 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 e' un po' anomalo!

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

no.out=anscombe[-3,]

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

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

#ristimo quindi il modello con il database risotto, 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 piu' outlier possono inficiare la nostra stima e portarci

# a conclusioni fuorvianti. Conviene percio' capire quando un dato e' da

# escludere o meno.

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

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

# la retta interpola perfettamente tutti i dati.

# 14. Modello con variabile esplicativa dummy ---------------------------------

# Il regressore e' una variabile binaria.

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

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

#Remark: read.csv() e' la funzioni in R per importare automaticamente un dataset,

#di cui specifico:

#header = T. Cioe' la prima riga avra' i nomi delle variabili e non

#verra' 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 e' pari a 20 o e' 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 puo' essere vista come la pendenza della retta di regressione,

# dal momento che "dummy.str" non e' una variabile continua, bensi'

# una variabile binaria (puo' 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 puo' valere sia 1, sia 0, dove il pedice "i"

# sta ad indicare l'i-esima osservazione) e' 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 e' pari a 0, a prescindere da quanto

# vale D[i], cioe': E(e|D[i])=0

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

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

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

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

# Piuttosto che come coefficiente angolare, beta1 puo' essere visto, dunque,

# come la differenza tra le medie condizionate.

# Quindi, se beta1 e' significativamente diverso da zero, cio' 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

# e' maggiore o uguale a 20.

# beta1 e' significativamente diverso da 0 in quanto la statistica test

# pari a 4 e' maggiore di 1.96 (ossia il valore che mette davanti a se' solo

# lo 0.025 della massa di probabilita')

# 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", pero'

# questa volta regrediamo la variabile risposta (testscr) su

# piu' variabili esplicative:

# str (rapporto studenti/insegnanti)

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

# per facilita' 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 significativita' del 5%, tutti i coefficienti sono

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