Codice Vittadini
Alessandro Risaro
Fase di setting dell’ ambiente R
Comando per ripulire l’environment
rm(list=ls())Caricamento delle librerie utili alla nostra analisi
library(car)
library(describedata)
library(skedastic)
library(psych)## Warning: package 'psych' was built under R version 4.0.5library(klaR)## Warning: package 'MASS' was built under R version 4.0.5library(olsrr)
library(sandwich) #importata in companies e anche in countries
#library(het.test) #importata in hartangel
library(DataCombine)#importata in hartangel
library(systemfit) #importata in countries,cigarettes e anche in hsb
library(lmtest) #importata in countries,cigarettes e anche in hsbFunzione white test del prof (meno restrittiva di quella di R) White test che valuta l’eteroschedasticità basandosi sul modello di regressione in cui come variabile indipendente si ha il quadrato dei residui e come variabili indipendenti ci sono i fittedvalues e i fitted values al quadrato. Questo test è meno restrittivo delprecedente. Infatti, il white test precdente si basava su unaregressione dei residui al quadrato, su variabili indipendenti, i loro quadrati e le loro interazioni.
white.test<-function(lmod){
u2<-lmod$residuals^2
y<-lmod$fitted
R2u<-summary(lm(u2~y+I(y^2)))$r.squared
LM<-length(y)*R2u
p.val<-1-pchisq(LM,2)
data.frame("Test Statistic"=LM, "P"=p.val)
}Importare i dati (definire il separatore)
nazioni<-read.csv(file.choose(), sep=";")
companies<-read.csv(file.choose(), sep=';')
hartangel<-read.csv(file.choose(), sep=",")
countries<-read.csv(file.choose(),sep=';')
cigarettes<-read.csv(file.choose(),sep=' ')
hsb<-read.csv(file.choose(),sep=',')Preprocessing e pulizia dei dati
In hartangel notiamo che le variabili ftheft e mtheft hanno dei valori nulli e che i dati non sono ordinati nel tempo. Eliminiamo tutte le righe che hanno almeno una covariata missed e procediamo ad ordinarli temporalmente
hartangel<-na.omit(hartangel)
hartangel<-hartangel[order(hartangel$year),]Consegna: Regressioni lineari multiple di packpc sulle covariate cpi, pop, income e tax per gli stati Arkansas (state=”AR”) e California (state=”CA”)
Dobbiamo fare del subsetting in quanto ci vogliamo riferire solo ai dati relativi ad un determinato stato
cigarettesca<-cigarettes[cigarettes$state=='CA',] #tutte le osservazioni che hanno come stato la california
cigarettesar<-cigarettes[cigarettes$state=='AR',] #tutte le osservazioni che hanno come stato l'arizonaAndiamo ad aggiungere il prefisso CA ai nomi delle variabili che si trovano nel nuovo dataset cigarettesca e il prefisso AR alle variabili che si trovano nel nuovo dataset cigarettesar
names(cigarettesca)<-paste0(names(cigarettesca),'_CA')
names(cigarettesar)<-paste0(names(cigarettesar),'_AR')Aandiamo ad eliminare alcune colonne che contengono delle informazioni ridondanti
cigarettesca$state_CA<-NULL
cigarettesca$ID_CA<-NULL
cigarettesar$state_AR<-NULL
cigarettesar$ID_AR<-NULLAndiamo ad unire i due dataset
cigarettesarca<-cbind(cigarettesca,cigarettesar)Andiamo ad unificare la colonna year in una sola così da fare un po’ di pulizia del dataset E facciamo anche la stessa cosa con la variabile cpi (che è una costante)
cigarettesarca$year<-cigarettesca$year_CA
cigarettesarca$year_CA<-NULL
cigarettesarca$year_AR<-NULL
cigarettesarca$cpi<-cigarettesca$cpi_CA
cigarettesarca$cpi_CA<-NULL
cigarettesarca$cpi_AR<-NULLStatistiche descrittive
Guardare i dati
View(nazioni)
View(cigarettes)
View(hsb)Notiamo che relig viene tratta come una variabile numerica, mentre è una variabile categorica quindi procediamo a renderla categorica
nazioni$relig<-as.factor(nazioni$relig)Notiamo che in cigarettes la variabile cpi è una costante in ogni anno Notiamo che in hsb il dataset è composto da 5 variablili stringa e 5 variabili numeriche e 5 variabili stringa, più la variabile id (fattorizziamo le variabili gender e prog in quanto sono le uniche richieste nell’esercizio)
hsb$gender<-as.factor(hsb$gender)
hsb$prog<-as.factor(hsb$prog) #questo comando farà si che il nostro lm le riconoscerà come delle variabili dummyCreiamo un vettore contenente tutte le variabili numeriche
var_num_nazioni<-c("densita","urbana","vitamas","vitafem","alfabet","pil")
var_num_companies<-c('assets','mark_val','sales','profits','cash','employ')
var_num_countries<-c('Area','Irrigated','Population','Under.14','Life.expectancy','Literacy.Rate','Unemployment','ISPs.million','Tvs.person','Railways','Airports')
var_num_cigarettes<-c('cpi','pop','packpc','income','tax','avgprs','taxs')
var_num_hsb<-c('read','write','math','science','socst')Calcoliamo alcune statistiche descrittive (possiamo decidere se mettere solo var_num o anche relig)
summary(nazioni)## nazione densita urbana vitafem
## Length:66 Min. : 2.00 Min. : 5.00 Min. :43.00
## Class :character 1st Qu.: 19.75 1st Qu.:49.50 1st Qu.:70.00
## Mode :character Median : 61.00 Median :64.50 Median :76.00
## Mean :100.15 Mean :62.18 Mean :72.74
## 3rd Qu.:122.25 3rd Qu.:75.00 3rd Qu.:79.00
## Max. :605.00 Max. :96.00 Max. :82.00
## vitamas alfabet pil relig
## Min. :41.00 Min. : 27.00 Min. : 208 1:41
## 1st Qu.:64.00 1st Qu.: 83.50 1st Qu.: 1412 2: 8
## Median :69.00 Median : 95.50 Median : 4464 3:17
## Mean :66.58 Mean : 87.58 Mean : 7303
## 3rd Qu.:73.00 3rd Qu.: 99.00 3rd Qu.:14048
## Max. :76.00 Max. :100.00 Max. :23474summary(nazioni[,var_num_nazioni])## densita urbana vitamas vitafem
## Min. : 2.00 Min. : 5.00 Min. :41.00 Min. :43.00
## 1st Qu.: 19.75 1st Qu.:49.50 1st Qu.:64.00 1st Qu.:70.00
## Median : 61.00 Median :64.50 Median :69.00 Median :76.00
## Mean :100.15 Mean :62.18 Mean :66.58 Mean :72.74
## 3rd Qu.:122.25 3rd Qu.:75.00 3rd Qu.:73.00 3rd Qu.:79.00
## Max. :605.00 Max. :96.00 Max. :76.00 Max. :82.00
## alfabet pil
## Min. : 27.00 Min. : 208
## 1st Qu.: 83.50 1st Qu.: 1412
## Median : 95.50 Median : 4464
## Mean : 87.58 Mean : 7303
## 3rd Qu.: 99.00 3rd Qu.:14048
## Max. :100.00 Max. :23474Vediamo la distribuzione della variabile categorica in un altro modo e anche per hsb
table(nazioni$relig)/nrow(nazioni)##
## 1 2 3
## 0.6212121 0.1212121 0.2575758table(hsb$gender)##
## female male
## 109 91table(hsb$gender)/nrow(data)## numeric(0)table(hsb$prog)##
## academic general vocational
## 105 45 50table(hsb$prog)/nrow(data)## numeric(0)Descriviamo le correlazioni tra le variabili numeriche
pairs.panels(nazioni[,var_num_nazioni])
Vediamo una forte correlazione tra le variabili vitamas e vitafem e anche tre vitafem e alfabet e tra vitamas alfabet vi sono correlazioni meno forti, ma evidenti Le variabili non sembrano distribuirsi normalmente, in particolare, la variabile densità è caratterizzata da una distribuzione assimmetrica positiva, vitamas, vitafem e alfabet sono caratterizzate da una distribuzione assimmetrica negativa, mentre pil sembra presentare due mode (quindi c’è un grande gruppo di paesi poveri e un grande numero di paesi più ricchi)
Andiamo a vedere come fare sia i summary che i pairs panels condizionati
summary(hsb[hsb$gender=='female',var_num_hsb])## read write math science socst
## Min. :28.00 Min. :35.00 Min. :33.00 Min. :29.0 Min. :26.00
## 1st Qu.:44.00 1st Qu.:50.00 1st Qu.:45.00 1st Qu.:44.0 1st Qu.:46.00
## Median :50.00 Median :57.00 Median :53.00 Median :50.0 Median :56.00
## Mean :51.73 Mean :54.99 Mean :52.39 Mean :50.7 Mean :52.92
## 3rd Qu.:57.00 3rd Qu.:62.00 3rd Qu.:58.00 3rd Qu.:58.0 3rd Qu.:61.00
## Max. :76.00 Max. :67.00 Max. :72.00 Max. :69.0 Max. :71.00summary(hsb[hsb$gender=='male',var_num_hsb])## read write math science
## Min. :31.00 Min. :31.00 Min. :35.00 Min. :26.00
## 1st Qu.:44.50 1st Qu.:41.00 1st Qu.:45.00 1st Qu.:46.00
## Median :52.00 Median :52.00 Median :52.00 Median :55.00
## Mean :52.82 Mean :50.12 Mean :52.95 Mean :53.23
## 3rd Qu.:63.00 3rd Qu.:59.00 3rd Qu.:59.50 3rd Qu.:61.00
## Max. :76.00 Max. :67.00 Max. :75.00 Max. :74.00
## socst
## Min. :26.00
## 1st Qu.:46.00
## Median :51.00
## Mean :51.79
## 3rd Qu.:61.00
## Max. :71.00Vediamo che per esempio nella lettura i maschi ottengono un punteggio medio più alto
par(mfrow=c(1,2))
pairs.panels(hsb[hsb$gender=='male',var_num_hsb])
pairs.panels(hsb[hsb$gender=='female',var_num_hsb])
## Statistiche descrittive per serie temporali In hartangel che era una serie temporale quando indaghiamo le correlazioni dobbiamo togliere la variabile years
pairs.panels(hartangel[,-1])
Andiamo a vedere l’andamento della nostra variabile dipendente nel tempo
plot(hartangel$year,hartangel$ftheft,ylab="ftheft",xlab='year',type='b')
abline(h=mean(hartangel$ftheft))
Vediamo che dal 1955 in poi inizia un trend crescente, mentre prima del ’55 si vedono fluttuazioni intorno al valore 20, inoltre in questo range di fluttuazioni si vedono valori che tendono ad essere uguali tra di loro che indicano la possibile presenza di autocorrelazione positiva.
Andiamo a verificare la presenza di una correlazione positiva attraverso un correlogramma
par(mfrow=c(2,1))
acf(hartangel$ftheft,main='autocorr. variabile dipendente')
pacf(hartangel$ftheft,main='autocorr. parziale variabile dipendente')
Il primo grafico indica l autocorrelazione della variabile dipendente a diversi lag temporali, vediamo che vi sono autocorrelazioni significative fino al quarto lag (ricordiamo l’ACF al lag 0 non va mai letta) Il secondo grafico ci segnala la presenza di autocorrelazione significativa di ordine 1 Quindi possiamo concludere che vi sia almeno un autocorrelazione positiva di primo ordine.
Modelli lineari classici multipli: diagnostica ipotesi e correlazione tra i residui
Soluzione eteroschedasticità apportando trasformazioni alle variabili dipendenti ed indipendenti
Iniziamo svolgendo una regressione lineare tra le variabili dipendenti e indipendenti richieste (alt+126= ~)
lin_nazioni<-lm(pil~alfabet, nazioni)
summary(lin_nazioni) ## vediamo che il model##
## Call:
## lm(formula = pil ~ alfabet, data = nazioni)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7461 -4800 -2668 5282 13912
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13689.5 3919.6 -3.493 0.000873 ***
## alfabet 239.7 44.0 5.448 8.69e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5826 on 64 degrees of freedom
## Multiple R-squared: 0.3168, Adjusted R-squared: 0.3061
## F-statistic: 29.68 on 1 and 64 DF, p-value: 8.686e-07Vediamo che la variabile indipendente alfabet risulta essere significativa, la statistica F (in questo caso uguale alla statistica t) indica che il modello nel suo complesso è significativo ed abbiamo un R^2 pari a 0.3 il quale indica uno scarso adattamento ai dati
Andiamo a verificare la presenza di eteroschedasticità
white.test(lin_nazioni)## Test.Statistic P
## 1 13.82682 0.0009943592white_lm(lin_nazioni)## # A tibble: 1 x 5
## statistic p.value parameter method alternative
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 13.8 0.000994 2 White's Test greaterVediamo che il test di white segnala la presenza di eteroschedasticità dei residui Verifichiamo la presenza di eteroschedasticità anche in maniera grafica
plot(lin_nazioni$fitted.values, rstandard(lin_nazioni), xlab='fitted',ylab='std.residuals') 
Notiamo che i residui hanno un andamento a ventaglio, quindi si può desumere anche graficamente di essere in presenza di eteroschedasticità.
Vediamo ora alcune trasformazioni applicate alle variabili dipendenti e indipendenti e come esse agiscono sulla bontà di adattamento del modello, sulla significatività delle variabili esplicative e sull’eteroschedasticità dei residui:
- trasformazione lin-log
linlog_nazioni<-lm(pil~I(log(alfabet)), nazioni)
summary(linlog_nazioni)##
## Call:
## lm(formula = pil ~ I(log(alfabet)), data = nazioni)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7303 -4929 -2475 5953 14331
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -57490 13869 -4.145 0.000102 ***
## I(log(alfabet)) 14566 3113 4.679 1.54e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6084 on 64 degrees of freedom
## Multiple R-squared: 0.2549, Adjusted R-squared: 0.2432
## F-statistic: 21.89 on 1 and 64 DF, p-value: 1.543e-05Il coefficente 14566 va interpretato come un incremento del 1% del tasso di alfabetizzazione porta ad un incremento di 145 dollari del PIL (quandoo abbiamo solo il log nella var ind dividiamo per 100) vediamo che il modello è ancora significativo, ma abbiamo perso di R^2 quindi possiamo dire che questo modello interpreta peggio i dati rispetto a quello lineare
Andiamo di nuvo a controllare la presenza di eteroschedasticità dei residui:
plot(linlog_nazioni$fitted.values, rstandard(linlog_nazioni), xlab='fitted',ylab='std.residuals', main="linlog")
white_lm(linlog_nazioni) ## # A tibble: 1 x 5
## statistic p.value parameter method alternative
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 14.6 0.000670 2 White's Test greaterVediamo che i residui si dispongono a ventaglio, quindi graficamente notiamo la presenza di eteroschedasticità. Che viene anche confermata dal White Test
- Trasformazioni loglog
loglog_nazioni<-lm((log(pil))~I(log(alfabet)), nazioni)
summary(loglog_nazioni) ##
## Call:
## lm(formula = (log(pil)) ~ I(log(alfabet)), data = nazioni)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6880 -0.6551 -0.1290 0.8736 2.3993
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.2166 2.0702 -4.452 3.48e-05 ***
## I(log(alfabet)) 3.9268 0.4647 8.450 5.20e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9081 on 64 degrees of freedom
## Multiple R-squared: 0.5273, Adjusted R-squared: 0.52
## F-statistic: 71.41 on 1 and 64 DF, p-value: 5.196e-12Vediamo come la variabile alfabet rimane significativa e la capacità del modello di adattarsi ai dati migliora di molto R^2=0.5 Il coefficente adesso si interpreta come elasticità, ovvero un rapporto tra variazioni percentuali (ovvero l’aumento del 1% del tasso di alfabetizzazione è associato ad un aumento del 3.9% del pilprocapite)
Andiamo a valutare la presenza di eteroschedasticità:
plot(loglog_nazioni$fitted.values, rstandard(loglog_nazioni), xlab='fitted',ylab='std.residuals', main="loglog")
white_lm(loglog_nazioni)## # A tibble: 1 x 5
## statistic p.value parameter method alternative
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 22.5 0.0000128 2 White's Test greaterDa questo grafico si vede un grande miglioramento dei residui, ma non posso dire nulla sull’eteroschedasticità, mentre il test di white è ancora significativo e quindi gli errori secondo il test sono ancora eteroschedastici.
- trasformazioni loglin
loglin_nazioni<-lm((log(pil))~alfabet, nazioni)
summary(loglin_nazioni)##
## Call:
## lm(formula = (log(pil)) ~ alfabet, data = nazioni)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6371 -0.5733 -0.1513 0.7705 1.7647
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.802183 0.556619 5.034 4.16e-06 ***
## alfabet 0.062221 0.006249 9.958 1.25e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8273 on 64 degrees of freedom
## Multiple R-squared: 0.6077, Adjusted R-squared: 0.6016
## F-statistic: 99.16 on 1 and 64 DF, p-value: 1.246e-14Vediamo come il modello rimane significativo e l’ R^2 migliora ancora arrivando a 0.6 Il coefficente di alfabet cambia interpretazione ovvero aumento percentuale del pil a fronte di un aumento unitario del tasso di alfabetizzazione (un aumento del 0.06 di alfabet porta ad un aumento del 1% del pil)
Andiamo a verificare se vi è ancora o meno la presenza di eteroschedasticità
plot(loglin_nazioni$fitted.values, rstandard(loglin_nazioni), xlab='fitted',ylab='std.residuals', main="loglin")
white_lm(loglin_nazioni)## # A tibble: 1 x 5
## statistic p.value parameter method alternative
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 3.62 0.163 2 White's Test greaterGraficamente vediamo che l’andamento dei residui sembra migliorato, in quanto sembrano essere più omoschedastici però non si capisce molto bene, notiamo che il p-value è alto e ciò porta all’accettazione dell’ipotesi nulla di omoschedasticità. Quindi per bontà di adattamento e per avere riportato gli errori ad essere omoschedastici questo modello sembra essere il migliore.
- trasformazioni polinomiali (in questo caso di secondo grado)
quad1_nazioni<-lm(pil~alfabet+I(alfabet^2), nazioni)
summary(quad1_nazioni)##
## Call:
## lm(formula = pil ~ alfabet + I(alfabet^2), data = nazioni)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7875 -4754 197 4711 13259
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24473.855 11053.547 2.214 0.030447 *
## alfabet -891.719 312.555 -2.853 0.005852 **
## I(alfabet^2) 7.678 2.103 3.650 0.000534 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5335 on 63 degrees of freedom
## Multiple R-squared: 0.4361, Adjusted R-squared: 0.4182
## F-statistic: 24.36 on 2 and 63 DF, p-value: 1.456e-08Vediamo che l’andamento quadratico di secondo grado è significativo in quanto sono significativi sia i due coefficenti, sia il test F Gli standard error così elevati dicono che potrebbe essere stata introdotta multicollinearità tra le variabili indipendenti In questo caso però potremmo avere introdotto multicollinearità tra le variabili esplicative, andiamo quindi a controllare
vif(quad1_nazioni)## alfabet I(alfabet^2)
## 60.17541 60.17541Infatti il Vif di entrambe le variabili esplicative risulta essere molto maggiore di 10, questo stabilisce la presenza di multicollinearità tra le due variabili
Per eliminare la multicollinearità possiamo centrare le variabili indipendenti: - procediamo in un primo momento a creare la variabile centrata, che non è altro che la variabile meno la sua media
nazioni$alfabet_centrata<-nazioni$alfabet-mean(nazioni$alfabet)- Andiamo a rifare la regressioni usando la variabile centrata ed il suo quadrato
quad2_nazioni<-lm(pil~alfabet_centrata+I(alfabet_centrata^2), nazioni)
summary(quad2_nazioni)##
## Call:
## lm(formula = pil ~ alfabet_centrata + I(alfabet_centrata^2),
## data = nazioni)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7875 -4754 197 4711 13259
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5263.832 862.127 6.106 7.01e-08 ***
## alfabet_centrata 453.013 70.978 6.382 2.35e-08 ***
## I(alfabet_centrata^2) 7.678 2.103 3.650 0.000534 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5335 on 63 degrees of freedom
## Multiple R-squared: 0.4361, Adjusted R-squared: 0.4182
## F-statistic: 24.36 on 2 and 63 DF, p-value: 1.456e-08Notiamo che abbiamo dei t-value molto più ampi, in quanto nel caso di multicollinearità gli standard error erano molto alti rispetto al loro valore reale
plot(quad1_nazioni$fitted.values, quad2_nazioni$fitted.values, xlab='fitted quad1 (non cent)', ylab='fitted quad 2 (centr)')
Vediamo che i valori fittati sono pressochè identici
Andiamo a testare l’eteroschedasticità, in questo caso dato che abbiamo più di una variabile dipendente dobbiamo aggiungere le interazuibu all’interno del white test
white_lm(quad2_nazioni, interaction=T)## # A tibble: 1 x 5
## statistic p.value parameter method alternative
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 16.0 0.00697 5 White's Test greaterVediamo che anche in questo caso gli errori sono eteroschedastici.
Verifica della presenza di multicollinearità
Verificare la presenza di multicollinearità tra le variabili esplicative è una delle prime cose che deve essere fatta quando vogliamo fare una regressione che abbia più di una variabile esplicativa
mod1_nazioni<-lm(pil~densita+urbana+vitafem+vitamas+alfabet+relig, nazioni)
summary(mod1_nazioni)##
## Call:
## lm(formula = pil ~ densita + urbana + vitafem + vitamas + alfabet +
## relig, data = nazioni)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8367.9 -4094.3 -328.8 3226.2 12376.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22033.709 5251.591 -4.196 9.45e-05 ***
## densita 5.148 5.838 0.882 0.38153
## urbana 32.363 50.847 0.636 0.52697
## vitafem 265.059 411.114 0.645 0.52164
## vitamas 70.108 408.057 0.172 0.86419
## alfabet 21.886 91.730 0.239 0.81226
## relig2 -2461.099 2170.794 -1.134 0.26157
## relig3 4821.352 1609.681 2.995 0.00403 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5168 on 58 degrees of freedom
## Multiple R-squared: 0.5127, Adjusted R-squared: 0.4539
## F-statistic: 8.718 on 7 and 58 DF, p-value: 2.839e-07Vediamo che il modello ha un test F significativo, ovvero esiste almeno un coefficiente che è diverso da 0,in questo caso l’unica variabile significativa è la relig3 ovvero la presenza della religione protestante
Testiamo la presenza di multicollinearità
-Matrice di correlazione anche tramite il comando pairs.panels come visto nella statistica descrittiva
- VIF(Varianza multifattoriale)
ols_vif_tol(mod1_nazioni)## Variables Tolerance VIF
## 1 densita 0.90762091 1.101782
## 2 urbana 0.38357395 2.607059
## 3 vitafem 0.02355899 42.446641
## 4 vitamas 0.02904083 34.434274
## 5 alfabet 0.18109586 5.521937
## 6 relig2 0.80631027 1.240217
## 7 relig3 0.81683246 1.224241vif(mod1_nazioni)## GVIF Df GVIF^(1/(2*Df))
## densita 1.101782 1 1.049658
## urbana 2.607059 1 1.614639
## vitafem 42.446641 1 6.515109
## vitamas 34.434274 1 5.868072
## alfabet 5.521937 1 2.349880
## relig 1.364066 2 1.080709Considerando la soglia 10 vediamo un problema di multicollinearità per la variabili vitafem e vitamas, in quanto molto probabilmente l’aspettativa di vita delle femmine e dei maschi dipende dagli stessi fattori, quindi se ne deve inserire solo una delle due variabili esplicative
- Condition index
ols_eigen_cindex(mod1_nazioni)## Eigenvalue Condition Index intercept densita urbana
## 1 5.8435403371 1.000000 4.146107e-04 0.007592516 1.010306e-03
## 2 1.0089347275 2.406614 3.158112e-06 0.007629420 2.211371e-06
## 3 0.5998283855 3.121221 8.239702e-05 0.353333515 3.465049e-04
## 4 0.4865830514 3.465449 8.788833e-04 0.504208618 3.569386e-03
## 5 0.0459604379 11.275763 9.605069e-02 0.090457517 5.179907e-01
## 6 0.0113479022 22.692390 5.388246e-01 0.019160893 4.330651e-01
## 7 0.0035538597 40.549722 3.487615e-01 0.001104670 1.145356e-02
## 8 0.0002512986 152.490574 1.498409e-02 0.016512851 3.256222e-02
## vitafem vitamas alfabet relig2 relig3
## 1 1.274152e-05 1.543974e-05 1.741487e-04 0.0031184172 0.006063044
## 2 2.262155e-07 1.443358e-07 6.405641e-06 0.4874862559 0.173563754
## 3 9.647968e-07 1.573052e-06 1.259138e-06 0.1407552812 0.461093602
## 4 3.156164e-05 3.762416e-05 3.111197e-04 0.2331720590 0.239504437
## 5 1.530420e-04 2.453426e-04 6.369822e-09 0.0003791654 0.055211194
## 6 1.642043e-03 2.080159e-03 1.935465e-01 0.0287327096 0.006094018
## 7 2.328528e-02 4.330982e-02 7.529275e-01 0.0770575063 0.046892873
## 8 9.748741e-01 9.543099e-01 5.303314e-02 0.0292986054 0.011577078Abbiamo ben due autovalori associati ad un condition index con autovalori maggiori di 30, uno dei due spiega pIù del 95% della vairbilità di vitamas e vitafem, quindi alche il condition index indica la presenza di multicollinearità tra queste due variabili e si procede ad omettere una fra vitamas e vitafem, oppure si potrebbe calcolarne la media.
Decidiamo di togliere dal modello vitamas e andiamo di nuovo a testare la presenza di multicollinearità
mod2_nazioni<-lm(pil~densita+urbana+vitafem+alfabet+relig, nazioni)
summary(mod2_nazioni)##
## Call:
## lm(formula = pil ~ densita + urbana + vitafem + alfabet + relig,
## data = nazioni)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8357.1 -4109.6 -374.6 3196.2 12354.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22005.133 5205.608 -4.227 8.34e-05 ***
## densita 5.282 5.739 0.920 0.36115
## urbana 31.374 50.103 0.626 0.53360
## vitafem 330.796 149.148 2.218 0.03042 *
## alfabet 20.868 90.782 0.230 0.81899
## relig2 -2541.140 2102.705 -1.209 0.23167
## relig3 4836.980 1593.836 3.035 0.00358 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5126 on 59 degrees of freedom
## Multiple R-squared: 0.5125, Adjusted R-squared: 0.4629
## F-statistic: 10.34 on 6 and 59 DF, p-value: 8.403e-08ols_vif_tol(mod2_nazioni)## Variables Tolerance VIF
## 1 densita 0.9239501 1.082310
## 2 urbana 0.3885518 2.573659
## 3 vitafem 0.1760520 5.680141
## 4 alfabet 0.1818551 5.498882
## 5 relig2 0.8452392 1.183097
## 6 relig3 0.8194492 1.220332ols_eigen_cindex(mod2_nazioni)## Eigenvalue Condition Index intercept densita urbana
## 1 4.868516568 1.000000 5.922456e-04 0.011466635 1.464889e-03
## 2 1.008607534 2.197036 4.728968e-06 0.007217311 7.144921e-07
## 3 0.598738281 2.851542 8.847540e-05 0.394498541 3.437302e-04
## 4 0.466282786 3.231273 1.490959e-03 0.470182240 5.523430e-03
## 5 0.044805166 10.423997 1.186160e-01 0.093336018 4.970642e-01
## 6 0.010806096 21.225792 4.398915e-01 0.022271511 4.437282e-01
## 7 0.002243569 46.583120 4.393161e-01 0.001027744 5.187486e-02
## vitafem alfabet relig2 relig3
## 1 1.357979e-04 2.499443e-04 4.864762e-03 0.009207055
## 2 2.216398e-06 7.865957e-06 5.154375e-01 0.171221687
## 3 7.852879e-06 1.944954e-06 1.324095e-01 0.437789146
## 4 3.844941e-04 5.427798e-04 2.559910e-01 0.271025070
## 5 1.860275e-03 1.533274e-04 8.343093e-05 0.049003446
## 6 1.742514e-02 2.755957e-01 4.395841e-02 0.003406401
## 7 9.801842e-01 7.234484e-01 4.725535e-02 0.058347194Vediamo come anche la variabile vitafem diventa significativa e il nostro R^2 è aumentato Per i vif vediamo come vitafem e alfabet sono quelle più correlate, ma è minore di 10, quindi va bene così ols_eigen_cindex(mod2) notiamo come ci siano dei valori alti (>30 dell’ultimo autovalore), decidiamo di proseguire così anche se sarebbe ragionevole togliere una delle variabili vitafem, alfabet Il fatto di togliere o non togliere una variabile dipende sempre dal ragionamento teorico e esplicativo del nostro modello
Generalmente dopo aver risolto il problema di multicollinearità nel modello si va ad indagare se vi è eteroschedasticità o meno.
Studio degli outliers e dei valori influenti
Costruisco gli scalari k e n, pari rispettivamente al numero di variabili + intercetta (utilizzate nel modello su cui voglio indagare gli outliers) e al numero delle osservazioni. Questi saranno utili per andare a costruire le soglie.
k=length(coef(mod2_nazioni))
k## [1] 7n=nrow(nazioni)
n## [1] 66Vado ora ad indagare graficamente la presenza di outlier a livello grafico:
par(mfrow=c(2,2))
# primo grafico: distanza di cook vs indice delle osservazioni
plot(mod2_nazioni,which=4)
abline(h=4/n)
#aggiungiamo ora il secondo plot che è covratio
plot(covratio(mod2_nazioni),ylab='covratio')
abline(h=1+3*k/n)
abline(h=1-3*k/n)
text(covratio(mod2_nazioni))
#aggiungiamo ora il terzo plot studentized resisiduals vs leverage (anche detti hat values)
#aggiungendo anche le soglie che per gli stud residuals sono +2 e -2
#aggiungiamo anche le soglie di leverage che sono 2k/n
plot(hatvalues(mod2_nazioni),rstudent(mod2_nazioni),xlab='leverage',ylab='studentized res')
abline(h=2)
abline(h=-2)
abline(v=2*k/n)
text(hatvalues(mod2_nazioni),rstudent(mod2_nazioni))
#l'ultimo grafico da aggiungere è quello dei dffits che ha soglie 2*sqrt(k/n)
plot(dffits(mod2_nazioni), ylab='DFFITS')
text(dffits(mod2_nazioni))
abline(h=2*sqrt(k/n))
abline(h=-2*sqrt(k/n))
Le soglie vengono usate come linee guida, ma la decisione di quali osservazioni effettivamente eliminare viene presa a braccio cercando di sintetizzare le informazioni date dai vari grafici, eliminerò quindi osservazioni che sono troppo ‘fuori’ in un grafico oppure osservazioni che vengono ritenute outliers su più grafici
La cook’s distance è una misura di influenza, ovvero di quanto un’osservazione influisca sulle stime dei parametri, vediamo che l’osservazione 5 ha un valore molto influente, notiamo anche che le osservazioni6 e 15 superano le soglie e bisogna prestare attenzione (anche se non superano le soglie) a le osservazioni 60,62 e anche la 57
Il covratio è un’altra misura di influenza, notiamo osservazioni fuori dalle soglie come la 50, la 66 e la 60 (che era già stata notata precedentemente)
I residui studentizzati ennesima misura di influenza, indica che la 66, la 50 (entrambe già evidenziate dal covratio) e la 5 (già evidenziata dalla distanza di coook sono osservazioni influenti)
I levarage indicano la presenza di outlier nelle osservazoni 60,12 e 67
Il DFFITS ci segnala soprattuto l’osservazione 5 che è molto fuori
La decisione di escludere valori viene fatta escludendo inanzitutto gli outlier e poi fare a braccio e trovare quali osservazioni influenti escludere.
Ora uso dbetas, la quale è una misura di influenza della singola osservazione sul singolo coefficente, ovvero escludo l’intercetta (di cui non si può fare a meno nel modello, o se si decide di escluderla bisogna comunque essere cauti) per andare a vedere quali outliers vi sono per ogni variabili indipendente
dfbetasPlots(mod2_nazioni) 
dfbetasPlots(mod2) #commenti al minuto 56 Vediamo come l’osservazione 5 e 60 compaiono spesso,quindi anche i dfbetas ci danno alcune informaizioni che però avevamo già ottenuto dai grafici precedenti
Prima di eliminare le osservazioni che presentano valori anomali voglio testare la normalità dei residui
shapiro.test(mod2_nazioni$residuals)##
## Shapiro-Wilk normality test
##
## data: mod2_nazioni$residuals
## W = 0.96234, p-value = 0.04284par(mfrow=c(1,1))
plot(mod2_nazioni,which=2)
Per il test di shapiro i residui risultano essere non normali (p-value<0.05) Il QQ plot ci segnala la presenza di valori estrremi come la 5, la 12 e la 60 (e soprattutto nella parte inziale ci sono degli scostamenti rispetto alla normale)
Procedo ad eliminare gli outliers e le osservazioni influenti, come criteri usiamo un laverage > 2k/n e per gli studentized residuals usiamo valore assolute studentized residuals >2 (manteniamo quelli che non oltrepassano quindi queste soglie)
nazioni2<-nazioni[hatvalues(mod2_nazioni)<=2*k/n & abs(rstudent(mod2_nazioni))<2,]
mod3_nazioni<-lm(pil~densita+urbana+vitafem+alfabet+relig, nazioni2)confrontiamo le due regressioni con gli outliers e senza outliers
summary(mod2_nazioni)##
## Call:
## lm(formula = pil ~ densita + urbana + vitafem + alfabet + relig,
## data = nazioni)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8357.1 -4109.6 -374.6 3196.2 12354.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22005.133 5205.608 -4.227 8.34e-05 ***
## densita 5.282 5.739 0.920 0.36115
## urbana 31.374 50.103 0.626 0.53360
## vitafem 330.796 149.148 2.218 0.03042 *
## alfabet 20.868 90.782 0.230 0.81899
## relig2 -2541.140 2102.705 -1.209 0.23167
## relig3 4836.980 1593.836 3.035 0.00358 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5126 on 59 degrees of freedom
## Multiple R-squared: 0.5125, Adjusted R-squared: 0.4629
## F-statistic: 10.34 on 6 and 59 DF, p-value: 8.403e-08summary(mod3_nazioni)##
## Call:
## lm(formula = pil ~ densita + urbana + vitafem + alfabet + relig,
## data = nazioni2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10367.0 -2959.9 -491.3 2948.4 10408.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -19632.552 5572.814 -3.523 0.000899 ***
## densita 14.073 5.682 2.477 0.016547 *
## urbana 37.717 45.273 0.833 0.408600
## vitafem 287.755 146.204 1.968 0.054390 .
## alfabet 3.442 91.571 0.038 0.970158
## relig2 -1043.898 1804.111 -0.579 0.565342
## relig3 5995.940 1502.033 3.992 0.000207 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4212 on 52 degrees of freedom
## Multiple R-squared: 0.5968, Adjusted R-squared: 0.5503
## F-statistic: 12.83 on 6 and 52 DF, p-value: 7.82e-09il modello senza outliers spiega meglio il pil rispetto a quello con gli outliers
Studio dell’eteroschedasticità e risoluzione tramite il metodo FGLS
Settiamo la regressione sulla quale verrà poi studiata l’eteroschedasticità
mod1_companies<-lm(mark_val~assets+sales+profits+cash+employ,companies)
summary(mod1_companies)##
## Call:
## lm(formula = mark_val ~ assets + sales + profits + cash + employ,
## data = companies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1228.23 -222.80 -50.87 251.70 1180.29
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.438e+02 1.030e+02 2.367 0.0213 *
## assets 2.066e-04 2.127e-02 0.010 0.9923
## sales -1.802e-02 7.445e-02 -0.242 0.8096
## profits 8.828e-02 1.176e+00 0.075 0.9404
## cash 3.787e+00 8.620e-01 4.393 4.82e-05 ***
## employ 9.399e+00 4.875e+00 1.928 0.0587 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 444.8 on 58 degrees of freedom
## Multiple R-squared: 0.7394, Adjusted R-squared: 0.7169
## F-statistic: 32.91 on 5 and 58 DF, p-value: 9.393e-16Andiamo ad indagare la presenza di eteroschedasticità inizialmente a livello grafico
par(mfrow=c(2,2)) # (2,2) per comodità nella rappresentazione
plot(mod1_companies$fitted,mod1_companies$residuals,xlab='fitted',ylab='residuals')
abline(h=0)
plot(mod1_companies$fitted,(mod1_companies$residuals)^2,xlab='fitted',ylab='residuals square')
La forma a ventaglio con la quale si distribuiscono i residui denota la presenza di eteroschedasticità Andiamo a veder il grafico dei residui rispetto ad ognuna delle variabili indipendenti
par(mfrow=c(3,2))
var_indipendenti<-c('assets','sales','profits','cash','employ')
for (i in var_indipendenti){
plot(companies[,i],mod1_companies$residuals,xlab=i,ylab='residuals')
}
Anche questo grafico segnala la presenza di eteroschedasticità
Andiamo ora ad indagare la presenza di eteroschedasticità attraverso alcuni appositi test
white.test(mod1_companies)## Test.Statistic P
## 1 9.406023 0.009067927white_lm(mod1_companies,interactions = T)## # A tibble: 1 x 5
## statistic p.value parameter method alternative
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 36.5 0.0135 20 White's Test greaterEntrambi segnalano la presenza di eteroschedasticità
Vediamo ora la risoluzione del problema di eteroschedasticità attraverso lo stimatore FGLS, la quale si divide in step: - Costruire un modello che ha i residui al quadrato come variabile dipendendente e come x le variabili indipendenti del modello di regressione eteroschedastico
aux_companies<-lm(I(mod1_companies$residuals^2)~assets+sales+profits+cash+employ,companies)
summary(aux_companies)##
## Call:
## lm(formula = I(mod1_companies$residuals^2) ~ assets + sales +
## profits + cash + employ, data = companies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -358259 -146663 -33223 46966 1198288
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 103852.98 64252.30 1.616 0.1114
## assets -21.87 13.26 -1.649 0.1046
## sales -10.26 46.44 -0.221 0.8259
## profits -571.74 733.26 -0.780 0.4387
## cash 1267.04 537.70 2.356 0.0219 *
## employ -243.62 3040.81 -0.080 0.9364
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 277400 on 58 degrees of freedom
## Multiple R-squared: 0.209, Adjusted R-squared: 0.1408
## F-statistic: 3.065 on 5 and 58 DF, p-value: 0.016- Andiamo a vedere le statistiche descrittive dei valori fitted del modello ausiliario
summary(aux_companies$fitted.values)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -80634 114204 156567 179298 211041 556892boxplot(aux_companies$fitted.values)
Vediamo che aux predice valori negativi dei residui al quadrato (quindi non è un buon modello perchè i residui al quadrato sono di per se positivi) Rispecifichiamo quindi aux_companies applicando ai residui al quadrato la funzione logaritmo, così da avere risultati che sono sempre positivi
aux_companies<-lm(I(log(mod1_companies$residuals^2))~assets+sales+profits+cash+employ,companies)
summary(aux_companies)##
## Call:
## lm(formula = I(log(mod1_companies$residuals^2)) ~ assets + sales +
## profits + cash + employ, data = companies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.6357 -0.6649 0.3790 1.2436 3.7077
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.027e+01 4.563e-01 22.508 <2e-16 ***
## assets -1.858e-04 9.420e-05 -1.973 0.0533 .
## sales 4.871e-05 3.298e-04 0.148 0.8831
## profits -5.172e-03 5.208e-03 -0.993 0.3248
## cash 6.335e-03 3.819e-03 1.659 0.1025
## employ 1.338e-02 2.160e-02 0.620 0.5379
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.97 on 58 degrees of freedom
## Multiple R-squared: 0.1926, Adjusted R-squared: 0.123
## F-statistic: 2.767 on 5 and 58 DF, p-value: 0.02617- Calcoliamo la standard deviation per ogni valore residuo, però prima devo ritrasformare dal logaritmo
companies$varianza<-exp(aux_companies$fitted.values)
companies$sd<-sqrt(companies$varianza)- Si costruisce il modello FGLS con le variabili trasformate e lo confrontiamo con quello precedente che aveva errori eteroschedastici
fgls1_companies<-lm(I(mark_val/sd)~0+I(1/sd)+I(assets/sd)+I(profits/sd)+I(sales/sd)+I(cash/sd)+I(employ/sd),companies)
summary(fgls1_companies)##
## Call:
## lm(formula = I(mark_val/sd) ~ 0 + I(1/sd) + I(assets/sd) + I(profits/sd) +
## I(sales/sd) + I(cash/sd) + I(employ/sd), data = companies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6940 -1.3034 -0.2434 1.1278 5.2908
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## I(1/sd) 243.716439 77.769231 3.134 0.00271 **
## I(assets/sd) -0.010719 0.009896 -1.083 0.28321
## I(profits/sd) 1.718887 1.339286 1.283 0.20444
## I(sales/sd) 0.051569 0.074081 0.696 0.48914
## I(cash/sd) 2.838162 1.067281 2.659 0.01011 *
## I(employ/sd) 4.518141 5.937544 0.761 0.44977
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.834 on 58 degrees of freedom
## Multiple R-squared: 0.8922, Adjusted R-squared: 0.8811
## F-statistic: 80.04 on 6 and 58 DF, p-value: < 2.2e-16summary(mod1_companies)##
## Call:
## lm(formula = mark_val ~ assets + sales + profits + cash + employ,
## data = companies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1228.23 -222.80 -50.87 251.70 1180.29
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.438e+02 1.030e+02 2.367 0.0213 *
## assets 2.066e-04 2.127e-02 0.010 0.9923
## sales -1.802e-02 7.445e-02 -0.242 0.8096
## profits 8.828e-02 1.176e+00 0.075 0.9404
## cash 3.787e+00 8.620e-01 4.393 4.82e-05 ***
## employ 9.399e+00 4.875e+00 1.928 0.0587 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 444.8 on 58 degrees of freedom
## Multiple R-squared: 0.7394, Adjusted R-squared: 0.7169
## F-statistic: 32.91 on 5 and 58 DF, p-value: 9.393e-16Vediamo che cash è rimasta significativa, anche se lo è molto meno, vediamo che il modello rimane significativo nel suo complesso e la bontà di afattamento del modello è cresciuta Andiamo a vedere se l’eteroschedasticità è stata risolta con il modello con stime FGLS
white.test(fgls1_companies)## Test.Statistic P
## 1 0.1750454 0.9161981Vediamo che il test di white conferma che sia stata eliminata l’omoschedasticità dei residui
Vi è anche un metodo alternativo per costruire gli stimatori FGLS (anche se il prof consiglia comunque di utilizzare il primo metodo)
fgls2_companies<-lm(mark_val~assets+sales+profits+cash+employ, weights=(1/companies$varianza),companies)
summary(fgls2_companies)##
## Call:
## lm(formula = mark_val ~ assets + sales + profits + cash + employ,
## data = companies, weights = (1/companies$varianza))
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -4.6940 -1.3034 -0.2434 1.1278 5.2908
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 243.716439 77.769231 3.134 0.00271 **
## assets -0.010719 0.009896 -1.083 0.28321
## sales 0.051569 0.074081 0.696 0.48914
## profits 1.718887 1.339286 1.283 0.20444
## cash 2.838162 1.067281 2.659 0.01011 *
## employ 4.518141 5.937544 0.761 0.44977
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.834 on 58 degrees of freedom
## Multiple R-squared: 0.6829, Adjusted R-squared: 0.6556
## F-statistic: 24.98 on 5 and 58 DF, p-value: 2.471e-13Vediamo che i risultati sono identici rispetto a fgls1 (l’unica cosa che cambia sono gli square)
Andiamo a vedere ora un metodo alternativo per affrontare l’eteroschedasticità (guardare slide prof), HC3 viene utilizzato quando il dataset contiene meno di 250 osservazioni
coeftest(mod1_companies, vcov=vcovHC(mod1_companies,type='HC3'))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.4378e+02 1.0965e+02 2.2232 0.03011 *
## assets 2.0661e-04 1.5334e-02 0.0135 0.98930
## sales -1.8019e-02 1.0774e-01 -0.1672 0.86776
## profits 8.8281e-02 1.9543e+00 0.0452 0.96413
## cash 3.7867e+00 1.4766e+00 2.5646 0.01294 *
## employ 9.3994e+00 7.8173e+00 1.2024 0.23410
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod1_companies)##
## Call:
## lm(formula = mark_val ~ assets + sales + profits + cash + employ,
## data = companies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1228.23 -222.80 -50.87 251.70 1180.29
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.438e+02 1.030e+02 2.367 0.0213 *
## assets 2.066e-04 2.127e-02 0.010 0.9923
## sales -1.802e-02 7.445e-02 -0.242 0.8096
## profits 8.828e-02 1.176e+00 0.075 0.9404
## cash 3.787e+00 8.620e-01 4.393 4.82e-05 ***
## employ 9.399e+00 4.875e+00 1.928 0.0587 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 444.8 on 58 degrees of freedom
## Multiple R-squared: 0.7394, Adjusted R-squared: 0.7169
## F-statistic: 32.91 on 5 and 58 DF, p-value: 9.393e-16Vediamo come i valori sono identici rispetto a quelli precedenti, ma gli std error sono più piccoli, quindi ci porteranno ad avere delle stime significative anche quando in realtà esse non lo sono
N.B: l’eteroschedasticità potrebbe anche essere dovuta ad una forma funzionale sbagliata, quindi un tentativo da fare potrebbe essere quello di provare varie forme funzionali, prima di ricorrere alle stime FGLS
Studio e soluzione dell’autocorrelazione dei residui
Andiamo inizialemente a costruire la nostra regressione e guardare com’è il nostro modello
mod1_hartangel<-lm(ftheft~partic+degrees+mtheft,hartangel)
summary(mod1_hartangel)##
## Call:
## lm(formula = ftheft ~ partic + degrees + mtheft, data = hartangel)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.8242 -3.2892 0.5658 2.5798 9.3956
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.06091 8.11357 -4.445 0.000111 ***
## partic 0.14143 0.02849 4.965 2.57e-05 ***
## degrees 0.53861 0.05380 10.012 4.45e-11 ***
## mtheft 0.05282 0.02279 2.318 0.027467 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.754 on 30 degrees of freedom
## Multiple R-squared: 0.9244, Adjusted R-squared: 0.9168
## F-statistic: 122.2 on 3 and 30 DF, p-value: < 2.2e-16Tutte le variabili sono significative, modello nel suo complesso è significativo, l’ R^2 è molto buono
Andiamo a vedere che gli errori sono omoschedastici
white_lm(mod1_hartangel)## # A tibble: 1 x 5
## statistic p.value parameter method alternative
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 4.88 0.559 6 White's Test greaterAndiamo a vedere un diagramma dei residui nel tempo
plot(hartangel$year,mod1_hartangel$residuals,xlab='year',ylab='residuals',type='b')
abline(h=0)
Vediamo che vi sono situazioni che confermano un’autocorrelazione positiva di ordine 1 dei residui
Andiamo a verificare questa evidenza grafica con i correlogrammi
par(mfrow=c(2,1))
acf(mod1_hartangel$residuals,main='autocorr. residuals')
pacf(mod1_hartangel$residuals,main='autocorr. parziale residuals')
Emerge l’autocorrelazione positiva di ordine uno sia per quanto riguarda le autocorrelazioni e le autocorrelazioni, mentre vediamo che anche l’autocorrelazione di quarto ordine vuole emergere, ma non supera il livello di significatività
Andiamo avedere con il test di darbinwatson se si ottengono ulteriori conferme riguardo la presenza di autocorrelazione (con max.lag=5 studiamo la presenza di autocorrelazione fino al quinto ordine)
durbinWatsonTest(mod1_hartangel, max.lag=5)## lag Autocorrelation D-W Statistic p-value
## 1 0.53952640 0.9051118 0.000
## 2 0.22105334 1.5208872 0.102
## 3 0.08943159 1.6523470 0.366
## 4 0.22930798 1.3375369 0.120
## 5 0.18871355 1.4154863 0.318
## Alternative hypothesis: rho[lag] != 0Vediamo che con il test di Durbin Watson risulta significativa l’autocorrelazione di primo ordine
3.Studio dell’ autocorrelazione Cerchiamo di risolvere il problema di autocorrelazione con la procedura di cochrane-orcutt, la quale si divide in diversi passi:
- Ottenere i residui di mod1_hartangel, per fare ciò andiamo a creare i residui ritardati attraverso la funzione datacombine, andiamo in fine ad andare a visualizzare i dati per verificare di avere realmente creato i residui ritardati
hartangel$u_hat<-mod1_hartangel$residuals
hartangel<-slide(data=hartangel,Var='u_hat',TimeVar='year',NewVar='u_hat_lag')##
## Lagging u_hat by 1 time units.View(hartangel) - regressione di uhat su u_hat_lag (ovvero regressione dei residui sui residui di lag1), in quanto vogliamo andare a ricercare il coefficente di autocorrelazione p
aux_hartangel<-lm(u_hat~u_hat_lag,hartangel)
summary(aux_hartangel) ##
## Call:
## lm(formula = u_hat ~ u_hat_lag, data = hartangel)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.9764 -2.3988 0.0085 2.3227 7.9746
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1119 0.6795 -0.165 0.87024
## u_hat_lag 0.5429 0.1503 3.612 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.903 on 31 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2961, Adjusted R-squared: 0.2734
## F-statistic: 13.04 on 1 and 31 DF, p-value: 0.001061Stimiamo il coefficente di autocorrelazione 0.5429
- Memorizziamo il coefficiente di autocorrelazione (che corrisponde al secondo elemento di aux_hartangel)
rho_hartangel<-aux_hartangel$coefficients[2]- Tramite la funzione slide andiamo a creare ftheft lag e anche tutte le altre variabili lag (ritardate)
hartangel<-slide(data=hartangel,Var='ftheft',TimeVar='year',NewVar='ftheft_lag')##
## Lagging ftheft by 1 time units.hartangel<-slide(data=hartangel,Var='mtheft',TimeVar='year',NewVar='mtheft_lag')##
## Lagging mtheft by 1 time units.hartangel<-slide(data=hartangel,Var='partic',TimeVar='year',NewVar='partic_lag')##
## Lagging partic by 1 time units.hartangel<-slide(data=hartangel,Var='degrees',TimeVar='year',NewVar='degrees_lag')##
## Lagging degrees by 1 time units.View(hartangel)- Andiamo a calcolare le variabili trasformate tramite il coefficente di autocorrelazione
hartangel$ftheft_t<-hartangel$ftheft-rho_hartangel*hartangel$ftheft_lag
hartangel$mtheft_t<-hartangel$mtheft-rho_hartangel*hartangel$mtheft_lag
hartangel$partic_t<-hartangel$partic-rho_hartangel*hartangel$partic_lag
hartangel$degrees_t<-hartangel$degrees-rho_hartangel*hartangel$degrees_lag
hartangel$interc_t<-1-rho_hartangel- Stimiamo il modello delle variabili trasformate (metto lo zero perchè devo togliere la costante del modello e al suo posto sarà inserito 1-rho) e lo confrontiamo con quello del modello con errori correlati
mod2_hartangel<-lm(ftheft_t~0+interc_t+partic_t+ degrees_t+mtheft_t,hartangel)
summary(mod2_hartangel)##
## Call:
## lm(formula = ftheft_t ~ 0 + interc_t + partic_t + degrees_t +
## mtheft_t, data = hartangel)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6504 -2.0545 -0.1914 2.0164 8.0233
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## interc_t -26.82350 10.27163 -2.611 0.01413 *
## partic_t 0.11814 0.03537 3.340 0.00231 **
## degrees_t 0.51676 0.06954 7.431 3.45e-08 ***
## mtheft_t 0.04257 0.03074 1.385 0.17672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.944 on 29 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9517, Adjusted R-squared: 0.945
## F-statistic: 142.8 on 4 and 29 DF, p-value: < 2.2e-16summary(mod1_hartangel)##
## Call:
## lm(formula = ftheft ~ partic + degrees + mtheft, data = hartangel)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.8242 -3.2892 0.5658 2.5798 9.3956
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.06091 8.11357 -4.445 0.000111 ***
## partic 0.14143 0.02849 4.965 2.57e-05 ***
## degrees 0.53861 0.05380 10.012 4.45e-11 ***
## mtheft 0.05282 0.02279 2.318 0.027467 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.754 on 30 degrees of freedom
## Multiple R-squared: 0.9244, Adjusted R-squared: 0.9168
## F-statistic: 122.2 on 3 and 30 DF, p-value: < 2.2e-16Vediamo che dopo aver corretto per l’autocorrelazione mtheft non è più significativo e il modello si adatta molto bene ai dati
Andiamo allora a vedere dal test di dw la presenza di autocorrelazione dei residui o meno a 5 lag
durbinWatsonTest(mod2_hartangel, max.lag=5)## lag Autocorrelation D-W Statistic p-value
## 1 0.09711383 1.795933 0.294
## 2 0.04064272 1.895891 0.710
## 3 -0.06849700 1.971350 0.936
## 4 0.21131125 1.358931 0.184
## 5 0.04234945 1.677867 0.816
## Alternative hypothesis: rho[lag] != 0Vediamo che il test conferma che abbiamo eliminato l’autocorrelazione dei residui
Controlliamo anche per via grafica di aver eliminato l’autocorrelazione dei residui
acf(mod2_hartangel$residuals,main='autocorr. residui di mod2')
pacf(mod2_hartangel$residuals,main='autocorr. parziale residui di mod 2')
vediamo anche come i correlogrammmi indicano che non vi sia più autocorrelazione
plot(hartangel$year[-1], mod2_hartangel$residuals,xlab='year',ylab='residui mod2',type='b')
abline(h=0)
Vediamo un secondo approccio, ovvero quello basato sul modello autoregressivo, order è un vettore con numero di AR, numero di differenziazioni, Numero di Ma. ML indica la massima verosomiglianza
#library('lmtest')
#mod3_hartangel<-arima(hartangel$ftheft, order=C(1,0,0),xreg=hartangel[,c('partic','degrees','mtheft')],method='ML')
#coeftest(mod3_hartangel) #abbiamo la significativit? di partic e degrees ma non pi? di mtheftStudio della normalità dei residui
Costruiamo la nostra regressione
m1_countries<-lm(Life.expectancy~ ISPs.million+ Irrigated+ Under.14 +Literacy.Rate,countries)
summary(m1_countries)##
## Call:
## lm(formula = Life.expectancy ~ ISPs.million + Irrigated + Under.14 +
## Literacy.Rate, data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.3828 -3.3426 0.5599 4.3173 17.5504
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.261e+01 1.462e+01 6.333 3.63e-07 ***
## ISPs.million -3.255e-02 1.861e-01 -0.175 0.862
## Irrigated 6.513e-06 1.493e-05 0.436 0.665
## Under.14 -9.566e-01 2.081e-01 -4.598 6.01e-05 ***
## Literacy.Rate 3.402e-02 1.125e-01 0.302 0.764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.867 on 33 degrees of freedom
## Multiple R-squared: 0.657, Adjusted R-squared: 0.6155
## F-statistic: 15.8 on 4 and 33 DF, p-value: 2.542e-07La normalità dei residui viene in un primo luogo indagata graficamente, attraverso l’uso di un QQ-plot
plot(m1_countries, which=2) # il prof consiglia di approfondire il which
abline(h=-2)
abline(h=2)
Vediamo in questo caso che vi sono valori che si discostano dalla normalità (probabilmente sono degli outlier), anche se gli errori sembrano approssimativamente distribuirsi in modo normale
Andiamo a verificare la normalità tramite l’uso di alcuni test:
ols_test_normality(m1_countries)## -----------------------------------------------
## Test Statistic pvalue
## -----------------------------------------------
## Shapiro-Wilk 0.9488 0.0811
## Kolmogorov-Smirnov 0.1373 0.4322
## Cramer-von Mises 2.6918 0.0000
## Anderson-Darling 0.7831 0.0384
## -----------------------------------------------I primi due test ci segnalano la presenza di normalità infatti accettano l’ipotesi nulla( N.B:il test di shapiro-wilk è molto sensibile a gli outliers)
Test F per restrizioni multiple sulle variabili singolarmente non significative
Andiamo ora a testare se vi siano delle restrizioni multiple, ovvero testare che siano uguali a zero i coefficienti delle variabili che non risultavano significative, ovvero urbana, alfabet e relig2 (nel dataset delle nazioni)
linearHypothesis(mod3_nazioni,c('urbana=0','relig2=0', 'alfabet=0'))## Linear hypothesis test
##
## Hypothesis:
## urbana = 0
## relig2 = 0
## alfabet = 0
##
## Model 1: restricted model
## Model 2: pil ~ densita + urbana + vitafem + alfabet + relig
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 55 943109645
## 2 52 922353689 3 20755956 0.3901 0.7606vediamo che questo test ci porta a non rigettare l’ipotesi nulla che i tre coefficienti siano uguali a zero
Creiamo quindi un nuovo dataset contenente solo le variabili ritenute significative
mod4_nazioni<-lm(pil~densita+vitafem, nazioni2)
summary(mod4_nazioni)##
## Call:
## lm(formula = pil ~ densita + vitafem, data = nazioni2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7403 -3964 -1016 3699 8994
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -26797.536 5368.682 -4.991 6.18e-06 ***
## densita 15.171 6.549 2.316 0.0242 *
## vitafem 436.936 71.089 6.146 8.78e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4889 on 56 degrees of freedom
## Multiple R-squared: 0.4148, Adjusted R-squared: 0.3939
## F-statistic: 19.85 on 2 and 56 DF, p-value: 3.046e-07MODELLO LINEARE CLASSICO MULTIVARIATO
Verifica delle assunzioni del modello lineare classico in ciascun modello (caso univariato)
Per la verifica delle assunzioni e di tutti i modi di risoluzione ci si rifà a quanto detto nella sezione precedente, adesso vedremo un breve codice che serve in maniera rapida ad avere uno sguardo su tutte le ipotesi del modello classico (per il primo modello)
par(mfrow=c(2,2))
k<-length(coef(m1_countries))
n<-length(m1_countries$fitted.values)
#residui studentizzati vs laverage
plot(hatvalues(m1_countries),rstudent(m1_countries),xlab='laverage',ylab='Stud. residuals')
abline(h=-2)
abline(h=2)
abline(h=0)
abline(v=2*k/n)
text(hatvalues(m1_countries),rstudent(m1_countries))
#andiamo ora a vedere la distanza di cook
plot(cooks.distance(m1_countries),type='h',ylab='Cook\'s Distance')
abline(h=4/n)
#eteroschedasticità: fitted vs residui standardizzati
plot(fitted.values(m1_countries),rstandard(m1_countries),xlba='fitted',ylab='std.residual')## Warning in plot.window(...): parametro grafico "xlba" non valido## Warning in plot.xy(xy, type, ...): parametro grafico "xlba" non valido## Warning in axis(side = side, at = at, labels = labels, ...): parametro grafico
## "xlba" non valido
## Warning in axis(side = side, at = at, labels = labels, ...): parametro grafico
## "xlba" non valido## Warning in box(...): parametro grafico "xlba" non valido## Warning in title(...): parametro grafico "xlba" non validoabline(h=-2)
abline(h=2)
abline(h=0)
#qq plot per la normalità
plot(m1_countries, which=2) # il prof consiglia di approfondire il which
abline(h=-2)
abline(h=2)
Notiamo che i residui standardizzati assumono una forma a ventaglio, quindi siamo in presenza di eteroschedasticità. Notiamo che le osservazioni 10, 34 e 9 (soprattutto) sono riconosciute come punti influenti, soprattutto l’osservazione 9, la quale emerge anche come osservazione influente per quanto riguarda la distanza di cool, mentre le osservazioni 32, 38 e 22 sono riconosciuti come valori outlier. Vediamo in questo caso che vi sono valori che si discostano dalla normalità (probabilmente sono degli outlier), anche se gli errori sembrano approssimativamente distribuirsi in modo normale.
Possiamo controllare quanto visto graficamente anche a livello di test:
#andremo a corroborare le seguenti idee con i relativi test
white.test(m1_countries)## Test.Statistic P
## 1 6.533431 0.03813147#è significativo quindi si può rifiutare l'ipotesi nulla di omoschedasticità
ols_test_normality(m1_countries)## -----------------------------------------------
## Test Statistic pvalue
## -----------------------------------------------
## Shapiro-Wilk 0.9488 0.0811
## Kolmogorov-Smirnov 0.1373 0.4322
## Cramer-von Mises 2.6918 0.0000
## Anderson-Darling 0.7831 0.0384
## -----------------------------------------------# i primi due test ci segnalano la presenza di normalità infatti accettano l'ipotesi nulla(N.B:il test di shapiro-wilk è molto sensibile a gli outliers)
vif(m1_countries)## ISPs.million Irrigated Under.14 Literacy.Rate
## 1.254729 1.045980 2.742422 2.583564#vediamo che non vi sono grossi problemi di muticollinearità, questo lo si poteva sapere anche dalla matrice di correlazioneDovremo togliere le osservazioni 9,10 e 34 (influenti) e gli outliers, dovremmmo poi risolvere gli errori omoschedastici con i Fgls (ma non sono l’obiettivo di questa esercitazione quiundi andiamo avanti)
Andiamo a dare uno sguardo alla possibile violazione delle ipotesi del modello lineare classico, anche per il secondo modello, cosi definito:
m2_countries<-lm(Unemployment~ ISPs.million+ Irrigated+ Under.14 +Literacy.Rate,countries)
summary(m2_countries)##
## Call:
## lm(formula = Unemployment ~ ISPs.million + Irrigated + Under.14 +
## Literacy.Rate, data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.3115 -5.2740 0.5191 5.1450 24.1047
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.725e+01 2.091e+01 -1.303 0.2015
## ISPs.million -3.646e-02 2.660e-01 -0.137 0.8918
## Irrigated -1.388e-05 2.134e-05 -0.650 0.5199
## Under.14 1.045e+00 2.974e-01 3.515 0.0013 **
## Literacy.Rate 1.497e-01 1.608e-01 0.931 0.3586
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.25 on 33 degrees of freedom
## Multiple R-squared: 0.4257, Adjusted R-squared: 0.3561
## F-statistic: 6.116 on 4 and 33 DF, p-value: 0.0008505Svolgiamo gli stessi passaggi che abbiamo svolto per il modello precedente
par(mfrow=c(2,2))
k<-length(coef(m2_countries))
n<-length(m2_countries$fitted.values)
#residui studentizzati vs laverage
plot(hatvalues(m2_countries),rstudent(m2_countries),xlab='laverage',ylab='Stud. residuals')
abline(h=-2)
abline(h=2)
abline(h=0)
abline(v=2*k/n)
text(hatvalues(m2_countries),rstudent(m2_countries))
#andiamo ora a vedere la distanza di cook
plot(cooks.distance(m2_countries),type='h',ylab='Cook\'s Distance')
abline(h=4/n)
#eteroschedasticità: fitted vs residui standardizzati
plot(fitted.values(m2_countries),rstandard(m2_countries),xlba='fitted',ylab='std.residual')## Warning in plot.window(...): parametro grafico "xlba" non valido## Warning in plot.xy(xy, type, ...): parametro grafico "xlba" non valido## Warning in axis(side = side, at = at, labels = labels, ...): parametro grafico
## "xlba" non valido
## Warning in axis(side = side, at = at, labels = labels, ...): parametro grafico
## "xlba" non valido## Warning in box(...): parametro grafico "xlba" non valido## Warning in title(...): parametro grafico "xlba" non validoabline(h=-2)
abline(h=2)
abline(h=0)
#qq plot per la normalità
plot(m2_countries, which=2) # il prof consiglia di approfondire il which
abline(h=-2)
abline(h=2)
#9,10,34 sono ancora dei valori influenti e abbiamo ancora 3 outlier che andrebbero tolti
#guardando la cook almeno l osservazione 9 deve venire tolta dal dataset
#sembra che abbiamo degli errori eteroschedastici
# e gli errori in questo caso sembrano distribuirsi normalmente
#andremo a corroborare le seguenti idee con i relativi test
white.test(m2_countries)## Test.Statistic P
## 1 12.93848 0.001550401#è significativo quindi si può rifiutare l'ipotesi nulla di omoschedasticità
ols_test_normality(m2_countries)## -----------------------------------------------
## Test Statistic pvalue
## -----------------------------------------------
## Shapiro-Wilk 0.9746 0.5291
## Kolmogorov-Smirnov 0.0976 0.8278
## Cramer-von Mises 2.587 0.0000
## Anderson-Darling 0.3648 0.4197
## -----------------------------------------------#vediamo che i test ci dicono (contro la nostra previsione grafica) che gli errori si distribuiscono normalmente
vif(m2_countries)## ISPs.million Irrigated Under.14 Literacy.Rate
## 1.254729 1.045980 2.742422 2.583564#vediamo che non vi sono grossi problemi di muticollinearità, questo lo si poteva sapere anche dalla matrice di correlazione
#dovremo togliere le osservazioni 9,10 e 34 (influenti) e gli outliers, dovremmmo poi risolvere gli errori eteroschedastici con i fgls (ma non sono l'obiettivo di questa esercitazione quindi andiamo avanti)Modello lineare multivariato e test d’ipotesi su restrizioni multiple
Iniziamo costruendo il modello multivariato
mv1_countries<-lm(cbind(Life.expectancy, Unemployment)~ISPs.million+ Irrigated+ Under.14 +Literacy.Rate,countries)
summary(mv1_countries)## Response Life.expectancy :
##
## Call:
## lm(formula = Life.expectancy ~ ISPs.million + Irrigated + Under.14 +
## Literacy.Rate, data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.3828 -3.3426 0.5599 4.3173 17.5504
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.261e+01 1.462e+01 6.333 3.63e-07 ***
## ISPs.million -3.255e-02 1.861e-01 -0.175 0.862
## Irrigated 6.513e-06 1.493e-05 0.436 0.665
## Under.14 -9.566e-01 2.081e-01 -4.598 6.01e-05 ***
## Literacy.Rate 3.402e-02 1.125e-01 0.302 0.764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.867 on 33 degrees of freedom
## Multiple R-squared: 0.657, Adjusted R-squared: 0.6155
## F-statistic: 15.8 on 4 and 33 DF, p-value: 2.542e-07
##
##
## Response Unemployment :
##
## Call:
## lm(formula = Unemployment ~ ISPs.million + Irrigated + Under.14 +
## Literacy.Rate, data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.3115 -5.2740 0.5191 5.1450 24.1047
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.725e+01 2.091e+01 -1.303 0.2015
## ISPs.million -3.646e-02 2.660e-01 -0.137 0.8918
## Irrigated -1.388e-05 2.134e-05 -0.650 0.5199
## Under.14 1.045e+00 2.974e-01 3.515 0.0013 **
## Literacy.Rate 1.497e-01 1.608e-01 0.931 0.3586
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.25 on 33 degrees of freedom
## Multiple R-squared: 0.4257, Adjusted R-squared: 0.3561
## F-statistic: 6.116 on 4 and 33 DF, p-value: 0.0008505Vediamo che otteniamo esattamente gli stessi risultati di guardare i due modelli univariati in modo separato
Proseguiamo ora svolgendo test d’ ipotesi su restrizioni multiple
Iniziamo stimando un modello nullo, ovvero un modello in cui compaia la sola costante e tutti i coefficenti siano nulli, in quanto vogliamo testare HO: tutti i coefficenti=0 in entrambe le equazioni
mvnull_countries<-lm(cbind(Life.expectancy, Unemployment)~1,countries)
anova(mv1_countries,mvnull_countries)## Analysis of Variance Table
##
## Model 1: cbind(Life.expectancy, Unemployment) ~ ISPs.million + Irrigated +
## Under.14 + Literacy.Rate
## Model 2: cbind(Life.expectancy, Unemployment) ~ 1
## Res.Df Df Gen.var. Pillai approx F num Df den Df Pr(>F)
## 1 33 64.344
## 2 37 4 102.953 0.74566 4.9043 8 66 8.865e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vediamo che la statistica F risulta significativa e ciò ci porta a rifiutare l’ipotesi nulla, quindi vi è almeno un coefficente che è significativo per almeno una variabile dipendente
Andiamo a prendere ora i coefficienti che non sono significativi in entrambe le equazioni e andiamo a testarli:
- prima ipotesi il coefficente Irrigated=0 in entrambe le equazioni
linearHypothesis(mv1_countries,hypothesis.matrix=c('Irrigated=0'))##
## Sum of squares and products for the hypothesis:
## Life.expectancy Unemployment
## Life.expectancy 11.78022 -25.10959
## Unemployment -25.10959 53.52120
##
## Sum of squares and products for error:
## Life.expectancy Unemployment
## Life.expectancy 2042.374 -2004.063
## Unemployment -2004.063 4174.011
##
## Multivariate Tests:
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0126660 0.2052553 2 32 0.8155
## Wilks 1 0.9873340 0.2052553 2 32 0.8155
## Hotelling-Lawley 1 0.0128285 0.2052553 2 32 0.8155
## Roy 1 0.0128285 0.2052553 2 32 0.8155Vediamo che non si rifiuta l’ipotesi nulla, ovvero che irrigated non sia significativo per entrambe le equazioni, quindi proseguiamo a rimuovere irrigated da entrambi i modelli
Costruiamo il nuovo modello di regressione multivariata escludendo irrigated
mv2_countries<-lm(cbind(Life.expectancy, Unemployment)~ISPs.million+ Under.14 +Literacy.Rate,countries)
summary(mv2_countries)## Response Life.expectancy :
##
## Call:
## lm(formula = Life.expectancy ~ ISPs.million + Under.14 + Literacy.Rate,
## data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.4030 -3.0755 0.4149 4.3508 17.5113
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94.07054 14.06813 6.687 1.12e-07 ***
## ISPs.million -0.02856 0.18360 -0.156 0.877
## Under.14 -0.97250 0.20237 -4.806 3.07e-05 ***
## Literacy.Rate 0.02451 0.10905 0.225 0.824
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.773 on 34 degrees of freedom
## Multiple R-squared: 0.655, Adjusted R-squared: 0.6246
## F-statistic: 21.52 on 3 and 34 DF, p-value: 5.412e-08
##
##
## Response Unemployment :
##
## Call:
## lm(formula = Unemployment ~ ISPs.million + Under.14 + Literacy.Rate,
## data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.8336 -4.9370 0.3016 4.5583 24.1477
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -30.34887 20.18195 -1.504 0.14188
## ISPs.million -0.04496 0.26340 -0.171 0.86548
## Under.14 1.07945 0.29032 3.718 0.00072 ***
## Literacy.Rate 0.17000 0.15644 1.087 0.28482
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.15 on 34 degrees of freedom
## Multiple R-squared: 0.4184, Adjusted R-squared: 0.3671
## F-statistic: 8.153 on 3 and 34 DF, p-value: 0.0003181Vediamo che la variabile ISPs.million non è significativa, andiamo a testare che il coefficente ISPs.million=0 in entrambe le equazioni
linearHypothesis(mv2_countries,hypothesis.matrix=c('ISPs.million=0'))##
## Sum of squares and products for the hypothesis:
## Life.expectancy Unemployment
## Life.expectancy 1.461820 2.301222
## Unemployment 2.301222 3.622621
##
## Sum of squares and products for error:
## Life.expectancy Unemployment
## Life.expectancy 2054.154 -2029.173
## Unemployment -2029.173 4227.532
##
## Multivariate Tests:
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0050029 0.0829631 2 33 0.92058
## Wilks 1 0.9949971 0.0829631 2 33 0.92058
## Hotelling-Lawley 1 0.0050281 0.0829631 2 33 0.92058
## Roy 1 0.0050281 0.0829631 2 33 0.92058Vediamo che non si rifiuta l’ipotesi nulla, ovvero che ISP.million non sia significativo per entrambe le equazioni, quindi proseguiamo a rimuovere ISP.million da entrambi i modelli
Costruiamo il nuovo modello di regressione multivariata escludendo ISP.million
mv3_countries<-lm(cbind(Life.expectancy, Unemployment)~ Under.14 +Literacy.Rate,countries)
summary(mv3_countries)## Response Life.expectancy :
##
## Call:
## lm(formula = Life.expectancy ~ Under.14 + Literacy.Rate, data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.3843 -2.8877 0.4936 4.4272 17.5317
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 93.80617 13.76903 6.813 6.67e-08 ***
## Under.14 -0.96430 0.19264 -5.006 1.58e-05 ***
## Literacy.Rate 0.02349 0.10732 0.219 0.828
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.664 on 35 degrees of freedom
## Multiple R-squared: 0.6548, Adjusted R-squared: 0.6351
## F-statistic: 33.2 on 2 and 35 DF, p-value: 8.245e-09
##
##
## Response Unemployment :
##
## Call:
## lm(formula = Unemployment ~ Under.14 + Literacy.Rate, data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.8864 -4.9807 0.4401 4.7070 24.1772
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -30.7650 19.7543 -1.557 0.128377
## Under.14 1.0924 0.2764 3.952 0.000358 ***
## Literacy.Rate 0.1684 0.1540 1.094 0.281580
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11 on 35 degrees of freedom
## Multiple R-squared: 0.4179, Adjusted R-squared: 0.3846
## F-statistic: 12.56 on 2 and 35 DF, p-value: 7.721e-05Vediamo ce il coefficente di Literacy.Rate non è significativo, andiamo a testare che il coefficente ISPs.million=0 in entrambe le equazioni
linearHypothesis(mv3_countries,hypothesis.matrix=c('Literacy.Rate=0'))##
## Sum of squares and products for the hypothesis:
## Life.expectancy Unemployment
## Life.expectancy 2.813344 20.16898
## Unemployment 20.168980 144.59227
##
## Sum of squares and products for error:
## Life.expectancy Unemployment
## Life.expectancy 2055.616 -2026.872
## Unemployment -2026.872 4231.155
##
## Multivariate Tests:
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0784870 1.447923 2 34 0.24919
## Wilks 1 0.9215130 1.447923 2 34 0.24919
## Hotelling-Lawley 1 0.0851719 1.447923 2 34 0.24919
## Roy 1 0.0851719 1.447923 2 34 0.24919Vediamo che non si rifiuta l’ipotesi nulla, ovvero che Literacy.Rate non sia significativo per entrambe le equazioni, quindi proseguiamo a rimuovere Literacy.Rate da entrambi i modelli
Costruiamo il nuovo modello di regressione multivariata escludendo Literacy.Rate
mv4_countries<-lm(cbind(Life.expectancy, Unemployment)~ Under.14,countries)
summary(mv4_countries)## Response Life.expectancy :
##
## Call:
## lm(formula = Life.expectancy ~ Under.14, data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.0250 -2.9740 0.4021 4.7701 17.8797
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 96.7036 3.7354 25.888 <2e-16 ***
## Under.14 -0.9969 0.1208 -8.255 8e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.562 on 36 degrees of freedom
## Multiple R-squared: 0.6543, Adjusted R-squared: 0.6447
## F-statistic: 68.15 on 1 and 36 DF, p-value: 7.998e-10
##
##
## Response Unemployment :
##
## Call:
## lm(formula = Unemployment ~ Under.14, data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.0690 -5.1599 -0.6156 4.7316 26.7526
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.9935 5.4463 -1.835 0.0748 .
## Under.14 0.8589 0.1761 4.878 2.18e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.02 on 36 degrees of freedom
## Multiple R-squared: 0.398, Adjusted R-squared: 0.3813
## F-statistic: 23.8 on 1 and 36 DF, p-value: 2.177e-05Andiamo adesso a vedere cos’è successo nei vari modelli per quanto riguarda la bontà di adattamento
summary(mv1_countries)$`Response Life.expectancy`$adj.r.squared## [1] 0.6154522summary(mv2_countries)$`Response Life.expectancy`$adj.r.squared## [1] 0.6246097summary(mv3_countries)$`Response Life.expectancy`$adj.r.squared## [1] 0.6350756summary(mv4_countries)$`Response Life.expectancy`$adj.r.squared## [1] 0.6447268summary(mv1_countries)$`Response Unemployment`$adj.r.squared## [1] 0.3561386summary(mv2_countries)$`Response Unemployment`$adj.r.squared## [1] 0.3670626summary(mv3_countries)$`Response Unemployment`$adj.r.squared## [1] 0.3846196summary(mv4_countries)$`Response Unemployment`$adj.r.squared## [1] 0.3812681Notiamo che avendo tolto delle variabili non significative il nostro indice di determinazione aggiustato (che tiene conto del numero delle variabili esplicative incluse nel modello) è aumentato
MODELLO SUR
Consegna: Regressioni lineari multiple di packpc sulle covariate cpi, pop, income e tax per gli stati Arkansas (state=”AR”) e California (state=”CA”)
N.B: la fase di preprocessing con subsetting è stata già effettuata ed esposta nella sezione di preprocessing
Andiamo adesso a stimare i modelli di regressione, questi due modelli hanno variabili diverse (es. income viene trattato diversamente per ognun dei due stati, l’unica variabile in comune è cpi)
ar1<-lm(packpc_AR~cpi+pop_AR+income_AR+tax_AR,cigarettesarca)
summary(ar1)##
## Call:
## lm(formula = packpc_AR ~ cpi + pop_AR + income_AR + tax_AR, data = cigarettesarca)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1297 -0.9481 -0.4056 0.5941 3.4628
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.195e+02 8.439e+01 4.971 0.00252 **
## cpi -2.577e+02 6.589e+01 -3.912 0.00788 **
## pop_AR -6.305e-05 4.495e-05 -1.403 0.21032
## income_AR 6.858e-06 1.869e-06 3.670 0.01046 *
## tax_AR -1.333e+00 5.979e-01 -2.229 0.06735 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.311 on 6 degrees of freedom
## Multiple R-squared: 0.9563, Adjusted R-squared: 0.9272
## F-statistic: 32.82 on 4 and 6 DF, p-value: 0.0003231Vediamo un F significativo un R^2 che spiega molto bene i dati e le variabili significative sono income e cpi pairs.panels(d1) #vediamo che ci sono delle correlazioni molto elevate, queste dipendono sia al campione molto rispetto sia alla
ca1<-lm(packpc_CA~cpi+pop_CA+income_CA+tax_CA,cigarettesarca)
summary(ca1)##
## Call:
## lm(formula = packpc_CA ~ cpi + pop_CA + income_CA + tax_CA, data = cigarettesarca)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99694 -0.33594 0.07137 0.30917 1.27674
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.788e+02 3.336e+01 8.358 0.000159 ***
## cpi -6.138e+01 1.290e+01 -4.756 0.003138 **
## pop_CA -4.806e-06 1.593e-06 -3.016 0.023506 *
## income_CA 3.840e-08 3.305e-08 1.162 0.289449
## tax_CA -1.120e-01 7.017e-02 -1.596 0.161627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7991 on 6 degrees of freedom
## Multiple R-squared: 0.9985, Adjusted R-squared: 0.9974
## F-statistic: 976.3 on 4 and 6 DF, p-value: 1.442e-08Vediamo che F è significativo, buona capacità di adattamento ai dati, variabili significative cpi e pop_ca
Andiamo a vedere le correlazioni tra i due modelli
pairs.panels(cigarettesarca)
Notiamo che ci sono delle correlazioni molto grandi, questo è dovuto al fatto che le variabili sono state costruite dalla stessa variabile
Ora dovremmo andare a vedere le ipotesi del modello lineare classico per ognuno dei due modelli (lo facciamo solo per uno in maniera del tutto didattica anche perchè non è lo scopo di questa parte)
vif(ar1)## cpi pop_AR income_AR tax_AR
## 206.94592 17.84578 298.48359 38.50857#vediamo che ci sono dei problemi di multicollinearità molto importanti
white.test(ar1)## Test.Statistic P
## 1 2.019625 0.3642873#abbiamo errori omoschedastici
ols_test_normality(ar1)## -----------------------------------------------
## Test Statistic pvalue
## -----------------------------------------------
## Shapiro-Wilk 0.9071 0.2253
## Kolmogorov-Smirnov 0.2383 0.4870
## Cramer-von Mises 1.4848 1e-04
## Anderson-Darling 0.4575 0.2121
## -----------------------------------------------#vediamo che abbiamo degli errori normaliAndiamo ora a vedere la correlazione tra i residui di ar1 e ca1
cor(resid(ar1),resid(ca1))## [1] 0.5666212Vediamo che c’è una correlazione del 56% dei residui che non è poco, quindi possiamo utilizzare un modello apposito che accomoda questa condizione di correlazione dei residui, ovvero la integra all’interno del modello stesso. Ovvero il modello SUR
Andiamo quindi a costruire un modello SUR
Dobbiamo inanzitutto memorizzare le nostre regressioni all’interno di una lista
elemento1<-packpc_AR~cpi+pop_AR+income_AR+tax_AR
elemento2<-packpc_CA~cpi+pop_CA+income_CA+tax_CA
sistema<-list(e1=elemento1,e2=elemento2)Entriamo nel vivo del modello sure tramite la libreria systemfit
sur_cigarettes<-systemfit(sistema,'SUR',data=cigarettesarca)
summary(sur_cigarettes)##
## systemfit results
## method: SUR
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 22 12 38.9895 1.24531 0.987932 0.998651
##
## N DF SSR MSE RMSE R2 Adj R2
## e1 11 6 34.79245 5.798741 2.408057 0.95255 0.920916
## e2 11 6 4.19707 0.699512 0.836368 0.99832 0.997199
##
## The covariance matrix of the residuals used for estimation
## e1 e2
## e1 5.34166 1.046491
## e2 1.04649 0.638571
##
## The covariance matrix of the residuals
## e1 e2
## e1 5.79874 1.676596
## e2 1.67660 0.699512
##
## The correlations of the residuals
## e1 e2
## e1 1.000000 0.832461
## e2 0.832461 1.000000
##
##
## SUR estimates for 'e1' (equation 1)
## Model Formula: packpc_AR ~ cpi + pop_AR + income_AR + tax_AR
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.63192e+02 7.08173e+01 5.12857 0.0021600 **
## cpi -2.49217e+02 5.84900e+01 -4.26084 0.0053163 **
## pop_AR -3.74634e-05 3.75376e-05 -0.99803 0.3567999
## income_AR 6.52862e-06 1.65535e-06 3.94394 0.0075902 **
## tax_AR -1.43426e+00 5.24627e-01 -2.73386 0.0340115 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.408057 on 6 degrees of freedom
## Number of observations: 11 Degrees of Freedom: 6
## SSR: 34.792445 MSE: 5.798741 Root MSE: 2.408057
## Multiple R-Squared: 0.95255 Adjusted R-Squared: 0.920916
##
##
## SUR estimates for 'e2' (equation 2)
## Model Formula: packpc_CA ~ cpi + pop_CA + income_CA + tax_CA
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.88314e+02 2.80560e+01 10.27636 4.9571e-05 ***
## cpi -6.77451e+01 1.20047e+01 -5.64320 0.0013272 **
## pop_CA -5.20450e-06 1.36478e-06 -3.81344 0.0088280 **
## income_CA 5.63196e-08 2.82632e-08 1.99268 0.0933686 .
## tax_CA -1.27173e-01 6.12950e-02 -2.07477 0.0833353 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.836368 on 6 degrees of freedom
## Number of observations: 11 Degrees of Freedom: 6
## SSR: 4.197071 MSE: 0.699512 Root MSE: 0.836368
## Multiple R-Squared: 0.99832 Adjusted R-Squared: 0.997199Vediamo che abbiamo degli R^2 molto simili alle stime precedenti e aumentano le variabili che sono significative, mentre i coefficienti non sono cambiati di molto rispetto a quelli precedenti
##Testare l’uguaglianza dI coefficentre tra due equazioni
Andiamo a testare l’uguaglianza del coefficiente cpi tra le due equazioni
linearHypothesis(sur_cigarettes,'e1_cpi=e2_cpi')## Linear hypothesis test (Theil's F test)
##
## Hypothesis:
## e1_cpi - e2_cpi = 0
##
## Model 1: restricted model
## Model 2: sur_cigarettes
##
## Res.Df Df F Pr(>F)
## 1 13
## 2 12 1 10.968 0.006204 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vediamo che il cpi è seignificativo, quindi rifiutiamo l’ipotesi nulla di uguaglianza di cpi tra le due equazioni
vi è un modo alternativo, uguale, per testare l’ipotesi nulla di uguaglianza dei coefficienti e si articola nei seguenti passaggi:
-estraggo i coefficienti del modello sur
coef(sur_cigarettes)## e1_(Intercept) e1_cpi e1_pop_AR e1_income_AR e1_tax_AR
## 3.631918e+02 -2.492166e+02 -3.746342e-05 6.528615e-06 -1.434257e+00
## e2_(Intercept) e2_cpi e2_pop_CA e2_income_CA e2_tax_CA
## 2.883141e+02 -6.774506e+01 -5.204497e-06 5.631957e-08 -1.271733e-01-vediamo che coef(m1) è di 10 elementi quindi possiamo anche creare una matrice di ipotesi di 10 elementi
lh1_cigarettes<-matrix(0,nrow=1,ncol=10) #creiamo una matrice di 0 di dimensioni 1x10-andiamo adesso a localizzare cpi cehe è per e1 al secondo elemento e per e2 al settimo elemento
lh1_cigarettes[2]<-1
lh1_cigarettes[7]<-1-andiamo a verificare l’ipotesi nulla
linearHypothesis(sur_cigarettes,lh1_cigarettes) #identico al linearHypothesis fatto in precedenza## Linear hypothesis test (Theil's F test)
##
## Hypothesis:
## e1_cpi + e2_cpi = 0
##
## Model 1: restricted model
## Model 2: sur_cigarettes
##
## Res.Df Df F Pr(>F)
## 1 13
## 2 12 1 32.961 9.291e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Consegna: Rimuovere variabili non significative dalla regressione al punto 3 e testare l’uguaglianza dei coefficienti di cpi con quelli ottenuti al punto 1.
Procediamo rivedendo quali erano le variabili non significative
summary(sur_cigarettes)##
## systemfit results
## method: SUR
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 22 12 38.9895 1.24531 0.987932 0.998651
##
## N DF SSR MSE RMSE R2 Adj R2
## e1 11 6 34.79245 5.798741 2.408057 0.95255 0.920916
## e2 11 6 4.19707 0.699512 0.836368 0.99832 0.997199
##
## The covariance matrix of the residuals used for estimation
## e1 e2
## e1 5.34166 1.046491
## e2 1.04649 0.638571
##
## The covariance matrix of the residuals
## e1 e2
## e1 5.79874 1.676596
## e2 1.67660 0.699512
##
## The correlations of the residuals
## e1 e2
## e1 1.000000 0.832461
## e2 0.832461 1.000000
##
##
## SUR estimates for 'e1' (equation 1)
## Model Formula: packpc_AR ~ cpi + pop_AR + income_AR + tax_AR
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.63192e+02 7.08173e+01 5.12857 0.0021600 **
## cpi -2.49217e+02 5.84900e+01 -4.26084 0.0053163 **
## pop_AR -3.74634e-05 3.75376e-05 -0.99803 0.3567999
## income_AR 6.52862e-06 1.65535e-06 3.94394 0.0075902 **
## tax_AR -1.43426e+00 5.24627e-01 -2.73386 0.0340115 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.408057 on 6 degrees of freedom
## Number of observations: 11 Degrees of Freedom: 6
## SSR: 34.792445 MSE: 5.798741 Root MSE: 2.408057
## Multiple R-Squared: 0.95255 Adjusted R-Squared: 0.920916
##
##
## SUR estimates for 'e2' (equation 2)
## Model Formula: packpc_CA ~ cpi + pop_CA + income_CA + tax_CA
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.88314e+02 2.80560e+01 10.27636 4.9571e-05 ***
## cpi -6.77451e+01 1.20047e+01 -5.64320 0.0013272 **
## pop_CA -5.20450e-06 1.36478e-06 -3.81344 0.0088280 **
## income_CA 5.63196e-08 2.82632e-08 1.99268 0.0933686 .
## tax_CA -1.27173e-01 6.12950e-02 -2.07477 0.0833353 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.836368 on 6 degrees of freedom
## Number of observations: 11 Degrees of Freedom: 6
## SSR: 4.197071 MSE: 0.699512 Root MSE: 0.836368
## Multiple R-Squared: 0.99832 Adjusted R-Squared: 0.997199Inseriamo le nostre nuove regressioni all’interno di due elementi
e1_cigarettes2<-packpc_AR~cpi+income_AR+tax_AR
e2_cigarettes2<-packpc_CA~cpi+pop_CA
sistema<-list(e1=e1_cigarettes2,e2=e2_cigarettes2)
m2_cigarettes<-systemfit(sistema,'SUR',data=cigarettesarca)
summary(m2_cigarettes)##
## systemfit results
## method: SUR
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 22 15 50.5545 4.65182 0.984353 0.995121
##
## N DF SSR MSE RMSE R2 Adj R2
## e1 11 7 43.73176 6.247395 2.499479 0.940358 0.914798
## e2 11 8 6.82269 0.852837 0.923492 0.997268 0.996585
##
## The covariance matrix of the residuals used for estimation
## e1 e2
## e1 6.079631 0.469357
## e2 0.469357 0.847833
##
## The covariance matrix of the residuals
## e1 e2
## e1 6.247395 0.822304
## e2 0.822304 0.852837
##
## The correlations of the residuals
## e1 e2
## e1 1.000000 0.356246
## e2 0.356246 1.000000
##
##
## SUR estimates for 'e1' (equation 1)
## Model Formula: packpc_AR ~ cpi + income_AR + tax_AR
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.95198e+02 3.29457e+01 8.96016 4.3905e-05 ***
## cpi -2.63278e+02 6.39618e+01 -4.11618 0.0044810 **
## income_AR 6.91440e-06 1.83716e-06 3.76363 0.0070413 **
## tax_AR -1.83126e+00 3.93425e-01 -4.65465 0.0023291 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.499479 on 7 degrees of freedom
## Number of observations: 11 Degrees of Freedom: 7
## SSR: 43.731764 MSE: 6.247395 Root MSE: 2.499479
## Multiple R-Squared: 0.940358 Adjusted R-Squared: 0.914798
##
##
## SUR estimates for 'e2' (equation 2)
## Model Formula: packpc_CA ~ cpi + pop_CA
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.52756e+02 1.34815e+01 18.74839 6.7657e-08 ***
## cpi -6.36196e+01 9.81623e+00 -6.48106 0.0001919 ***
## pop_CA -3.17627e-06 8.74255e-07 -3.63311 0.0066550 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.923492 on 8 degrees of freedom
## Number of observations: 11 Degrees of Freedom: 8
## SSR: 6.822694 MSE: 0.852837 Root MSE: 0.923492
## Multiple R-Squared: 0.997268 Adjusted R-Squared: 0.996585Andiamo a testare l’uguaglianza dei coefficienti di cpi del modello sure con quelli ottenuti al punto 1
coef(ar1)## (Intercept) cpi pop_AR income_AR tax_AR
## 4.194824e+02 -2.577434e+02 -6.304494e-05 6.858242e-06 -1.332888e+00coef(ar1)[2]## cpi
## -257.7434coef(ca1)[2]## cpi
## -61.37899Andiamo quindi ora a testare l’ipotesi nulla di uguaglianza dei coefficienti
linearHypothesis(sur_cigarettes,c('e1_cpi=-257.7434','e2_cpi=-61.37899 '))## Linear hypothesis test (Theil's F test)
##
## Hypothesis:
## e1_cpi = - 257.7434
## e2_cpi = - 61.37899
##
## Model 1: restricted model
## Model 2: sur_cigarettes
##
## Res.Df Df F Pr(>F)
## 1 14
## 2 12 2 0.18 0.8375Non rifiutiamo l’ipotesi nulla che i coefficenti stimati tramite il modello sur siano uguali a quelli inizialmente stimati tramite il modello ols
Le stime sure saranno tanto più efficenti di quelle OLS tanto più saranno alte le correlazioni tra i residui delle due equazioni
Modello sur in hsb
Costruamo le due regressioni di interesse
m1_hsb<-lm(science~math+prog+gender,hsb)
summary(m1_hsb)##
## Call:
## lm(formula = science ~ math + prog + gender, data = hsb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.6780 -3.9337 -0.1527 4.2994 22.8762
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.39716 3.79110 3.534 0.000511 ***
## math 0.69548 0.06534 10.645 < 2e-16 ***
## proggeneral 3.27165 1.41948 2.305 0.022229 *
## progvocational 0.56657 1.46593 0.386 0.699556
## gendermale 2.11313 1.07664 1.963 0.051103 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.575 on 195 degrees of freedom
## Multiple R-squared: 0.4264, Adjusted R-squared: 0.4147
## F-statistic: 36.25 on 4 and 195 DF, p-value: < 2.2e-16m2_hsb<-lm(write~read+prog+gender,hsb)
summary(m2_hsb)##
## Call:
## lm(formula = write ~ read + prog + gender, data = hsb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.1576 -4.3184 0.3107 4.6172 15.9053
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.06767 3.06588 10.133 < 2e-16 ***
## read 0.49127 0.05316 9.241 < 2e-16 ***
## proggeneral -1.67434 1.28609 -1.302 0.194491
## progvocational -4.53669 1.30794 -3.469 0.000644 ***
## gendermale -5.36496 0.99021 -5.418 1.76e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.958 on 195 degrees of freedom
## Multiple R-squared: 0.472, Adjusted R-squared: 0.4612
## F-statistic: 43.58 on 4 and 195 DF, p-value: < 2.2e-16Andiamo ad indagare la correlazione tra i residui delle due equazioni (due modelli)
cor(resid(m1_hsb),resid(m2_hsb))## [1] 0.1773396cor(hsb$write,hsb$science)## [1] 0.5704416Vediamo che tra i residui dei due modelli vi è una correlazione del 17% e tra le variabili dipendenti delle due equazioni vi è una correlazione del 57%, quindi si confutano le ipotesi alla base del modello lineare multivariato classico, andiamo quindi a creare un modello SUR
e1_hsb<-science~math+prog+gender
e2_hsb<-write~read+prog+gender
syst<-list(e1=e1_hsb,e2=e2_hsb)
sur_hsb<-systemfit(syst,'SUR',data=hsb)
summary(sur_hsb)##
## systemfit results
## method: SUR
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 400 390 20700.9 2668.83 0.446297 0.400312
##
## N DF SSR MSE RMSE R2 Adj R2
## e1 200 195 11225.11 57.5647 7.58714 0.424575 0.412771
## e2 200 195 9475.83 48.5940 6.97094 0.469998 0.459126
##
## The covariance matrix of the residuals used for estimation
## e1 e2
## e1 57.37741 9.34646
## e2 9.34646 48.41076
##
## The covariance matrix of the residuals
## e1 e2
## e1 57.5647 11.3345
## e2 11.3345 48.5940
##
## The correlations of the residuals
## e1 e2
## e1 1.000000 0.214306
## e2 0.214306 1.000000
##
##
## SUR estimates for 'e1' (equation 1)
## Model Formula: science ~ math + prog + gender
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.3386075 3.7529222 4.35357 2.1607e-05 ***
## math 0.6433571 0.0646401 9.95291 < 2.22e-16 ***
## proggeneral 2.9211667 1.4180387 2.06000 0.040726 *
## progvocational 0.0285594 1.4626432 0.01953 0.984442
## gendermale 2.1482484 1.0766251 1.99535 0.047396 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.587139 on 195 degrees of freedom
## Number of observations: 200 Degrees of Freedom: 195
## SSR: 11225.112722 MSE: 57.564681 Root MSE: 7.587139
## Multiple R-Squared: 0.424575 Adjusted R-Squared: 0.412771
##
##
## SUR estimates for 'e2' (equation 2)
## Model Formula: write ~ read + prog + gender
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.6081188 3.0355240 11.07160 < 2.22e-16 ***
## read 0.4456005 0.0525946 8.47237 5.7732e-15 ***
## proggeneral -1.9679990 1.2851261 -1.53137 0.12729952
## progvocational -4.9923735 1.3056647 -3.82363 0.00017687 ***
## gendermale -5.3097555 0.9901616 -5.36251 2.3082e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.970941 on 195 degrees of freedom
## Number of observations: 200 Degrees of Freedom: 195
## SSR: 9475.833735 MSE: 48.594019 Root MSE: 6.970941
## Multiple R-Squared: 0.469998 Adjusted R-Squared: 0.459126Andiamo a vedere che l’ R^2 sono rimasti più o meno come erano prima, vediamo che rispetto a m1 non ci sono grandi cambiamenti nei coefficenti ed emerge la significatività anche del sesso maschile, invece m1 sembra rimasto uguale
N.B:ricordiamo che la lettura delle variabili dummy avviene per esempio in questo modo, chi è maschio (gendermale) ha voti più bassi in lettura rispoetto alle femmine a parita di altre condizioni(le femmine è la variabile che non compare), stessa cosa per le variabili dummy a 3 variabili
Test significatività di prog in entrambe le equazioni, testiamo quindi prog=0 ovvero testare che non vi sia differenza tra i tipi di istruzione scelta rispetto alle variabili dipendenti
linearHypothesis(sur_hsb,c('e1_proggeneral=0','e2_proggeneral=0','e1_progvocational=0','e2_progvocational=0'))## Linear hypothesis test (Theil's F test)
##
## Hypothesis:
## e1_proggeneral = 0
## e2_proggeneral = 0
## e1_progvocational = 0
## e2_progvocational = 0
##
## Model 1: restricted model
## Model 2: sur_hsb
##
## Res.Df Df F Pr(>F)
## 1 394
## 2 390 4 5.0682 0.0005435 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vediamo ch questo test è significativo, ovvero esiste almeno un coefficente che è diverso da 0
Test di uguaglianza dei coefficienti di gender tra equazioni
linearHypothesis(sur_hsb,c('e1_gendermale=e2_gendermale'))## Linear hypothesis test (Theil's F test)
##
## Hypothesis:
## e1_gendermale - e2_gendermale = 0
##
## Model 1: restricted model
## Model 2: sur_hsb
##
## Res.Df Df F Pr(>F)
## 1 391
## 2 390 1 31.672 3.493e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vediamo come atteso che sono differenti