cat("
<style>
h1 {
  font-size: 5rem;
}

h2 {
  font-size: 4rem;
}

h3 {
  font-size: 3rem;
}

h4 {
  font-size: 2rem;
}

.tocify-header {
  font-size: 2rem;
}

.tocify-header {
  font-size: 2rem;
}

.tocify-subheader {
  font-size: 1.5rem;
}
</style>
")

Przygotowanie danych

library(tidyverse)
library(worldmet)
library(car)
library(MASS)
library(leaps)
library(earth)
library(Boruta)
library(caret)
library(hydroGOF)
library(ggplot2)

library(relaimpo)
library(ggfortify)
library(DT)

Dane pobrano z bazy danych NOAA. Zawierają one dane meteorologiczne dla stacji znajdującej się na Śnieżce. Regresja będzie badana dla danych z roku 2020. Dane podzielono na grupę walidacyjną oraz testową.

#dane_Sniezka <- importNOAA(code = "125100-99999", year = c(2019,2020), hourly = T, quiet = FALSE)

dane_Sniezka <- readRDS("dane_Sniezka.rds")

dane_Sniezka_wybrane_kolumny <- dane_Sniezka[,c(7,8,9,11,12,13,15,24)]
dane_Sniezka_wybrane_kolumny$cl_1[is.na(dane_Sniezka_wybrane_kolumny$cl_1)] <- 0
dane_gotowe <- dane_Sniezka_wybrane_kolumny[complete.cases(dane_Sniezka_wybrane_kolumny),]
colnames(dane_gotowe) <- c("Predkosc_wiatru","Kierunek_wiatru","Temperatura","Widocznosc","Temperatura_punktu_rosy","Wilgotnosc","Pokrycie_chmur","Opady_1h")

dane_do_modelowania <- dane_gotowe[1:4287,]
dane_do_testowania_modelu <- dane_gotowe[4287:8574,]

Modele

Model dla Temperatury w zależności od reszty wartości

m_all=lm(Temperatura~.,data=dane_do_modelowania, na.action = na.exclude)
summary(m_all)

## 
## Call:
## lm(formula = Temperatura ~ ., data = dane_do_modelowania, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6637 -0.7127  0.0089  0.4685 16.1974 
## 
## Coefficients:
##                           Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)              2.640e+01  2.065e-01  127.847  < 2e-16 ***
## Predkosc_wiatru         -9.510e-03  3.915e-03   -2.429   0.0152 *  
## Kierunek_wiatru          6.200e-04  2.543e-04    2.439   0.0148 *  
## Widocznosc              -4.059e-05  1.503e-06  -27.010  < 2e-16 ***
## Temperatura_punktu_rosy  9.558e-01  4.208e-03  227.134  < 2e-16 ***
## Wilgotnosc              -2.690e-01  1.929e-03 -139.490  < 2e-16 ***
## Pokrycie_chmur           5.428e-02  1.005e-02    5.402 6.96e-08 ***
## Opady_1h                 4.381e-02  9.790e-03    4.474 7.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.473 on 4279 degrees of freedom
## Multiple R-squared:  0.9367, Adjusted R-squared:  0.9366 
## F-statistic:  9046 on 7 and 4279 DF,  p-value: < 2.2e-16

Model dla Temperatury w zależności od reszty wartości wraz z ich interakcjami

m_all_with_all=lm(Temperatura~.*.,data=dane_do_modelowania)
summary(m_all_with_all)

## 
## Call:
## lm(formula = Temperatura ~ . * ., data = dane_do_modelowania)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5672 -0.1163  0.0274  0.1310  4.6851 
## 
## Coefficients:
##                                           Estimate Std. Error t value Pr(>|t|)
## (Intercept)                              1.277e+01  4.379e-01  29.170  < 2e-16
## Predkosc_wiatru                         -1.842e-01  1.697e-02 -10.854  < 2e-16
## Kierunek_wiatru                          1.072e-02  9.453e-04  11.341  < 2e-16
## Widocznosc                               6.737e-05  7.622e-06   8.839  < 2e-16
## Temperatura_punktu_rosy                  2.928e-01  8.727e-03  33.551  < 2e-16
## Wilgotnosc                              -1.248e-01  4.465e-03 -27.943  < 2e-16
## Pokrycie_chmur                          -2.330e-01  4.685e-02  -4.974 6.80e-07
## Opady_1h                                -5.387e-02  1.212e-01  -0.444  0.65672
## Predkosc_wiatru:Kierunek_wiatru          1.086e-04  2.427e-05   4.474 7.88e-06
## Predkosc_wiatru:Widocznosc               5.296e-07  1.153e-07   4.594 4.47e-06
## Predkosc_wiatru:Temperatura_punktu_rosy -1.531e-04  3.253e-04  -0.471  0.63795
## Predkosc_wiatru:Wilgotnosc               1.623e-03  1.503e-04  10.795  < 2e-16
## Predkosc_wiatru:Pokrycie_chmur          -2.083e-03  7.207e-04  -2.890  0.00387
## Predkosc_wiatru:Opady_1h                -9.915e-04  7.416e-04  -1.337  0.18132
## Kierunek_wiatru:Widocznosc              -3.162e-08  6.677e-09  -4.736 2.26e-06
## Kierunek_wiatru:Temperatura_punktu_rosy  2.971e-05  1.789e-05   1.661  0.09676
## Kierunek_wiatru:Wilgotnosc              -1.228e-04  9.630e-06 -12.755  < 2e-16
## Kierunek_wiatru:Pokrycie_chmur           8.714e-05  4.881e-05   1.785  0.07427
## Kierunek_wiatru:Opady_1h                -6.525e-06  6.410e-05  -0.102  0.91893
## Widocznosc:Temperatura_punktu_rosy       2.950e-06  9.251e-08  31.888  < 2e-16
## Widocznosc:Wilgotnosc                   -7.339e-07  7.628e-08  -9.621  < 2e-16
## Widocznosc:Pokrycie_chmur                1.244e-06  2.760e-07   4.505 6.80e-06
## Widocznosc:Opady_1h                     -3.974e-07  5.347e-07  -0.743  0.45739
## Temperatura_punktu_rosy:Wilgotnosc       7.320e-03  6.466e-05 113.210  < 2e-16
## Temperatura_punktu_rosy:Pokrycie_chmur  -5.394e-03  7.581e-04  -7.116 1.30e-12
## Temperatura_punktu_rosy:Opady_1h        -2.033e-03  9.114e-04  -2.231  0.02576
## Wilgotnosc:Pokrycie_chmur                2.340e-03  4.532e-04   5.163 2.54e-07
## Wilgotnosc:Opady_1h                      6.592e-04  1.214e-03   0.543  0.58711
## Pokrycie_chmur:Opady_1h                  2.066e-03  1.928e-03   1.072  0.28398
##                                            
## (Intercept)                             ***
## Predkosc_wiatru                         ***
## Kierunek_wiatru                         ***
## Widocznosc                              ***
## Temperatura_punktu_rosy                 ***
## Wilgotnosc                              ***
## Pokrycie_chmur                          ***
## Opady_1h                                   
## Predkosc_wiatru:Kierunek_wiatru         ***
## Predkosc_wiatru:Widocznosc              ***
## Predkosc_wiatru:Temperatura_punktu_rosy    
## Predkosc_wiatru:Wilgotnosc              ***
## Predkosc_wiatru:Pokrycie_chmur          ** 
## Predkosc_wiatru:Opady_1h                   
## Kierunek_wiatru:Widocznosc              ***
## Kierunek_wiatru:Temperatura_punktu_rosy .  
## Kierunek_wiatru:Wilgotnosc              ***
## Kierunek_wiatru:Pokrycie_chmur          .  
## Kierunek_wiatru:Opady_1h                   
## Widocznosc:Temperatura_punktu_rosy      ***
## Widocznosc:Wilgotnosc                   ***
## Widocznosc:Pokrycie_chmur               ***
## Widocznosc:Opady_1h                        
## Temperatura_punktu_rosy:Wilgotnosc      ***
## Temperatura_punktu_rosy:Pokrycie_chmur  ***
## Temperatura_punktu_rosy:Opady_1h        *  
## Wilgotnosc:Pokrycie_chmur               ***
## Wilgotnosc:Opady_1h                        
## Pokrycie_chmur:Opady_1h                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6095 on 4258 degrees of freedom
## Multiple R-squared:  0.9892, Adjusted R-squared:  0.9891 
## F-statistic: 1.395e+04 on 28 and 4258 DF,  p-value: < 2.2e-16

Wykresy rozrzutów dla wszystkich danych

scatterplotMatrix(dane_do_modelowania, spread=FALSE, lty.smooth=2)

Regresja krokowa wsteczna

Można tu wyczytać, że największą wartość dla modelu ma temperatura punktu rosy, a najmniejszą prędkość wiatru (przy pozostawienu tylko tych wartości w modelu).

m_AIC=stepAIC(m_all,direction = "backward")

## Start:  AIC=3329.09
## Temperatura ~ Predkosc_wiatru + Kierunek_wiatru + Widocznosc + 
##     Temperatura_punktu_rosy + Wilgotnosc + Pokrycie_chmur + Opady_1h
## 
##                           Df Sum of Sq    RSS     AIC
## <none>                                   9285  3329.1
## - Predkosc_wiatru          1        13   9298  3333.0
## - Kierunek_wiatru          1        13   9298  3333.0
## - Opady_1h                 1        43   9329  3347.1
## - Pokrycie_chmur           1        63   9348  3356.2
## - Widocznosc               1      1583  10868  4002.0
## - Wilgotnosc               1     42221  51506 10672.0
## - Temperatura_punktu_rosy  1    111946 121231 14341.6

summary(m_AIC)

## 
## Call:
## lm(formula = Temperatura ~ Predkosc_wiatru + Kierunek_wiatru + 
##     Widocznosc + Temperatura_punktu_rosy + Wilgotnosc + Pokrycie_chmur + 
##     Opady_1h, data = dane_do_modelowania, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6637 -0.7127  0.0089  0.4685 16.1974 
## 
## Coefficients:
##                           Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)              2.640e+01  2.065e-01  127.847  < 2e-16 ***
## Predkosc_wiatru         -9.510e-03  3.915e-03   -2.429   0.0152 *  
## Kierunek_wiatru          6.200e-04  2.543e-04    2.439   0.0148 *  
## Widocznosc              -4.059e-05  1.503e-06  -27.010  < 2e-16 ***
## Temperatura_punktu_rosy  9.558e-01  4.208e-03  227.134  < 2e-16 ***
## Wilgotnosc              -2.690e-01  1.929e-03 -139.490  < 2e-16 ***
## Pokrycie_chmur           5.428e-02  1.005e-02    5.402 6.96e-08 ***
## Opady_1h                 4.381e-02  9.790e-03    4.474 7.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.473 on 4279 degrees of freedom
## Multiple R-squared:  0.9367, Adjusted R-squared:  0.9366 
## F-statistic:  9046 on 7 and 4279 DF,  p-value: < 2.2e-16

Najlepsza konfiguracja występuje przy uwzględnienu wszystkich parametrów.

m_AIC$anova

Regresja krokowa LEAPS

Z poniższego wykresu wynika, że najlepszym parametrem jest temperatura punktu rosy, ponieważ występuje ona w największej liczbie modeli o największym R². Wysoki wpływ na model ma też wilgotność i widoczność. Aby uzyskać poziom R²= 94%, można wykorzystać do modelu wcześniej podane parametry.

mLEAPS=regsubsets(Temperatura~., data=dane_do_modelowania, nbest=5)
#mLEAPS
#summary(mLEAPS)
plot(mLEAPS,scale="r2")

m_Leaps=lm(Temperatura~Temperatura_punktu_rosy+Widocznosc+Wilgotnosc, data=dane_do_modelowania, na.action = na.exclude)

Algorytm Earth

Poniższy wykres wskazuje, że najważniejszym parametrem jest temperatura punktu rosy, a kolejnym wilgotność. Na podstawie tych parametrów można utworzyć model.

mEarth=earth(Temperatura~.,dane_do_modelowania)
#plot(mEarth)
ev=evimp(mEarth)
plot(ev)

m_Earth=lm(Temperatura~Temperatura_punktu_rosy+Wilgotnosc, data=dane_do_modelowania, na.action = na.exclude)

Algorytm Boruta

Z wykresu stworzonego na bazie algorytmu Boruta wynika, że najważniejszymi parametrami kolejno są temperatura punktu rosy, wilgotność i kierunek wiatru. Pokrycie chmur i opady mają najmniejszy wpływ na model. Na podstawie najważniejszych parametrów można utworzyć model.

mBoruta=Boruta(Temperatura~.,dane_do_modelowania)
boruta_signif <- names(mBoruta$finalDecision[mBoruta$finalDecision %in% c("Confirmed", "Tentative")])  
#print(boruta_signif) 
plot(mBoruta, cex.axis=.9, las=2, xlab="")

m_Boruta=lm(Temperatura~Temperatura_punktu_rosy+Wilgotnosc+Kierunek_wiatru, data=dane_do_modelowania, na.action = na.exclude)

Algorytmy pakietu relaimpo

Algorytmy używają bootstrapowe przedziały ufności dla oceny relatywnej wagi parametru. Boot.relimp używa pakietu R, aby wykonać rzeczywisty bootstrap wybranych metryk, podczas gdy booteval.relimp ocenia wyniki i podaje przedziały ufności.

Metryka LMG

W metryce lmg największą wartość dla modelu wnoszą kolejno temperatura punktu rosy, wilgotność oraz widoczność. Najmniejszą ma kierunek wiatru. Na podstawie najważniejszych parametrów można utworzyć model.

plot(booteval.relimp(boot,sort=TRUE))

Metryka first

W metryce first największą wartość dla modelu wnoszą kolejno temperatura punktu rosy, prędkość wiatru oraz widoczność. Najmniejszą ma opad. Na podstawie najważniejszych parametrów można utworzyć model.

plot(booteval.relimp(boot,sort=TRUE))

Metryka pratt

W metryce pratt największą wartość dla modelu wnoszą kolejno temperatura punktu rosy, wilgotność oraz widoczność. Najmniejszą ma kierunek wiatru. Na podstawie najważniejszych parametrów można utworzyć model.

plot(booteval.relimp(boot,sort=TRUE))

Walidacja krzyżowa w caret

Control <- trainControl(method = "cv", number = 100, verboseIter = F)
mCV <- caret::train(Temperatura~.,
                    data=dane_do_modelowania,
                    method = "lm",
                    metric = 'MAE', 
                    trControl = Control,
                    na.action = na.omit,
                    tuneLength = 15
)

mCV

## Linear Regression 
## 
## 4287 samples
##    7 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (100 fold) 
## Summary of sample sizes: 4243, 4243, 4244, 4244, 4244, 4243, ... 
## Resampling results:
## 
##   RMSE     Rsquared   MAE      
##   1.38418  0.9383278  0.8776339
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

summary(mCV)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6637 -0.7127  0.0089  0.4685 16.1974 
## 
## Coefficients:
##                           Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)              2.640e+01  2.065e-01  127.847  < 2e-16 ***
## Predkosc_wiatru         -9.510e-03  3.915e-03   -2.429   0.0152 *  
## Kierunek_wiatru          6.200e-04  2.543e-04    2.439   0.0148 *  
## Widocznosc              -4.059e-05  1.503e-06  -27.010  < 2e-16 ***
## Temperatura_punktu_rosy  9.558e-01  4.208e-03  227.134  < 2e-16 ***
## Wilgotnosc              -2.690e-01  1.929e-03 -139.490  < 2e-16 ***
## Pokrycie_chmur           5.428e-02  1.005e-02    5.402 6.96e-08 ***
## Opady_1h                 4.381e-02  9.790e-03    4.474 7.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.473 on 4279 degrees of freedom
## Multiple R-squared:  0.9367, Adjusted R-squared:  0.9366 
## F-statistic:  9046 on 7 and 4279 DF,  p-value: < 2.2e-16

Diagnostyka modelu

Kryterium informacyjne Akaike

Na podstawie kryterium AIC dwa modele wyróżniają się na tle innych. Najbardziej dokładny jest model oparty na wszystkich parametrach wraz z ich interakcjami, natomiast najmniej dokładny model z pakietu relaimpo z metryką first. Warto zwrócić też uwagę na modele, które wybrały tylko część parametrów, a porawdziły sobie równie dobrze jak model ze wszystkimi parametrami. Są to modele Leaps, ralaimpo z metryką lmg / pratt, Boruta oraz Earth. Poniżej zostaną zaprezentoowane wykresy i dokładniejsze dane na temat tych modeli.

datatable(AIC(m_all,m_all_with_all,m_AIC,m_Leaps,m_Earth,m_Boruta,m_relaimpo_first,m_relaimpo_lmg,m_relaimpo_pratt))

Wszystkie dane

aa=predict(m_all)
plot(dane_do_modelowania$Temperatura,aa)

ggof(aa,dane_do_modelowania$Temperatura)

Wszystkie dane z interakcjami

bb=predict(m_all_with_all)
plot(dane_do_modelowania$Temperatura,bb)

ggof(bb,dane_do_modelowania$Temperatura)

mAIC

cc=predict(m_AIC)
plot(dane_do_modelowania$Temperatura,cc)

ggof(cc,dane_do_modelowania$Temperatura)

LEAPS

dd=predict(m_Leaps)
plot(dane_do_modelowania$Temperatura,dd)

ggof(dd,dane_do_modelowania$Temperatura)

Earth

ee=predict(m_Earth)
plot(dane_do_modelowania$Temperatura,ee)

ggof(ee,dane_do_modelowania$Temperatura)

Boruta

ff=predict(m_Boruta)
plot(dane_do_modelowania$Temperatura,ff)

ggof(ff,dane_do_modelowania$Temperatura)

Relaimpo

LMG

gg=predict(m_relaimpo_lmg)
plot(dane_do_modelowania$Temperatura,gg)

ggof(gg,dane_do_modelowania$Temperatura)

first

hh=predict(m_relaimpo_first)
plot(dane_do_modelowania$Temperatura,hh)

ggof(hh,dane_do_modelowania$Temperatura)

pratt

ii=predict(m_relaimpo_pratt)
plot(dane_do_modelowania$Temperatura,ii)

ggof(ii,dane_do_modelowania$Temperatura)

mCV

jj=predict(mCV)
plot(dane_do_modelowania$Temperatura,jj)

ggof(jj,dane_do_modelowania$Temperatura)

Testowanie i diagnostyka wybranych modeli

Do testowania wybrano modele utworzone na podstawie kilku wybranych parametrów:

1)temperatura punktu rosy, wilgotność, widoczność (Leaps, relaimpo - lmg / pratt,mCV)

2)temperatura punktu rosy, wilgotność, kierunek wiatru (Boruta)

3)temperatura punktu rosy, wilgotność (Earth)

Diagnostyka modeli z pakietu ggfortify

Residuals vs fitted sprawdza zależność liniową. Pozioma linia wskazuje na dobrą zależność. Normal Q-Q sprawdza, czy reszty mają rozkład normalny. Punkty pokrywające się z linią przerywaną wskazują na rozkład normalny. Scale-Location służy do sprawdzania jednorodności wariancji reszt (homoskedastyczności). Pozioma linia z równo rozłożonymi punktami jest dobrym wskaźnikiem homoskedastyczności. Residuals vs Leverage służy do identyfikacji wartości ekstremalnych mogących wpływać na wyniki regresji.

Leaps, relaimpo - lmg / pratt, mCV

autoplot(m_Leaps)

Boruta

autoplot(m_Boruta)

Earth

autoplot(m_Earth)

### Testowanie

Leaps, relaimpo - lmg / pratt, mCV

mt1=lm(Temperatura~Temperatura_punktu_rosy+Wilgotnosc+Widocznosc,data=dane_do_testowania_modelu, na.action = na.exclude)
t1=predict(mt1)

dane_t1=data.frame(dane_do_testowania_modelu$Temperatura,t1)

ggplot(dane_t1,aes(x=dane_do_testowania_modelu.Temperatura,y=t1))+
  geom_point()+
  stat_smooth(method="lm",col="red",se=TRUE)

ggof(t1,dane_do_testowania_modelu$Temperatura)

Boruta

mt2=lm(Temperatura~Temperatura_punktu_rosy+Wilgotnosc+Kierunek_wiatru,data=dane_do_testowania_modelu, na.action = na.exclude)
t2=predict(mt2)

dane_t2=data.frame(dane_do_testowania_modelu$Temperatura,t2)

ggplot(dane_t2,aes(x=dane_do_testowania_modelu.Temperatura,y=t2))+
  geom_point()+
  stat_smooth(method="lm",col="red",se=TRUE)

ggof(t2,dane_do_testowania_modelu$Temperatura)

Earth

mt3=lm(Temperatura~Temperatura_punktu_rosy+Wilgotnosc,data=dane_do_testowania_modelu, na.action = na.exclude)
t3=predict(mt3)

dane_t3=data.frame(dane_do_testowania_modelu$Temperatura,t3)

ggplot(dane_t3,aes(x=dane_do_testowania_modelu.Temperatura,y=t3))+
  geom_point()+
  stat_smooth(method="lm",col="red",se=TRUE)

ggof(t3,dane_do_testowania_modelu$Temperatura)

Wnioski

Najskuteczniejszym modelem bez względu na liczbę parametrów jest model ze wszystkimi danymi wraz z ich interakcjami. Bardzo dobrze działa również model z małą liczbą parametrów - temperaturą punktu rosy, wilgotnością oraz widocznością, które są również najistotniejszymi parametrami według większości metod. Temperatura punktu rosy jednak nie powinna być stosowana do modelu z tego względu, że wyznacza się ją właśnie z temperatury, natomiast wilgotność i widoczność może być w pewien sposób skorelowana z temperaturą, dzięki czemu ma taki wysoki wpływ przy analizie regresji.

Regresja

Dominik Kapusta

15 04 2021

Przygotowanie danych

Modele

Model dla Temperatury w zależności od reszty wartości

Model dla Temperatury w zależności od reszty wartości wraz z ich interakcjami

Wykresy rozrzutów dla wszystkich danych

Regresja krokowa wsteczna

Regresja krokowa LEAPS

Algorytm Earth

Algorytm Boruta

Algorytmy pakietu relaimpo

Metryka LMG

Metryka first

Metryka pratt

Walidacja krzyżowa w caret

Diagnostyka modelu

Kryterium informacyjne Akaike

Wszystkie dane

Wszystkie dane z interakcjami

mAIC

LEAPS

Earth

Boruta

Relaimpo

LMG

first

pratt

mCV

Testowanie i diagnostyka wybranych modeli

Diagnostyka modeli z pakietu ggfortify

Leaps, relaimpo - lmg / pratt, mCV

Boruta

Earth

Leaps, relaimpo - lmg / pratt, mCV

Boruta

Earth

Wnioski