Transcript Scarica file

energia<-read.table("C:\\Users\\Public\\Documents\\DOCrosaria\\lezioni\\Data
Mining\\corso2013\\dati-energia.txt",header=T)
energia<-read.table("dati-energia.txt",header=T)
energia
dim(energia)
attach(energia)
> names(energia)
[1] "Socie" "Poteninst" "NImpianti" "Ndip"
"Fatturato"
[6] "Produzione" "Ricavo" "Provincia" "ConvBanche" "Fgiur"
Pacchetti o librerie da caricare per la REGRESSIONE
1# pacchetto {PERregress}
2# pacchetto {tseries}
3# pacchetto{moments}
4# pacchetto{lmtest} #opzionale
5# pacchetto{leaps}
6#pacchetto{car}
boxplot(Fatturato)
regressione1<- lm(Fatturato~ NImpianti+Produzione+Ndip+Ricavo+Poteninst)
summary(regressione1)
plot (Fatturato ~ Ndip + Produzione)
plot (Fatturato ~ Ndip)# scatter-plot regressione semplice
plot (Fatturato ~ Produzione) # scatter-plot regressione semplice
abline (regressione1$coeff[1],regressione1$coeff[3]) #tracciamo la retta della regressione
boxplot(regressione1$residuals)
hist(regressione1$residuals)
plot(density(regressione1$residuals))
plot(fitted(regressione1), residuals(regressione1))
stdResiduals1<-residuals(regressione1)/sqrt(var(residuals(regressione1)))
plot(fitted(regressione1), stdResiduals1)
sort(stdResiduals1)
VERIFICA NORMALITA’
normPlot (regressione1$residuals )
# pacchetto {PERregress}
jarque.bera.test (stdResiduals1)
# pacchetto {tseries}
skewness (stdResiduals1)
# pacchetto{moments}
kurtosis (stdResiduals1)
# pacchetto{moments}
VERIFICA OMOSCHEDASTICITA’
ncvTest(regressione1) # pacchetto{car} equivale a
e graficamente per visualizzare lo “spread”
spreadLevelPlot(regressione1)
OPPURE
spl(regressione1)
test di Breusch-Pagan
Ipotesi di omoschedasticità (equivale a ncvtest )
library(lmtest)
Il test di Breusch-Pagan, largamente utilizzato in econometria per verificare l’ipotesi di omoschedasticità,
applica ai residui gli stessi concetti della regressione lineare.
Per effettuare il test di Breusch Pagan in R è necessario caricare il package lmtest ove è contenuta la
funzione bptest().
help("bptest")
bptest(regressione1)
bptest(stdResiduals1~1)
Ipotesi di Multicollinearità per la funzione VIF o copiare il seguente codice
oppure caricare la libreria "car" (companion to Apllied Regression) dal CRAN di R
library(car)
help("vif")
vif(regressione1)
vif(regressione1) # variance inflation factors
sqrt(vif(regressione1)) > 2 # problemi
> x= NImpianti + Produzione + Ndip + Ricavo
VERIFICA LINEARITA’
crPlots(regressione1)
VALORI ANOMALI ED INFLUENTI
influence.measures(regressione1)
si calcolano per ciascuna variabile alcune statistiche “diagnostiche DFBETAS, DFFITS,
covariance ratios, Cook's distances e gli elementi diagonali della matrice di proiezionehat matrix-.
I casi influenti sono marcati con un asterisco.
Grafico:
plot(hatvalues(lm.out))
AUTOMATICA SELEZIONE DEL MIGLIOR MODELLO DI REGRESSIONE MULTIPLA
library(leaps)
# The leaps function provides the all-subsets analysis:
# Printing the 2 best models of each size, using the Cp criterion:
leaps(x = cbind(Produzione, NumImpianti, Numdipendenti), y = Fatturato, nbest=2, method="Cp")
# Printing the 2 best models of each size, using the adjusted R^2 criterion:
leaps(x
=
cbind(Produzione,
nbest=2,method="adjr2")
NumImpianti,
Numdipendenti),
y
=
Fatturato,
# The TRUE and FALSE values in the output refer to which variables are in each
# candidate model and which variables are out.
################################################################################
#####**
Per esempio con tre variabili indipendenti si ottengono 5 modelli nel primo è
presente solo la seconda variabile, nel secondo la prima variabile, nel terzo le
prime due variabili, nel quarto modello la prima e la terza, e infine nel quinto
modello tutte e tre le variabili. Il valore di R2 corretto aiuta nella scelta
del modello
leaps(x = cbind(Produzione, NumImpianti, Numdipendenti), y = Fatturato, nbest=2,
method="adjr2")
$which
1 2 3
1 FALSE TRUE FALSE
1 TRUE FALSE FALSE
2 TRUE TRUE FALSE
2 TRUE FALSE TRUE
3 TRUE TRUE TRUE
$label
[1] "(Intercept)" "1"
"2"
"3"
$size
[1] 2 2 3 3 4
$adjr2
[1] 0.8607942 0.7661865 0.9183009 0.8927873 0.9346892
help("leaps")
Per esempio con due variabili indipendenti
leaps(x = cbind(Produzione, NumImpianti), y = Fatturato, nbest=2, method="adjr2")
$which
1 2
1 FALSE TRUE
1 TRUE FALSE
2 TRUE TRUE
$label
[1] "(Intercept)" "1"
"2"
$size
[1] 2 2 3
$adjr2
[1] 0.8607942 0.7661865 0.9183009