Pozyskanie danych

Do przeprowadzenia analizy wybrano dane meteorologiczne. Do ich uzyskania zastosowano funkcję importNOAA z pakietu worldmet. Do wykonania analizy wybrano dane z jednej ze stacji w Krakowie - Balice, z lat 2004-2019. Poniżej widoczne jest wczytanie danych oraz zapisanie ich do pliku.

setwd("C:/PDS/cw3")
getMeta()
metaPL <- getMeta(country = "PL")

importNOAA(code = "125660-99999", year = 2004:2019) %>%
  as.tbl() -> meteo_kra
meteo_kra %>%
  select(date, ws, wd, air_temp, atmos_pres, visibility, dew_point, RH) -> meteo_kra
save(meteo_kra, file = "meteo_kra.Rdata")

Uśrednianie danych i sprawdzenie kompletności

Za pomocą funkcji timeAverage utworzono zmienne zawierające informacje o danych średnio-miesięcznych oraz średnio-dobowych. Nastepnie utworzono wykres szeregu czasowego dla danych dobowych w celu sprawdzenia kompletnosci danych. Jak widać, żadna substancja nie posiada dużych braków danych, więc możliwe było uwzględnienie ich w dalszej analizie.

load("meteo_kra.Rdata")
dane_mies <- meteo_kra %>% timeAverage(avg.time = "month")
dane_day <- meteo_kra %>% timeAverage(avg.time = "day")
timePlot(mydata = dane_day, pollutant = c("ws", "wd", "air_temp", "atmos_pres", "visibility", "dew_point", "RH"), y.relation = "free")

Do przeprowadzenia modelowania oraz diagnostyki wykorzystano dane miesięczne, które podzielono na dwa zestawy - treningowe i testowe (2019 rok). Jako zmienną prognozowaną wybrano wilgotność wzgledną (RH). Natomiast jej predykotry to pozostałe parametry:

ws - prędkość wiatru
wd - kierunek wiatru
air_temp - temperatura powietrza
atmos_pres - ciśnienie atmosferyczne
visibility - widoczność
dew_point - temperatura punktu rosy

treningowe <- dane_mies %>% filter(year(date) < 2019)
testowe <- dane_mies %>% filter(year(date) == 2019)

Za pomocą funkcji as_tsibble zamieniono dane na obiekty tsibble, które przechowują szeregi czasowe.

tsdane <- dane_mies %>% 
  mutate(date = yearmonth(date)) %>% 
  as_tsibble(index = date)
tstren <- treningowe %>% 
  mutate(date = yearmonth(date)) %>% 
  as_tsibble(index = date)
tstest <- testowe %>% 
  mutate(date = yearmonth(date)) %>% 
  as_tsibble(index = date)

Analiza danych

Wykresy korelacji

Przed budową modeli zbadano korelacje wszystkich zmiennych. Między zmiennymi air_temp i dew_point występuje zbyt duża korelacja (0,98), wiec postanowiono wykluczyć temperaturę punktu rosy z dalszej analizy. Predyktory najbardziej skorelowane ze zmienną prognozowaną to: temperatura powietrza, widoczność oraz ciśnienie atmosferyczne.

treningowe %>%
  GGally::ggpairs(columns = 2:7)

M=cor(treningowe[,-1])
corrplot(M, method = "circle")

Wykresy rozrzutu

Poniżej widoczne są wykresy rozrzutu dla poszczególnych zmiennych.

scatterplotMatrix(treningowe[,-1], spread=FALSE, lty.smooth=2, col = carPalette())

W celu wstępnego zbadania korelacji między prognozowaną zmienną a przykladowymi predyktorami utworzono wykresy rozrzutu. Na poniższym wykresie widać dość wysoką korelację ujemną dla RH i temperatury powietrza (-0.68) oraz widoczności (-0.64).

ggplot(data = treningowe, aes(x = RH, y = air_temp)) +
  geom_point(color = "firebrick4") +
  geom_smooth(color = "dark blue") +
  labs(title = "Zależność między wilgotnoscią a temperatura") +
  xlab("wilgotnosc wzgledna") +
  ylab("temperatura powietrza")

cor(tstren[,"air_temp"], tstren[,"RH"], method = "pearson")

##                  RH
## air_temp -0.6778432

ggplot(data = treningowe, aes(x = RH, y = visibility)) +
  geom_point(color = "springgreen4") +
  geom_smooth(color = "purple") +
  labs(title = "Zależność między wilgotnoscią a widocznoscia") +
  xlab("wilgotnosc wzgledna") +
  ylab("widocznosc")

cor(tstren[,"visibility"], tstren[,"RH"], method = "pearson")

##                    RH
## visibility -0.6375119

Wykres czasowy

Na poniższym wykresie widoczna jest sezonowość zmiennej prognozowanej.

autoplot(tstren, RH) +
  labs(title = "Szereg czasowy dla zmiennej RH",
       x = "Czas", 
       y = "RH")

Regresja liniowa prosta

Zmienna prognozowana i temperatura powietrza

Przetestowano model regresji liniowej dla zmiennej prognozowanej i temperatury powietrza. Współczynnik determinacji jest stosunkowo niski (0.46).

r1=lm(RH~air_temp, treningowe)
summary(r1)

## 
## Call:
## lm(formula = RH ~ air_temp, data = treningowe)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.856  -4.124   1.036   3.918  12.281 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 85.52193    0.64311   133.0   <2e-16 ***
## air_temp    -0.66595    0.05414   -12.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.634 on 178 degrees of freedom
## Multiple R-squared:  0.4595, Adjusted R-squared:  0.4564 
## F-statistic: 151.3 on 1 and 178 DF,  p-value: < 2.2e-16

hist(r1$residuals, main = "Histogram reszt",
     xlab = "Reszty", ylab = "Czestotliowsc")

Wykres zaleznosci modelu od zmiennej prognozowanej

Poniżej widoczny jest wykres zależności modelu od wilgotności względnej. Można na nim zauważyć dobre wpasowanie danych w model.

m1 = fitted(r1)
RH = treningowe$RH
plot(RH,m1, main = "Wykres zaleznosci modelu od RH",
     ylab = "Model")

cor(treningowe[,"RH"], m1, method = "pearson")

##         [,1]
## RH 0.6778432

Analiza reszt

W dobrym modelu suma reszt jest równa 0, dlatego postanowiono sprawdzić to założenie dla wykonanego modelu. Jak widać poniżej suma oraz średnia są bardzo zbliżone do wartości 0.

res1 = r1$residuals
sum(res1)

## [1] -6.078471e-14

mean(res1)

## [1] -3.379434e-16

Za pomocą testu Shapiro-Wilka sprawdzono, czy reszty mają rozkład normalny. W tym przypadku wartość p-value jest mniejsza od 0.05, wiec warunek nie został spełniony. Hipotezę alternatywną potwierdził również test Jarque-Bera.

shapiro.test(res1)

## 
##  Shapiro-Wilk normality test
## 
## data:  res1
## W = 0.97252, p-value = 0.001265

jarque.test(res1)

## 
##  Jarque-Bera Normality Test
## 
## data:  res1
## JB = 9.2148, p-value = 0.009978
## alternative hypothesis: greater

hist(res1, main = "Histogram reszt dla pierwszego modelu",
     xlab = "Reszty", ylab = "Czestotliwosc")

Za pomocą testu Durbin-Watson zbadano autokorelację reszt. W przypadku badanego modelu odrzucono hipotezę zerową, która twierdzi, że korelacja między resztami nie wystepuje i przyjęto hipotezę alternatywną. Potwierdza to poniższy wykres.

dwtest(r1)

## 
##  Durbin-Watson test
## 
## data:  r1
## DW = 0.84941, p-value = 2.729e-15
## alternative hypothesis: true autocorrelation is greater than 0

acf(res1,lag.max = 100)

Za pomocą funkcji gvlma zbadano skośność, kurtozę oraz heteroskedastyczność wariancji (wewnątrz wariancji wszystkie błędy są równe). Jak widać poniżej, tylko kurtoza, heteroskedastyczność oraz link function spełniają założenia.

gvlma(r1)

## 
## Call:
## lm(formula = RH ~ air_temp, data = treningowe)
## 
## Coefficients:
## (Intercept)     air_temp  
##      85.522       -0.666  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = r1) 
## 
##                        Value  p-value                   Decision
## Global Stat        1.029e+01 0.035851 Assumptions NOT satisfied!
## Skewness           9.214e+00 0.002402 Assumptions NOT satisfied!
## Kurtosis           7.479e-04 0.978182    Assumptions acceptable.
## Link Function      2.753e-01 0.599773    Assumptions acceptable.
## Heteroscedasticity 7.975e-01 0.371827    Assumptions acceptable.

Wykresy diagnostyczne dla modelu

ols_plot_diagnostics(r1)

Wpasowanie danych w model

Jak widać na poniższym wykresie model dość dobrze wpasował dane. Ekstrema znajdują się w tych samych przedziałach czasowych, ale model znacząco ich nie doszacował.

ggof(m1, treningowe$RH)

Zmienna prognozowana i widoczność

Przetestowano model regresji liniowej dla zmiennej prognozowanej i widoczności. Współczynnik determinacji jest stosunkowo niski (0.41).

r2=lm(RH~visibility, treningowe)
summary(r2)

## 
## Call:
## lm(formula = RH ~ visibility, data = treningowe)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.174  -3.777   1.131   4.531   9.926 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 89.2789994  0.9866019   90.49   <2e-16 ***
## visibility  -0.0007000  0.0000634  -11.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.905 on 178 degrees of freedom
## Multiple R-squared:  0.4064, Adjusted R-squared:  0.4031 
## F-statistic: 121.9 on 1 and 178 DF,  p-value: < 2.2e-16

hist(r2$residuals, main = "Histogram reszt",
     xlab = "Reszty", ylab = "Czestotliowsc")

Wykres zaleznosci modelu od zmiennej prognozowanej

Poniżej widoczny jest wykres zależności modelu od wilgotności względnej. Można na nim zauważyć niezbyt dobre wpasowanie danych w model.

m2 = fitted(r2)
vis = treningowe$visibility
plot(RH,m2, main = "Wykres zaleznosci modelu od RH",
     ylab = "Model")

Analiza reszt

Jak widać poniżej suma oraz średnia są bardzo zbliżone do wartości 0.

res2 = r2$residuals
sum(res2)

## [1] -1.073724e-13

mean(res2)

## [1] -5.962919e-16

shapiro.test(res2)

## 
##  Shapiro-Wilk normality test
## 
## data:  res2
## W = 0.94905, p-value = 4.661e-06

jarque.test(res2)

## 
##  Jarque-Bera Normality Test
## 
## data:  res2
## JB = 16.836, p-value = 0.0002209
## alternative hypothesis: greater

hist(res2, main = "Histogram reszt dla pierwszego modelu",
     xlab = "Reszty", ylab = "Czestotliwosc")

dwtest(r2)

## 
##  Durbin-Watson test
## 
## data:  r2
## DW = 0.72809, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

acf(res2,lag.max = 100)

Za pomocą funkcji gvlma zbadano skośność, kurtozę oraz heteroskedastyczność wariancji (wewnątrz wariancji wszystkie błędy są równe). Jak widać poniżej, tylko kurtoza spełnia wymagane założenia.

gvlma(r2)

## 
## Call:
## lm(formula = RH ~ visibility, data = treningowe)
## 
## Coefficients:
## (Intercept)   visibility  
##     89.2790      -0.0007  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = r2) 
## 
##                       Value   p-value                   Decision
## Global Stat        47.06223 1.480e-09 Assumptions NOT satisfied!
## Skewness           16.77699 4.204e-05 Assumptions NOT satisfied!
## Kurtosis            0.05862 8.087e-01    Assumptions acceptable.
## Link Function      24.57785 7.137e-07 Assumptions NOT satisfied!
## Heteroscedasticity  5.64878 1.747e-02 Assumptions NOT satisfied!

Wykresy diagnostyczne dla modelu

ols_plot_diagnostics(r2)

Wpasowanie danych w model

Jak widać na poniższym wykresie model słabo wpasował dane. Ekstrema znajdują się w tych samych przedziałach czasowych, ale model znacząco ich nie doszacował.

ggof(m2, treningowe$RH)

Regresja wieloraka

Poniżej widoczny jest wykres zależności modelu od wszystkich predyktorów. Można na nim zauważyć dobre wpasowanie danych w model. Współczynnik determinacji wynosi 0.54, czyli jest wikszy niż w przypadku regresji liniowej prostej.

rw1=lm(formula = RH~ ., data = treningowe[,-1], na.action = na.exclude)
summary(rw1)

## 
## Call:
## lm(formula = RH ~ ., data = treningowe[, -1], na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.221  -2.648   0.997   3.810  10.587 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.858e+01  1.082e+02   0.911 0.363679    
## ws          -1.729e+00  7.845e-01  -2.203 0.028883 *  
## wd          -1.049e-04  3.781e-03  -0.028 0.977906    
## air_temp    -5.399e-01  8.096e-02  -6.668 3.31e-10 ***
## atmos_pres  -4.898e-03  1.058e-01  -0.046 0.963140    
## visibility  -3.043e-04  8.075e-05  -3.769 0.000224 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.261 on 174 degrees of freedom
## Multiple R-squared:  0.5394, Adjusted R-squared:  0.5262 
## F-statistic: 40.75 on 5 and 174 DF,  p-value: < 2.2e-16

mw1=fitted(rw1)
plot(RH,mw1, main = "Wykres zaleznosci modelu od RH",
     ylab = "Model")

plot(rw1)

Sprawdzanie założeń

Jak widać poniżej suma oraz średnia są bardzo zbliżone do wartości 0.

sum(rw1$residuals)

## [1] 6.27276e-15

mean(rw1$residuals)

## [1] 3.504141e-17

shapiro.test(rw1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  rw1$residuals
## W = 0.96449, p-value = 0.0001554

jarque.test(rw1$residuals)

## 
##  Jarque-Bera Normality Test
## 
## data:  rw1$residuals
## JB = 16.518, p-value = 0.0002589
## alternative hypothesis: greater

hist(rw1$residuals, main = "Histogram reszt dla modelu regresji wielorakiej",
     xlab = "Reszty", ylab = "Czestotliwosc")

dwtest(rw1)

## 
##  Durbin-Watson test
## 
## data:  rw1
## DW = 0.98133, p-value = 1.069e-12
## alternative hypothesis: true autocorrelation is greater than 0

acf(rw1$residuals,lag.max = 100)

gvlma(rw1)

## 
## Call:
## lm(formula = RH ~ ., data = treningowe[, -1], na.action = na.exclude)
## 
## Coefficients:
## (Intercept)           ws           wd     air_temp   atmos_pres   visibility  
##  98.5798957   -1.7285978   -0.0001049   -0.5398595   -0.0048979   -0.0003043  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = rw1) 
## 
##                     Value   p-value                   Decision
## Global Stat        19.176 7.258e-04 Assumptions NOT satisfied!
## Skewness           15.386 8.765e-05 Assumptions NOT satisfied!
## Kurtosis            1.132 2.873e-01    Assumptions acceptable.
## Link Function       0.258 6.115e-01    Assumptions acceptable.
## Heteroscedasticity  2.400 1.213e-01    Assumptions acceptable.

Wykresy diagnostyczne dla modelu

ols_plot_diagnostics(rw1)

Wpasowanie danych w model

Jak widać na poniższym wykresie model stosunkowo dobrze wpasował dane. Ekstrema znajdują się w tych samych przedziałach czasowych, model w dużej czesci dobrze je dopasował.

ggof(mw1, treningowe$RH)

Regresja wieloraka z interakcjami

Poniżej widoczny jest wykres zależności modelu od wszystkich predyktorów z wszystkimi ich interakcjami. Można na nim zauważyć dobre wpasowanie danych w model. Współczynnik determinacji wynosi 0.66.

rw2=lm(RH~.*.,data=treningowe)
summary(rw2)

## 
## Call:
## lm(formula = RH ~ . * ., data = treningowe)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.3271  -2.8678   0.6078   3.4003   9.6508 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            5.342e+02  1.235e+03   0.433   0.6658  
## date                  -6.156e-07  1.099e-06  -0.560   0.5764  
## ws                     8.089e+01  1.916e+02   0.422   0.6735  
## wd                     1.154e+00  1.157e+00   0.997   0.3202  
## air_temp              -7.586e-02  2.166e+01  -0.004   0.9972  
## atmos_pres            -4.206e-01  1.206e+00  -0.349   0.7277  
## visibility            -5.314e-03  3.463e-02  -0.153   0.8783  
## date:ws                4.390e-10  7.015e-09   0.063   0.9502  
## date:wd               -2.242e-11  3.575e-11  -0.627   0.5315  
## date:air_temp          4.083e-10  7.858e-10   0.520   0.6040  
## date:atmos_pres        5.974e-10  1.075e-09   0.556   0.5792  
## date:visibility        2.274e-12  1.204e-12   1.888   0.0608 .
## ws:wd                  7.357e-03  8.624e-03   0.853   0.3949  
## ws:air_temp           -2.062e-01  1.550e-01  -1.330   0.1854  
## ws:atmos_pres         -8.048e-02  1.896e-01  -0.424   0.6718  
## ws:visibility         -1.637e-04  2.054e-04  -0.797   0.4265  
## wd:air_temp            2.336e-04  8.579e-04   0.272   0.7858  
## wd:atmos_pres         -1.132e-03  1.135e-03  -0.998   0.3200  
## wd:visibility          5.382e-07  1.028e-06   0.524   0.6013  
## air_temp:atmos_pres   -2.732e-04  2.131e-02  -0.013   0.9898  
## air_temp:visibility    7.462e-07  1.268e-05   0.059   0.9532  
## atmos_pres:visibility  1.750e-06  3.418e-05   0.051   0.9592  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.774 on 158 degrees of freedom
## Multiple R-squared:  0.6555, Adjusted R-squared:  0.6097 
## F-statistic: 14.32 on 21 and 158 DF,  p-value: < 2.2e-16

mw2=predict(rw2)

plot(RH,mw2, main = "Wykres zaleznosci modelu od RH",
     ylab = "Model")

plot(rw2)

Sprawdzanie założeń

Jak widać poniżej suma oraz średnia są bardzo zbliżone do wartości 0.

sum(rw2$residuals)

## [1] 1.387779e-15

mean(rw2$residuals)

## [1] 8.023096e-18

shapiro.test(rw2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  rw2$residuals
## W = 0.97496, p-value = 0.002502

jarque.test(rw2$residuals)

## 
##  Jarque-Bera Normality Test
## 
## data:  rw2$residuals
## JB = 10.858, p-value = 0.004387
## alternative hypothesis: greater

hist(rw2$residuals, main = "Histogram reszt dla modelu regresji wielorakiej z interakcjami",
     xlab = "Reszty", ylab = "Czestotliwosc")

dwtest(rw2)

## 
##  Durbin-Watson test
## 
## data:  rw2
## DW = 1.2152, p-value = 5.525e-09
## alternative hypothesis: true autocorrelation is greater than 0

acf(rw1$residuals,lag.max = 100)

Za pomocą funkcji bptest zbadano homoskedastyczność wariancji (wewnątrz wariancji wszystkie błędy są równe). Jak widac ponizej, p-value jest większe od 0.05, więc warunek jest spełniony.

bptest(rw2)

## 
##  studentized Breusch-Pagan test
## 
## data:  rw2
## BP = 18.639, df = 21, p-value = 0.6083

Wykresy diagnostyczne dla modelu

ols_plot_diagnostics(rw2)

Wpasowanie danych w model

Jak widać na poniższym wykresie model stosunkowo dobrze wpasował dane. Ekstrema znajdują się w tych samych przedziałach czasowych, model w dużej części dobrze je dopasował.

ggof(mw2, treningowe$RH)

Regresja krokowa LEAPS

Z poniższego wykresu wynika, że najlepszym parametrem jest temperatura powietrza, ponieważ występuje ona w największej liczbie modeli o największym R2. Duży wpływ na model ma też widoczność. Do utworzenia modelu wykorzystano temperaturę powietrza oraz widoczność, ponieważ te dwa predyktory dają wartość R2 równą 0.53. Dodanie innych parametrów tylko nieznacznie zwiększyłoby tę wartość.

rLeaps=regsubsets(RH~., data=treningowe[,-1], nbest=5)
plot(rLeaps,scale="r2")

m_Leaps=lm(RH~air_temp+visibility, data=treningowe, na.action = na.exclude)
mLeaps = predict(m_Leaps)
plot(RH,mLeaps, main = "Wykres zaleznosci modelu od RH",
     ylab = "Model")

Sprawdzanie założeń

Jak widać poniżej suma oraz średnia są bardzo zbliżone do wartości 0.

sum(m_Leaps$residuals)

## [1] 5.467848e-15

mean(m_Leaps$residuals)

## [1] 3.162015e-17

Za pomocą testu Shapiro-Wilka sprawdzono, czy reszty mają rozkład normalny. W tym przypadku wartość p-value jest mniejsza od 0.05, więc warunek nie został spełniony. Hipotezę alternatywną potwierdził również test Jarque-Bera.

shapiro.test(m_Leaps$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  m_Leaps$residuals
## W = 0.95747, p-value = 2.921e-05

jarque.test(m_Leaps$residuals)

## 
##  Jarque-Bera Normality Test
## 
## data:  m_Leaps$residuals
## JB = 18.126, p-value = 0.0001159
## alternative hypothesis: greater

hist(m_Leaps$residuals, main = "Histogram reszt dla modelu krokowej LEAPS",
     xlab = "Reszty", ylab = "Czestotliwosc")

Za pomocą testu Durbin-Watson zbadano autokorelację reszt. W przypadku badanego modelu odrzucono hipotezę zerową, która twierdzi, że korelacja między resztami nie występuje i przyjęto hipotezę alternatywną. Potwierdza to poniższy wykres.

dwtest(m_Leaps)

## 
##  Durbin-Watson test
## 
## data:  m_Leaps
## DW = 0.86742, p-value = 4.197e-15
## alternative hypothesis: true autocorrelation is greater than 0

acf(m_Leaps$residuals,lag.max = 100)

gvlma(m_Leaps)

## 
## Call:
## lm(formula = RH ~ air_temp + visibility, data = treningowe, na.action = na.exclude)
## 
## Coefficients:
## (Intercept)     air_temp   visibility  
##  88.7631390   -0.4480569   -0.0003735  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = m_Leaps) 
## 
##                      Value   p-value                   Decision
## Global Stat        24.0243 7.898e-05 Assumptions NOT satisfied!
## Skewness           17.3452 3.117e-05 Assumptions NOT satisfied!
## Kurtosis            0.7806 3.770e-01    Assumptions acceptable.
## Link Function       3.6238 5.696e-02    Assumptions acceptable.
## Heteroscedasticity  2.2747 1.315e-01    Assumptions acceptable.

Wykresy diagnostyczne dla modelu

ols_plot_diagnostics(m_Leaps)

Wpasowanie danych w model

Jak widać na poniższym wykresie model stosunkowo dobrze wpasował dane. Ekstrema znajdują się w tych samych przedziałach czasowych, ale model w dużej części ich nie doszacował.

ggof(mLeaps, treningowe$RH)

Dodatkowa metoda - regresja dynamiczna

Jako dodatkową metodę wybrano regresję dynamiczną - model ARIMA + regresja. Z raportu wynika, że wartość kryterium AIC wynosi 905, co oznacza, że jest to najniższa wartość z wszystkich testowanych modeli. Poniższy Time plot pokazuje, że model bardzo dobrze dopasował dane. Na poniższym wykresie rozrzutu także widać bardzo dobre wpasowanie danych w model.

m_dyn <- tstren %>%
  model(ARIMA(RH ~ visibility + air_temp + ws + wd + atmos_pres)) %>% 
  report()

## Series: RH 
## Model: LM w/ ARIMA(1,0,0)(0,1,2)[12] errors 
## 
## Coefficients:
##          ar1     sma1    sma2  visibility  air_temp       ws      wd
##       0.2845  -0.9365  0.0893      -4e-04   -0.4337  -0.1735  0.0017
## s.e.  0.0832   0.0887  0.0938       1e-04    0.1392   0.5143  0.0021
##       atmos_pres  intercept
##          -0.2883     0.2899
## s.e.         NaN     0.1085
## 
## sigma^2 estimated as 10.92:  log likelihood=-442.69
## AIC=905.38   AICc=906.78   BIC=936.62

dyn_dop <- bind_rows(m_dyn %>% 
            augment() %>% 
            mutate(nazwa = "m_dyn"))

dyn_dop %>% 
  ggplot(aes(x = date)) +
  geom_line(aes(y = .fitted, color = "model")) +
  geom_line(aes(y =RH, color = "dane")) +
  labs(color = "Reprezentacja", x = "Rok", y = "fitted", 
       title = "Timpe plot")

dyn_dop %>% 
  ggplot(aes(x = RH, y = .fitted, color = nazwa)) +
  ggtitle("Wykres rozrzutu (PM10 vs fitted)") +
  geom_point() +
  geom_abline(intercept = 0, slope = 1)

Na podstawie poniższej tabeli można zauważyć, że powyższy model ma bardzo wysoki współczynnik korelacji, a średnia oraz mediana reszt są zblizone do 0, co oznacza, że model regresji dynamicznej jest bardzo dobrym modelem.

dyn_dop %>% 
  as.data.frame() %>% 
  group_by(nazwa) %>% 
  summarise(kor = cor(RH, .fitted), 
            res_mean = mean(.resid),
            res_med = median(.resid)) %>% 
  knitr::kable()

nazwa	kor	res_mean	res_med
m_dyn	0.915657	0.1151379	0.3082575

Na podstawie poniższego wykresu reszt można zauważyć, że reszty nie są autoskorelowane, co potwierdził test autokorelacji, a ich rozkład nie jest normalny, co potwierdził test Shapiro-Wilka (p_value = 0.02).

m_dyn %>% gg_tsresiduals() + ggtitle("Wykres reszt dla modelu regresji dynamicznej")

shapiro.test(m_dyn%>% residuals() %>% pull(.resid))

## 
##  Shapiro-Wilk normality test
## 
## data:  m_dyn %>% residuals() %>% pull(.resid)
## W = 0.98144, p-value = 0.01694

m_dyn %>% augment() %>% features(.resid, ljung_box, lag = 20, dof = 2)

## # A tibble: 1 x 3
##   .model                                                   lb_stat lb_pvalue
##   <chr>                                                      <dbl>     <dbl>
## 1 ARIMA(RH ~ visibility + air_temp + ws + wd + atmos_pres)    27.9    0.0637

Prognoza

Przeprowadzono prognozę dla modelu regresji dynamicznej na danych testowych. Poniższy wykres pokazuje, że model bardzo dobrze wpasował dane na poziomie istotności 95%. Poniżej widoczny jest również raport błędów. Są one niewielkie, co ponownie sugeruje, że jest to bardzo dobry model.

forecast <- m_dyn %>% forecast(tstest)
dop_dyn <- bind_rows(accuracy(m_dyn),
            forecast %>% accuracy(tstest))
dop_dyn

## # A tibble: 2 x 10
##   .model         .type     ME  RMSE   MAE     MPE  MAPE    MASE   RMSSE     ACF1
##   <chr>          <chr>  <dbl> <dbl> <dbl>   <dbl> <dbl>   <dbl>   <dbl>    <dbl>
## 1 ARIMA(RH ~ vi~ Trai~  0.115  3.11  2.35  0.0405  3.00   0.568   0.568 -0.00732
## 2 ARIMA(RH ~ vi~ Test  -3.41   4.40  3.63 -4.74    5.01 NaN     NaN     -0.110

forecast %>% 
  autoplot(tstren, level = 95) + xlab("Year") +
  ggtitle("Prognoza dla modelu regresji dynamicznej")

Podsumowanie

tab

## # A tibble: 1 x 2
##   model AIC   
##   <chr> <chr> 
## 1 m_dyn 905.38

datatable(AIC(r1,r2,rw1,rw2,m_Leaps))

Porównując modele za pomocą kryterium AIC, dla którego im mniejsza wartość tym lepszy jest model, zauważono, że nalepiej wypadły: model dynamiczny oraz model regresji wielorakiej z wszystkimi interakcjami. Najmniej dokładne modele to te oparte na jednym predyktorze. Model regresji krokowej LEAPS, mimo że brał pod uwagę tylko dwa parametry, jest niewiele gorszy od modelu regresji wielorakiej. Trzeba też zauważyć, że żaden model nie spełnił kryterium normalności reszt. Natomiast kryterium autokorelacji spełnił tylko model dynamiczny. Miał on rownież najwyższy współczynnik korelacji, dlatego można założyć, że jest on najlepszy sposród analizowanych modeli.

Regresja

Przetwarzanie Danych Środowiskowych

Patrycja Majdak

10 czerwiec, 2021

Pozyskanie danych

Uśrednianie danych i sprawdzenie kompletności

Analiza danych

Wykresy korelacji

Wykresy rozrzutu

Wykres czasowy

Regresja liniowa prosta

Zmienna prognozowana i temperatura powietrza

Wykres zaleznosci modelu od zmiennej prognozowanej

Analiza reszt

Wykresy diagnostyczne dla modelu

Wpasowanie danych w model

Zmienna prognozowana i widoczność

Wykres zaleznosci modelu od zmiennej prognozowanej

Analiza reszt

Wykresy diagnostyczne dla modelu

Wpasowanie danych w model

Regresja wieloraka

Sprawdzanie założeń

Wykresy diagnostyczne dla modelu

Wpasowanie danych w model

Regresja wieloraka z interakcjami

Sprawdzanie założeń

Wykresy diagnostyczne dla modelu

Wpasowanie danych w model

Regresja krokowa LEAPS

Sprawdzanie założeń

Wykresy diagnostyczne dla modelu

Wpasowanie danych w model

Dodatkowa metoda - regresja dynamiczna

Prognoza

Podsumowanie