PROJEKT REGRESJE
Marcin Borkiewicz
8 06 2020
Korelacja Najwieksza korelacja zachodzi pomiedzy PM10 i PM2.5 poniewarz maj różne srednice.
korelacja <- cor(dane)
corrplot.mixed(korelacja)
Model regresji wielokrotna dla PM10
mod1<-lm(PM10~.,data=dane,na.action = na.exclude)
ols_plot_diagnostics(mod1)
Wykres rozrzutu
p1<-scatterplotMatrix(dane, spread=FALSE, lty.smooth=2)
Kryterium informacyjne Akaikego jeżeli warto jest nizsza to model jest dokadniejszy
mod2<-lm(PM10~PM2.5+RH+SO2*NO2*NOx*O3,data=dane)
mod3<-lm(PM10~.*.,data=dane)
AIC(mod1,mod2,mod3)## df AIC
## mod1 8 36233.10
## mod2 19 35291.51
## mod3 23 35044.88Algorytm LEAPS
mAIC <- stepAIC(mod3, direction="backward")## Start: AIC=20859.17
## PM10 ~ (RH + PM2.5 + SO2 + NO2 + NOx + O3) * (RH + PM2.5 + SO2 +
## NO2 + NOx + O3)
##
## Df Sum of Sq RSS AIC
## - PM2.5:NO2 1 8.3 321747 20857
## - RH:PM2.5 1 23.5 321762 20858
## <none> 321738 20859
## - PM2.5:O3 1 281.2 322019 20862
## - SO2:O3 1 524.2 322262 20865
## - PM2.5:NOx 1 669.2 322407 20868
## - PM2.5:SO2 1 732.3 322471 20869
## - NO2:O3 1 741.8 322480 20869
## - SO2:NO2 1 1020.1 322758 20873
## - RH:SO2 1 1149.1 322887 20875
## - RH:O3 1 1377.3 323116 20879
## - RH:NO2 1 1623.0 323361 20882
## - NOx:O3 1 3761.4 325500 20915
## - RH:NOx 1 11525.5 333264 21033
## - SO2:NOx 1 11634.7 333373 21035
## - NO2:NOx 1 11810.7 333549 21037
##
## Step: AIC=20857.3
## PM10 ~ RH + PM2.5 + SO2 + NO2 + NOx + O3 + RH:PM2.5 + RH:SO2 +
## RH:NO2 + RH:NOx + RH:O3 + PM2.5:SO2 + PM2.5:NOx + PM2.5:O3 +
## SO2:NO2 + SO2:NOx + SO2:O3 + NO2:NOx + NO2:O3 + NOx:O3
##
## Df Sum of Sq RSS AIC
## - RH:PM2.5 1 17.3 321764 20856
## <none> 321747 20857
## - PM2.5:O3 1 335.2 322082 20861
## - SO2:O3 1 566.8 322313 20864
## - PM2.5:SO2 1 741.7 322488 20867
## - NO2:O3 1 743.7 322490 20867
## - PM2.5:NOx 1 997.3 322744 20871
## - RH:SO2 1 1170.9 322918 20874
## - RH:O3 1 1502.8 323249 20879
## - RH:NO2 1 1797.5 323544 20883
## - SO2:NO2 1 2024.5 323771 20887
## - NOx:O3 1 3768.4 325515 20914
## - RH:NOx 1 12242.7 333989 21042
## - NO2:NOx 1 12936.7 334683 21052
## - SO2:NOx 1 15056.4 336803 21084
##
## Step: AIC=20855.57
## PM10 ~ RH + PM2.5 + SO2 + NO2 + NOx + O3 + RH:SO2 + RH:NO2 +
## RH:NOx + RH:O3 + PM2.5:SO2 + PM2.5:NOx + PM2.5:O3 + SO2:NO2 +
## SO2:NOx + SO2:O3 + NO2:NOx + NO2:O3 + NOx:O3
##
## Df Sum of Sq RSS AIC
## <none> 321764 20856
## - PM2.5:O3 1 425.3 322189 20860
## - SO2:O3 1 550.9 322315 20862
## - NO2:O3 1 733.0 322497 20865
## - PM2.5:SO2 1 845.4 322609 20867
## - PM2.5:NOx 1 1037.7 322802 20870
## - RH:SO2 1 1276.5 323040 20873
## - RH:O3 1 1485.6 323250 20877
## - RH:NO2 1 1790.7 323555 20881
## - SO2:NO2 1 2025.3 323789 20885
## - NOx:O3 1 3804.0 325568 20912
## - RH:NOx 1 12766.0 334530 21048
## - NO2:NOx 1 12939.0 334703 21051
## - SO2:NOx 1 15242.5 337006 21085mLEAPS<-regsubsets(PM10~.,data=dane,nbest=10)
p<-plot(mLEAPS,scale="r2")
Algorytm EARTH
mEarth<-earth(PM10~.,dane)
ev<-evimp (mEarth)
plot(ev)
plot(mEarth)
Algorytm BORUTA
mBoruta<-Boruta(PM10~.,dane)
boruta_signif <- names(mBoruta$finalDecision[mBoruta$finalDecision %in% c("Confirmed", "Tentative")])
plot(mBoruta, cex.axis=.8, las=3, xlab="")
Walicaja krzywoliniowa
Control <- trainControl(method = "cv", number = 10, verboseIter = F)
mCV <- caret::train(PM10~.,
data=dane,
method = "lm",
metric = 'MAE',
trControl = Control,
na.action = na.omit,
tuneLength = 15
)
wCV<-predict(mCV)
plot(dane$PM10,wCV,ylab="wartosci predykcji", xlab="wartosc PM10",main="")
gof(wCV,dane$PM10)## [,1]
## ME 0.00
## MAE 5.49
## MSE 82.14
## RMSE 9.06
## NRMSE % 35.90
## PBIAS % 0.00
## RSR 0.36
## rSD 0.93
## NSE 0.87
## mNSE 0.68
## rNSE 0.88
## d 0.96
## md 0.83
## rd 0.97
## cp 0.16
## r 0.93
## R2 0.87
## bR2 0.84
## KGE 0.91
## VE 0.85ggof(wCV,dane$PM10,ylab="wartosci", xlab="nr obserwacji",main="")
plot(dane$PM10-wCV,ylab="wartosci", xlab="nr obserwacji",main="RÄ‚Å‚znica PM10 i wartosci predykcji")
plot(abs(dane$PM10-wCV),ylab="wartosci", xlab="nr obserwacji",main="Wartoci bezwgledne z różnica PM10 i wartosci predykcji")
plot(sort(abs(dane$PM10-wCV)),ylab="wartosci", xlab="nr obserwacji",main="Posortowane wartosci bezwzgledne")
Regresja
reg1<-lm(dane$PM10~dane$PM2.5)
ggplot(dane,aes(x = PM10, y = PM2.5)) +
geom_point() +
stat_smooth(method = "lm", col = "red")
Przedzial ufnosci
confint(reg1)## 2.5 % 97.5 %
## (Intercept) 12.6559248 13.5334823
## dane$PM2.5 0.9736836 0.9999065Badanie rozkadu PM2.5
RH<-dane[,1]
PM10<-dane[,2]
PM2.5<-dane[,3]
SO2<-dane[,4]
NO2<-dane[,5]
NOx<-dane[,6]
O3<-dane[,7]
plot(PM2.5)
shapiro.test(PM2.5[1:2000])##
## Shapiro-Wilk normality test
##
## data: PM2.5[1:2000]
## W = 0.82021, p-value < 2.2e-16Badanie rozkadu O3
plot(O3)
shapiro.test(O3[1:2000])##
## Shapiro-Wilk normality test
##
## data: O3[1:2000]
## W = 0.96429, p-value < 2.2e-16Wykres reszt
regresja<-lm(PM2.5~PM10+RH+SO2+NO2+NOx+O3)
reszty<-resid(regresja)
reszty_sort<-sort(reszty)
boxplot(reszty)
plot(reszty_sort)
Wykres modelu predykcji PM2.5
model<-fitted(regresja)
model_pre<-predict(regresja)
plot(model_pre)
Test modelu predykcji
testowy=PM2.5~model
plot(testowy)
Inflacja wariacji VIF
vif(regresja)## PM10 RH SO2 NO2 NOx O3
## 2.769157 2.464370 1.821476 4.294749 3.588674 4.132260Wykres modelu regresji PM2.5
plot(regresja)
Walidacja zalozen modelu regresji PM2.5
gvlma(regresja)##
## Call:
## lm(formula = PM2.5 ~ PM10 + RH + SO2 + NO2 + NOx + O3)
##
## Coefficients:
## (Intercept) PM10 RH SO2 NO2 NOx
## -20.80688 0.71343 0.26127 1.17993 -0.12596 -0.02046
## O3
## -0.02231
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = regresja)
##
## Value p-value Decision
## Global Stat 104529.1 0 Assumptions NOT satisfied!
## Skewness 3994.2 0 Assumptions NOT satisfied!
## Kurtosis 99921.2 0 Assumptions NOT satisfied!
## Link Function 168.6 0 Assumptions NOT satisfied!
## Heteroscedasticity 445.1 0 Assumptions NOT satisfied!Model regresji wielorakiej dla PM2.5 metoda najmniejszych kwadratów
ols_plot_diagnostics(regresja)
