Tesina sul caso di studio della regressione in R

Il p-value viene utilizzato per una valutazione più veloce ed è valido per entrambi i modelli e

rispetto all’ipotesi nulla è una misura del grado di disaccordo e più è piccolo il p-value e maggiore

sarà l’evidenza contro l’ipotesi nulla.

è minore di α si può rifiutare l’ipotesi nulla e considerare le variabili

Mentre se il p-value

significative. 2

Dal programma R ci viene offerto tramite la grafica di output gli asterischi (*) che si trovano

accanto alla variabili stimate e più è grande il numero di asterischi e più sarà maggiore la

significatività delle stesse.

Commento sullo svolgimento del lavoro e Analisi conoscitiva dei dati

Dataset inviato dalla professoressa Scaccia riguardante sei diverse variabili osservate su 52

professori di un college americano.

Le variabili sono indicate di seguito:

# sx = sesso, codificato come 1 per femmina e 0 per maschio;

# rk = ruolo, codificato come 1 per ricercatore, 2 per professore associato, 3 per professore

ordinario;

# yr = numero di anni da cui si riveste il ruolo attuale;

# dg = titolo di studio, codificato come 1 per dottorato e 0 per master;

# yd = numero di anni dal conseguimento del titolo di studio;

# sl = salario annuo, in dollari.

Prima di incominciare a costruire il modello di regressione multipla, voliamo studiare la relazione

che legano la 5 variabili presenti nel dataset alla variabile salario.

Cominciamo con dei grafici a dispersione dei dati.

In R, innanzitutto ho caricato il dataset:

> library(foreign)

> dati = read.dta("salary.dta")

> dati

Poi ho reso disponibili le variabili:

> attach(dati) 3

> reg1=lm(sl ~ sx)

> summary(reg1)

> names(reg1)

> plot(sx,sl)

> confint(reg1) 4

> reg2=lm(sl ~ rk)

> summary(reg2)

> names(reg2)

> plot(rk,sl)

> confint(reg2) 5

> reg3=lm(sl ~ yr)

> summary(reg3)

> names(reg3)

> plot(yr,sl)

> abline(reg3$coe,col="red")

> confint(reg3) 6

> reg4=lm(sl ~ dg)

> summary(reg4)

> names(reg4)

> plot(dg,sl)

> abline(reg4$coe,col="red")

> confint(reg4) 7

> reg5=lm(sl ~ yd)

> summary(reg5)

> names(reg5)

> plot(yd,sl)

> abline(reg5$coe,col="red")

> confint(reg5) 8

Commento:

Per quanto riguarda le regressioni lineare semplice non hanno molto senso, mentre i grafici sono

molto significativi.

Nel primo grafico abbiamo due box plot che ci dicono la relazione tra il salario e il sesso per quanto

riguarda le donne c’è un dato anomalo in quanto non possiamo controllare se è un errore o se è la

verità, lo trattiamo come valore anomalo che è.

Non possiamo dedurre se ci sia una discriminazione salariale legata al sesso in quanto può

dipendere da una terza varabile che non è stata presa in considerazione.

Nel secondo grafico abbiamo tre box plot che ci dicono la relazione tra il salario e il ruolo e

possiamo dedurre che per quanto riguarda gli assistenti abbiamo la presenza di due dati anomali.

Nel terzo grafico abbiamo una regressione lineare che ci dice la relazione tra il salario e il numero

di anni da cui si riveste il ruolo attuale ed è significativa in quanto abbiamo una correlazione

positiva cioè maggiore di zero.

Nel quarto grafico abbiamo una regressione lineare che ci dice la relazione tra il salario e il titolo di

studio non è molto significativa in quanto abbiamo una correlazione negativa cioè minore di zero.

Ed è la relazione che meno ci interessa .

Ed infine abbiamo il quinto grafico con una regressione lineare che ci dice la relazione tra il salario

e il numero di anni dal conseguimento del titolo di studio significativa in quanto abbiamo una

correlazione positiva cioè maggiore di zero.

Inoltre si possono riassumere i dati attraverso degli indici descrittivi.

> summary(dati) 9

Osserviamo la distribuzione univariata.

> plot(density(sl))

Commento:

Apparte il grafico del salario gli altri non sono significativi. In quanto non si hanno una funzione di

densità per una variabile discreta (anche nel caso dicotomico) ma una funzione di probabilità.

Effettuiamo un test per vedere se tali correlazioni siano significativamente diverse da zero.

> cor.test(yr,sl)

In questo caso r= 0.700669 e il suo p-value=7.341e-09 ovvero non significativo per alfa= 0.05.

> cor.test(dg,sl)

In questo caso r= -0.06972576 e il suo p-value=0.06233 ovvero non significativo per alfa= 0.05. 10

> cor.test(yd,sl)

In questo caso r= 0.6748542 e il suo p-value=4.102e-08 ovvero non significativo per alfa= 0.05.

Si vuole studiare come varia il salario in funzione delle altre variabili.

In base ai dati precedenti , ho costruito un modello di regressione lineare multiplo con:

 Salario : variabile dipendente Y;

 Sesso : variabile indipendente X;

 Ruolo : variabile indipendente X;

 Numero di anni da cui si riveste il ruolo attuale: variabile indipendente X;

 Titolo di studio: variabile indipendente X;

 Numero di anni dal conseguimento del titolo di studio: variabile indipendente X.

Questo è il modello che voglio stimare:

= β β + Ruoloβ + N.a.d.c.s.r.i.l.r.a.β + N.d.a.d.c.d.t.d.s.β

Salario + Sesso + Titoli di studioβ .

0 1 2 3 4 5

Successivamente, si procede con la stima del modello:

> reg=lm(sl ~ factor(sx)+factor(rk)+yr+dg+yd)

> reg 11

Con la funzione summary otteniamo il risultato della regressione multipla:

> summary(reg)

Nell’output sopra , troviamo i valori relativi ai residui: valore minimo (Min) uguale a -4045.2,

valore massimo (Max) uguale a 9193.1, mediana (Median) uguale a -361.5 e 1° (1Q) uguale a

-1094.7 e 3° (Q3) uguale a 813.2 sono i quartili . dei parametri β relativi a

Abbiamo la parte relativa ai coefficienti, la prima colonna è quella

intercept (intercetta) ,sx,rk,yr ,dg e yd .

La colonna Estimate è quella relativa alle stime dei singoli parametri β, e come possiamo notare, il

valore dell’intercetta (β ) e dei coefficienti di regressione è la seguente:

0

β = 17134.66;

0

β = 1166.37;

1

β = 5292.36+11118.76;

2

β = 476.31;

3

β = -1388.61;

4

β = -124.57.

5

Pertanto il modello di regressione multipla è il seguente:

Y = 17134.66 + 1166.37x + 16411,12x + 476.31x -1388.61x -124.57x .

1 2 3 4 5

Da questo modello si può osservare che i regressori dg e yd hanno un valore negativo, quindi

significa che essi influenzano negativamente il salario, mentre gli altri regressori hanno effetti

positivi sul salario. β;

La colonna Std. Error è relativa alla stima delle singole varianze dei parametri la colonna dei t

value fornisce i valori dati dalla statistica test, che forniscono informazioni sulla significatività dei

singoli parametri e la colonna dei p-value, i cui valori sono importanti per capire la significatività

dei parametri.

Nel nostro caso abbiamo valori dei p-value bassi soprattutto per le variabili rk associate , rk Full, yr

con *** asterischi. Nel complesso sono tutte utili nello spiegare il comportamento della Y. 12

Poi troviamo il Residual standard error pari a 2398 on 45 degrees of freedom, che indica la stima

C’è il Multiple R-

della varianza dei residui (con 45 gradi di libertà dati da n-p-1=52-6-1=45).

squared ovvero R2, è un indicatore molto utile per capire quanta devianza nei dati osservati è stata

colta dal modello, nell’output è 0.855 , essendo un valore molto alto vuol dire che il grado di

accostamento del modello ai dati osservati è molto alto. Si può far riferimento anche all’R2

aggiustato ovvero Adjusted R-squared, si tratta di un R2 corretto in base ai gradi di libertà,

nell’output pari a 0.8357, sempre un valore alto che dà conferma di quanto detto per l’R2 , ma

leggermente più basso perché corretto in base ai gradi di libertà.

L’ultima riga del nostro output è relativa all’ANOVA ( analisi della varianza), che ci permette di

capire se le variabili X contribuiscono a spiegare la variabile Y. Quindi è utile per valutare il nostro

modello di regressione.

Troviamo la F-statistic, sarebbe la F calcolata nel seguente modo [(SSR/p)/(SSE/(n-p-1)],

trattandosi di un modello di regressione lineare multiplo, nell’output assume il valore di 44.24 su 6

e 45 gradi di libertà (p=6 e n-p-1=45).

L’ANOVA ci fornisce anche il p-value < 2.2e-16, si tratta di un p-value bassissimo.

d’ipotesi dell’ANOVA :

Considerando la batteria

H0 :modello di regressione con solo l’intercetta

H1:modello di regressione completo.

bassissimo, allora rifiuto H0 e accetto H1. Questo conferma l’importanza delle

Per un p-value

variabili esplicative nello spiegare la variabile dipendente.

Costruzione dell’intervallo stimato con livello di confidenza (1-α) pari a 0.95 con α=0.05:

> confint(reg) 13

Per valutare il rispetto delle ipotesi sottostanti al modello è possibile far riferimento ai 4 grafici dei

residui. Per la costruzione dei grafici su R si procede nel seguente modo:

> par(mfrow=c(2,2))

> plot(reg)

Oppure:

> layout(matrix(1:4,2,2))

> plot(reg)

Il primo grafico (Residuals vs Fitted) è il grafico dei residui rispetto ai valori previsti, dimostra una

situazione di variabilità costante, che ne conferma sia l’omoschedasticità e sia che la media degli

errori è nulla. E tre osservazioni (2,20,24) sono lontane dai valori stimati.

Nel Normal Q-Q viene dimostrata la normalità degli errori e si evidenzia una quasi perfetta

aderenza delle osservazioni alla linea della distribuzione normale, perché sono tutti concentrati sulla

linea Q-Q ad eccezione di 2 e 24 che sono fuori.

Lo scale-location ci permette di notare che la distribuzione degli errori è piuttosto concentrata sul

range dei valori previsti.

Il grafico residuals vs leverage mostra se ci sono dati anomali che riescono ad influenzare la stima

del modello. Sulla trama di punti vengono sovrapposte delle curve di livello per la distanza di Cook.

A piccole distanze significa che l’assenza di una determinata osservazione non influenza il modello,

invece grandi distanze indicano la presenza di dati anomali. 14

Adesso osserviamo il problema della elevata multicollinearità.

Siamo in presenza di elevata multicollinearità quando due o più regressori sono fortemente

correlati. In presenza di multicollinearità, alcune stime potrebbero presentare errori standard molto

grandi ed essere quindi molto imprecise.

Proviamo a vedere se tra i regressori