PDS_projekt
Dominik Kapusta, Filip Szot
1 06 2021
Pozyskanie danych
Dane dotyczące jakości powietrza pobrano z bazy GIOŚ dla lat 2007-2019. Zostały one odpowiednio wyselekcjonowane dla województwa mazowieckiego z pomiarem 1-godzinnym. Po przeglądnięciu dostępnych stacji pomiarowych zdecydowano się na wybór MzWarWokalna (dla lat 2007-2014 pod nazwą MzWarszUrsynow). Następnie pobrano dane meteorologiczne z bazy NOAA. Na poniższej mapie znaleziono najbliższą do MzWarWokalna stację meteorologiczną - OKECIE o kodzie 123750-99999 i pobrano dla niej dane. Stacje mają różne UTC, więc dla stacji zawierające dane meteorologiczne dodano godzinę. Oba zbiory danych zostały połączone funkcją inner_join.
library(worldmet)
library(tidyverse)
library(giosimport)
library(devtools)
library(tidyr)
library(naniar)
library(tsibble)
library(imputeTS)
library(feasts)
library(openair)
library(GGally)
library(fpp3)
library(ggpubr)
library(gridExtra)
library(purrr)
library(gtools)
library(DT)
# kat_dost <- "C:/Users/domka/Desktop/AGH/IVsemestr/PDS/czescResz/duzy_projekt_R"
#
# meta <- gios_metadane(type = "stacje",
# download = T,
# path = kat_dost,
# mode = "wb")
#
# dplyr::glimpse(meta)
#
# gios_vis(data = meta %>% filter(is.na(data.zamkniecia)))
#
# stanowiska <- gios_metadane(type = "stanowiska",
# download = F,
# path = kat_dost,
# mode = "wb")
#
# dplyr::glimpse(stanowiska)
#
# identyfikacja <- stanowiska %>%
# pull(kod.stacji)
#
# gios_vis(data = meta %>%
# filter(is.na(data.zamkniecia),
# kod.stacji %in% identyfikacja))
#
# statystyki <- gios_metadane(type = "statystyki",
# download = F,
# path = kat_dost,
# mode = "wb")
#
# names(statystyki)
#
# zrodlo %>%
# knitr::kable()
#
# pliki_7_19 <- map2(.x = as.list(zrodlo[8:20,1]),
# .y = as.list(zrodlo[8:20,2]),
# .f = gios_download,
# path = kat_dost,
# mode = "wb")
#
# names(pliki_7_19) <- paste0("R",zrodlo[8:20,2])
# names(pliki_7_19)
#
# pliki_7_19["R2016"]
#
# meta_mazowieckie <- meta %>%
# filter(wojewodztwo == "MAZOWIECKIE") %>%
# filter(data.zamkniecia >= "2007-01-01" | is.na(data.zamkniecia))
#
# stanowiska_mazowieckie <- stanowiska %>%
# filter(wojewodztwo == "MAZOWIECKIE") %>%
# filter(data.zamkniecia >= "2007-01-01" | is.na(data.zamkniecia))
#
# dplyr::glimpse(stanowiska_mazowieckie)
#
# stanowiska_mazowieckie_1g <- stanowiska %>%
# filter(wojewodztwo == "MAZOWIECKIE") %>%
# filter(data.zamkniecia >= "2007-01-01" | is.na(data.zamkniecia)) %>%
# filter(czas.usredniania == "1-godzinny")
#
# lista_folderow <- c(2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019)
#
# pliki_7_19_1g <- map(.x = pliki_7_19,
# .f = ~ .[str_detect(., pattern = "1g")]) %>%
# flatten_chr()
#
# pliki_1g <- map_df(.x = pliki_7_19_1g,
# .f = gios_read,
# czas_mu = "1g",
# path = kat_dost)
#
# kody_stacji_mazowieckie <- unique(stanowiska_mazowieckie_1g$kod.stacji)
#
# dane_1g_MzWarWokalna <- pliki_1g%>%
# filter(kod == "MzWarWokalna")
#
# dane_1g_MzWarszUrsynow <- pliki_1g%>%
# filter(kod == "MzWarszUrsynow")
#
#
#
# dane_1g_MzWarWokalna <- dane_1g_MzWarWokalna %>%
# spread(key=sub, value=obs)
#
# dane_1g_MzWarszUrsynow <- dane_1g_MzWarszUrsynow[!dane_1g_MzWarszUrsynow$sub=="1g",]
#
# dane_1g_MzWarszUrsynow <- dane_1g_MzWarszUrsynow %>%
# spread(key=sub, value=obs)
#
# dane_1g_MzWarszUrsynow<-subset(dane_1g_MzWarszUrsynow,select=-c(NOx,Nox))
# dane_1g_MzWarWokalna<-subset(dane_1g_MzWarWokalna,select=-c(NOx,NO))
#
# d1 <- rbind(dane_1g_MzWarWokalna, dane_1g_MzWarszUrsynow)
#
# dane_war_meteo = importNOAA(code = "123750-99999", year = 2007:2019, hourly = TRUE, quiet = FALSE)
load("d1.rda")
load("d_meteo.rda")
getMeta(country="PL", end.year="current", returnMap=TRUE)Sprawdzenie kompletności danych
Na początku usunięto kolumny z niepotrzebnymi danymi oraz tymi, które zawierały widoczne braki danych.Za pomocą funkcji vis_miss z pakietu naniar stworzony został wykres przedstawiający braki danych dla wyselekconowanych kolumn i miejsca, w których się one znajdują. Usunięte zostaną kolumny, których liczba NA > 5%. Są to C6H6, O3, PM10, PM2.5, SO2, ceil_hgt, cl_1, cl, cl_1_height.
dane_of<-subset(dane,select=-c(atmos_pres,cl_2,cl_3,cl_2_height,cl_3_height,precip_12,precip_6,precip,station,kod,code,latitude,longitude,elev))
vis_miss(dane_of, warn_large_data = F)
Po usunięciu wybranych kolumn braki danych w całym zbiorze są na poziomie 1% oraz żadna kolumna nie posiada więcej niż 4% NA. Do prognozowania został wybrany dwutlenek azout(NO2), ponieważ jako jednyny posiada kompletność danych powyżej 5% ze wszystkich zanieczyszczeń.
dane_of<-subset(dane_of,select=-c(PM10,PM2.5,O3,C6H6,SO2,ceil_hgt,cl_1,cl,cl_1_height))
vis_miss(dane_of, warn_large_data = F)
Kolejnym krokiem jest pozbycie się NA z danych. W tym celu skorzystano z dopasowania sezonowego wraz z metodą LOCF.
dane_gotowe <- na.seadec(
dane_of,
algorithm = "interpolation",
find_frequency = TRUE,
maxgap = Inf,
)
vis_miss(dane_gotowe, warn_large_data = F)
Szeregi czasowe i ich wizualizacja
Na początku uśredniono dane do średniomiesięcznych i przekonwertowano je na obiekt typu tsibble, tworząc szereg czasowy.
tsdane_month <- dane_gotowe %>%
mutate(month = yearmonth(date)) %>%
group_by(month) %>%
summarise(NO2 = mean(NO2), ws=mean(ws), wd=mean(wd), air_temp=mean(air_temp), visibility=mean(visibility), dew_point=mean(dew_point),RH=mean(RH)) %>%
as_tsibble(key = NULL, index = month)Dla poniższego wykresu widoczna jest pewna sezonowość - wzrosty i spadki występują cyklicznie.
autoplot(tsdane_month, NO2) +
ggtitle("Stężenie NO2 w latach 2007-2019") +
ylab("Stężenie NO2") +
xlab("Rok")
Z poniższego wykresu można zauważyć wzrosty stężenia NO2 w okresach zimowych i spadki w sezonie letnim. Na tle innych Wyróżnia się wykres rok 2007 z wysokim spadkiem w maju i wyraźnym wzrostem w kolejnych miesiącach.
gg_season(tsdane_month, NO2, labels = "both") +
ylab("NO2") +
xlab("miesiac")+
ggtitle("Wykres sezonowy: Stezenie sredniomiesieczne NO2 w latach 2007-2019")
Z tego wykresu można wyczytać, że największe średnie stężenie NO2 w latach 2007-2019 występuje w lutym, natomiast najmniejsze w czerwcu. W okresach letnich widać delikatną tendencję wzrostową, natomiast w październiku, listopadzie i grudniu występuje tendencja spadkowa.
tsdane_month %>%
gg_subseries(y = NO2, period = "year") +
labs(y = "NO2", title = "Subseries: Stezenie sredniomiesieczne NO2 w latach 2007-2019")
Występują dane poza zakresem przerywanych linii, zatem seria prawdopodobnie nie jest białym szumem.
ACF(tsdane_month, lag_max = 12) %>% autoplot() + labs(title = "Bialy szum - seria danych tsdane_month")
Korelacja
Patrząc na korelacje poszczególnych składników, temperatura jest bardzo skorelowana z widocznościa, temperaturą punktu rosy oraz wilgotnością. Z tego powodu została ona usunięta z analizy regresji, ponieważ wysoka korealcja pomiędzy składnikami może zaburzyć poprawność tworzonych modeli.
Sprawdzanie liczby danych dla poszczególnych lat
Na początku sprawdzono liczbę pomiarów poszczególnych latach. Od razu można zauważyć, że w latach 2008, 2012 i 2016 jest o około 25 danych więcej niż w pozostałych. Jest to spowodowane tym, że te lata mają dodatkowo jeden dzień, bo są przestępne. Nie licząc tego przypdaku nie widać dużych różnic w liczbie danych.
dane_tt <- subset(dane_gotowe, select = -c(air_temp))
dane_tt$years <- format(dane_tt$date, "%Y")
dane_tt %>%
group_by(years) %>%
summarise(liczba_pomiarow = n()) -> dane_tt_summary
ggplot(dane_tt_summary, aes(x=years, y=liczba_pomiarow, color=liczba_pomiarow)) +
geom_point(size=3, show.legend = FALSE) +
theme_minimal() +
coord_flip() +
xlab("rok") +
ylab("liczba pomiarów")
Podział danych na treningowe i testowe
Dane podzielono na dwie grupy - treningowe w zakresie lat 2007-2018 oraz testowe - rok 2019.
Modele
Regresja liniowa prosta
Pierwszy model jaki zostanie stworzony będzie bazował na regresji liniowej prostej, czyli do modelu wykorzystany będzie parametr, który miał największą korelację z NO2 - temperatura_punktu_rosy. R2 dla tego modelu wynosi 0.27.
dane_train %>%
ggplot(aes(x = NO2, y = dew_point)) +
geom_point() +
geom_smooth(method = "lm") +
ggtitle("Regresja liniowe NO2~dew_point")
## Series: NO2
## Model: TSLM
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.5937 -2.8616 -0.1427 2.3361 11.2236
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.06818 0.43899 57.104 < 2e-16 ***
## dew_point -0.40139 0.05417 -7.409 1.04e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.161 on 142 degrees of freedom
## Multiple R-squared: 0.2788, Adjusted R-squared: 0.2737
## F-statistic: 54.9 on 1 and 142 DF, p-value: 1.0384e-11Regresja wieloraka
Kolejny model zostanie zbudowany z zastosowaniem regresji wielorakiej. Do tego modelu wykorzystane będą wszystkie parametry. R2 w tym przypadku wynosi 0.40. Można zauważyć, że dew_point i ws mają największy wpływ na model.
fit_multi <- dane_train %>%
model(TSLM(formula = NO2 ~ ws + wd + visibility + dew_point + RH))
fit_multi %>% report()## Series: NO2
## Model: TSLM
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.36147 -2.43577 -0.05819 1.29585 10.64222
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.548e+01 4.995e+00 7.104 5.90e-11 ***
## ws -3.786e+00 6.947e-01 -5.450 2.25e-07 ***
## wd -1.045e-02 1.246e-02 -0.838 0.403
## visibility 8.942e-05 7.964e-05 1.123 0.263
## dew_point -5.577e-01 7.030e-02 -7.934 6.55e-13 ***
## RH 5.863e-02 4.527e-02 1.295 0.197
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.777 on 138 degrees of freedom
## Multiple R-squared: 0.4224, Adjusted R-squared: 0.4015
## F-statistic: 20.19 on 5 and 138 DF, p-value: 4.4398e-15Model na podstawie sezonowości i trendu
Ten model jest zbudowany tylko na podstawie trendu i sezonowości na podstawie naszych danych. W tym przypadku R2 wynosi 0.31.
fit_season_trend <- dane_train %>%
model(TSLM(NO2 ~ trend() + season()))
fit_season_trend %>% report()## Series: NO2
## Model: TSLM
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.3701 -2.4357 -0.4431 2.1567 10.6740
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.476719 1.291940 18.172 < 2e-16 ***
## trend() 0.022157 0.008158 2.716 0.007498 **
## season()year2 2.395169 1.655539 1.447 0.150352
## season()year3 0.770742 1.655599 0.466 0.642320
## season()year4 -1.788297 1.655700 -1.080 0.282089
## season()year5 -5.340537 1.655840 -3.225 0.001589 **
## season()year6 -6.072993 1.656021 -3.667 0.000355 ***
## season()year7 -6.203109 1.656242 -3.745 0.000269 ***
## season()year8 -5.118895 1.656503 -3.090 0.002444 **
## season()year9 -1.866675 1.656805 -1.127 0.261941
## season()year10 -0.578194 1.657146 -0.349 0.727717
## season()year11 -0.377452 1.657527 -0.228 0.820219
## season()year12 0.063502 1.657949 0.038 0.969506
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.055 on 131 degrees of freedom
## Multiple R-squared: 0.368, Adjusted R-squared: 0.3101
## F-statistic: 6.358 on 12 and 131 DF, p-value: 8.5456e-09Model z serią Fouriera
Dla K=2 R2 ma największą wartość równą 0.33.
## Series: NO2
## Model: TSLM
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.9188 -2.2852 -0.4733 2.0088 10.5255
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.457670 0.668830 32.082 < 2e-16 ***
## trend() 0.022285 0.008009 2.783 0.00615 **
## fourier(K = 2)C1_12 3.637644 0.469570 7.747 1.83e-12 ***
## fourier(K = 2)S1_12 -0.422088 0.470452 -0.897 0.37118
## fourier(K = 2)C2_12 -0.799249 0.469570 -1.702 0.09099 .
## fourier(K = 2)S2_12 1.141618 0.469706 2.430 0.01636 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.984 on 138 degrees of freedom
## Multiple R-squared: 0.3575, Adjusted R-squared: 0.3342
## F-statistic: 15.36 on 5 and 138 DF, p-value: 5.4576e-12Wybór predyktorów
Do wyboru predyktorów wypisano wszystkie możliwe kombinacje parametrów i sprawdzono zbudowany z nich model. Najlepszą kombinacją okazał się model z temperaturą punktu rosy oraz prędkością wiatru z R2 na poziomie 0.40 oraz najniższą wartością CV i AIC.
map(.x = 1:5,
.f = ~gtools::combinations(5, .x, colnames(dane_train)[-c(1,2)], repeats.allowed = F)) -> zmienne
list(mod1 = NO2 ~ dew_point,
mod2 = NO2 ~ wd,
mod3 = NO2 ~ ws,
mod4 = NO2 ~ visibility,
mod5 = NO2 ~ RH,
mod6 = NO2 ~ dew_point + wd,
mod7 = NO2 ~ dew_point + ws,
mod8 = NO2 ~ dew_point + visibility,
mod9 = NO2 ~ wd + ws,
mod10 = NO2 ~ wd + visibility,
mod11 = NO2 ~ RH + visibility,
mod12 = NO2 ~ RH + wd,
mod13 = NO2 ~ RH + ws,
mod14 = NO2 ~ RH + dew_point,
mod15 = NO2 ~ ws + visibility,
mod16 = NO2 ~ ws + wd + visibility,
mod17 = NO2 ~ ws + wd + dew_point,
mod18 = NO2 ~ ws + wd + RH,
mod19 = NO2 ~ ws + visibility + dew_point,
mod20 = NO2 ~ ws + visibility + RH,
mod21 = NO2 ~ ws + dew_point + RH,
mod22 = NO2 ~ wd + visibility + dew_point,
mod23 = NO2 ~ wd + visibility + RH,
mod24 = NO2 ~ wd + dew_point + RH,
mod25 = NO2 ~ visibility + dew_point + RH,
mod26 = NO2 ~ ws + wd + visibility + dew_point,
mod27 = NO2 ~ ws + wd + dew_point + RH,
mod28 = NO2 ~ ws + wd + visibility + RH,
mod29 = NO2 ~ ws + dew_point + visibility + RH,
mod30 = NO2 ~ dew_point + wd + RH + visibility,
mod31 = NO2 ~ dew_point + wd + RH + visibility+ws) -> formy
map(formy, as.formula) -> formy
out <- map(.x = formy,
.f = ~model(dane_train,
TSLM(formula = .x)) %>%
glance()
) %>%
do.call(rbind, .) %>%
select(.model, adj_r_squared, CV, AIC, AICc, BIC)
out %>%
mutate(.model = formy) %>%
arrange(AIC) %>%
knitr::kable(digits = 2)| .model | adj_r_squared | CV | AIC | AICc | BIC |
|---|---|---|---|---|---|
| NO2 ~ dew_point + ws | 0.40 | 14.44 | 387.01 | 387.30 | 398.89 |
| NO2 ~ ws + dew_point + RH | 0.40 | 14.69 | 388.43 | 388.87 | 403.28 |
| NO2 ~ ws + wd + dew_point | 0.40 | 14.60 | 388.65 | 389.08 | 403.50 |
| NO2 ~ ws + visibility + dew_point | 0.40 | 14.54 | 388.74 | 389.17 | 403.59 |
| NO2 ~ ws + dew_point + visibility + RH | 0.40 | 14.76 | 389.34 | 389.95 | 407.16 |
| NO2 ~ ws + wd + dew_point + RH | 0.40 | 14.85 | 389.92 | 390.53 | 407.74 |
| NO2 ~ ws + wd + visibility + dew_point | 0.40 | 14.70 | 390.35 | 390.96 | 408.17 |
| NO2 ~ dew_point + wd + RH + visibility + ws | 0.40 | 14.90 | 390.61 | 391.43 | 411.40 |
| NO2 ~ dew_point + wd | 0.28 | 17.47 | 414.57 | 414.86 | 426.45 |
| NO2 ~ dew_point | 0.27 | 17.53 | 414.59 | 414.76 | 423.50 |
| NO2 ~ wd + dew_point + RH | 0.28 | 17.70 | 415.57 | 416.01 | 430.42 |
| NO2 ~ RH + dew_point | 0.27 | 17.80 | 415.96 | 416.25 | 427.84 |
| NO2 ~ wd + visibility + dew_point | 0.27 | 17.64 | 416.48 | 416.92 | 431.33 |
| NO2 ~ dew_point + visibility | 0.27 | 17.70 | 416.54 | 416.82 | 428.41 |
| NO2 ~ dew_point + wd + RH + visibility | 0.28 | 17.81 | 416.68 | 417.29 | 434.50 |
| NO2 ~ visibility + dew_point + RH | 0.27 | 17.95 | 417.42 | 417.85 | 432.27 |
| NO2 ~ ws + wd + RH | 0.13 | 21.09 | 441.88 | 442.31 | 456.73 |
| NO2 ~ RH + wd | 0.13 | 21.18 | 441.92 | 442.21 | 453.80 |
| NO2 ~ ws + wd + visibility + RH | 0.13 | 21.16 | 442.72 | 443.33 | 460.54 |
| NO2 ~ wd + visibility + RH | 0.13 | 21.31 | 443.21 | 443.64 | 458.05 |
| NO2 ~ ws + visibility + RH | 0.11 | 21.80 | 446.48 | 446.92 | 461.33 |
| NO2 ~ RH + ws | 0.10 | 21.92 | 446.99 | 447.27 | 458.87 |
| NO2 ~ RH | 0.09 | 22.02 | 447.18 | 447.35 | 456.09 |
| NO2 ~ RH + visibility | 0.09 | 22.02 | 447.36 | 447.65 | 459.24 |
| NO2 ~ wd + visibility | 0.09 | 21.87 | 447.53 | 447.82 | 459.41 |
| NO2 ~ ws + wd + visibility | 0.10 | 21.81 | 447.85 | 448.28 | 462.70 |
| NO2 ~ ws + visibility | 0.08 | 22.22 | 449.94 | 450.23 | 461.82 |
| NO2 ~ visibility | 0.07 | 22.33 | 450.00 | 450.17 | 458.91 |
| NO2 ~ wd | 0.04 | 23.22 | 455.31 | 455.48 | 464.22 |
| NO2 ~ wd + ws | 0.03 | 23.39 | 457.15 | 457.44 | 469.03 |
| NO2 ~ ws | -0.01 | 24.20 | 461.43 | 461.61 | 470.34 |
Regresja nieliniowa
Stworzony został model, w którym została zlogarytmowana zmienna prognozowana, a jako predyktory zostały wybrane temperatura punktu rosy i prędkość wiatru. W tym przypadku wartość R2 jest mniejsza niż w modelu bez użycia logarytmu. Patrząc na wykres NO2 w czasie trudne byłoby stworzenie w tym przypadku węzłów, ponieważ nie występują widoczne możliwe zmiany kierunku linii trendu.
fit_wykl_dew_point_ws <- dane_train %>%
model(TSLM(formula = log(NO2)~dew_point+ws))
fit_wykl_dew_point_ws %>% report()## Series: NO2
## Model: TSLM
## Transformation: log(NO2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.847191 -0.093199 0.002779 0.092154 0.439508
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.816899 0.118752 32.142 <2e-16 ***
## dew_point -0.024701 0.002623 -9.416 <2e-16 ***
## ws -0.163162 0.031198 -5.230 6e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1723 on 141 degrees of freedom
## Multiple R-squared: 0.3865, Adjusted R-squared: 0.3778
## F-statistic: 44.42 on 2 and 141 DF, p-value: 1.0959e-15Model ETS
Jest to model tworzony przy użyciu metody wygładzenia wykładniczego, gdzie prognozy są ważonymi średnimi z poprzednich obserwacji, a wagi maleją wykładniczo z wiekiem obserwacji. Model ten ma minimalizować AIC, które w tym przypadku wynosi 1129.
## Series: NO2
## Model: ETS(M,N,A)
## Smoothing parameters:
## alpha = 0.1905546
## gamma = 0.0001001481
##
## Initial states:
## l s1 s2 s3 s4 s5 s6 s7
## 21.0828 1.950899 2.212831 1.819557 0.1658485 -3.184585 -4.57458 -4.204506
## s8 s9 s10 s11 s12
## -3.421376 -0.06813557 2.870246 4.712573 1.721227
##
## sigma^2: 0.0309
##
## AIC AICc BIC
## 1128.759 1132.509 1173.307Model Arima
Model Arima ma na celu opisanie autokorelacji danych. AIC w tym przypadku wyniosło 832.
## Series: NO2
## Model: ARIMA(1,0,0)(1,0,0)[12] w/ mean
##
## Coefficients:
## ar1 sar1 constant
## 0.4429 0.1609 10.7801
## s.e. 0.0773 0.0918 0.3452
##
## sigma^2 estimated as 18.27: log likelihood=-412.25
## AIC=832.49 AICc=832.78 BIC=844.37Model regresji dynamicznej
Dla NO2 i dew_point widać pewne podobieństwa, więc model regresji dynamicznej został stworzonny na podstawie temperatury punktu rosy.
tsdane_monthk %>%
gather("var", "value", NO2, dew_point) %>%
ggplot(aes(x = month, y = value)) +
geom_line() +
facet_grid(vars(var), scales = "free_y") +
xlab("data") + ylab(NULL) +
ggtitle("Miesięczne zmiany na MzWarWokalna dla NO2 i dew_point")
## Series: NO2
## Model: LM w/ ARIMA(1,0,0) errors
##
## Coefficients:
## ar1 dew_point I(dew_point^2) intercept
## 0.3627 -0.4403 0.0046 24.9259
## s.e. 0.0789 0.0963 0.0085 0.6395
##
## sigma^2 estimated as 15.31: log likelihood=-398.82
## AIC=807.64 AICc=808.08 BIC=822.49Dynamiczna regresja harmoniczna
Model przetestowano dla K ={1,..,6}. Dla naszych danych najlepiej wypadły przypadeki K = 4 i K = 6, ponieważ osiągnęły one najlepsze wyniki w testach ME, RMSE itd. W tym przypadku do dalszych analiz wybierzemy model z K = 4.
fit_ar_fiur <- dane_tr %>%
model(
`K = 1` = ARIMA(log(NO2) ~ fourier(K = 1) + PDQ(0,0,0)),
`K = 2` = ARIMA(log(NO2) ~ fourier(K = 2) + PDQ(0,0,0)),
`K = 3` = ARIMA(log(NO2) ~ fourier(K = 3) + PDQ(0,0,0)),
`K = 4` = ARIMA(log(NO2) ~ fourier(K = 4) + PDQ(0,0,0)),
`K = 5` = ARIMA(log(NO2) ~ fourier(K = 5) + PDQ(0,0,0)),
`K = 6` = ARIMA(log(NO2) ~ fourier(K = 6) + PDQ(0,0,0))
)
fit_ar_fiur %>%
forecast(h = "1 year") %>%
autoplot(dane_tr, size = 1, level = 80) +
facet_wrap(vars(.model), ncol = 2) +
guides(colour = FALSE) +
geom_label(
aes(x = yearmonth("2019 Jan"), y = 36, label = paste0("AICc = ", format(AICc))),
data = glance(fit_ar_fiur)
)
Model ws, trend i season
Model ten został stworzony testując różne kombinacje. Prędkość wiatru to jeden z parametrów, który był słabo skorelowany ze stężeniem NO2, ale dobrze go wyznaczał przy testach. Sezonowość i trend dobrze rozbudowują ten model, a temperatura punktu rosy nie była tu potrzebna, ponieważ nie zwiększała znacząco R2.
fit_dodatkowy <- dane_train %>%
model(TSLM(formula = NO2~log(ws)+season()+trend()))
fit_dodatkowy %>% report()## Series: NO2
## Model: TSLM
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4046 -2.0509 -0.1337 1.7387 9.5047
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.997503 4.043400 11.376 < 2e-16 ***
## log(ws) -15.853491 2.727660 -5.812 4.51e-08 ***
## season()year2 1.332620 1.491864 0.893 0.373369
## season()year3 0.498761 1.481413 0.337 0.736902
## season()year4 -3.494139 1.509570 -2.315 0.022199 *
## season()year5 -8.114593 1.555904 -5.215 7.04e-07 ***
## season()year6 -8.846387 1.556023 -5.685 8.20e-08 ***
## season()year7 -10.120636 1.627394 -6.219 6.33e-09 ***
## season()year8 -9.994006 1.702454 -5.870 3.42e-08 ***
## season()year9 -5.609382 1.615629 -3.472 0.000702 ***
## season()year10 -3.435276 1.561454 -2.200 0.029570 *
## season()year11 -2.256307 1.517236 -1.487 0.139405
## season()year12 -1.290035 1.500952 -0.859 0.391659
## trend() 0.017210 0.007345 2.343 0.020645 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.627 on 130 degrees of freedom
## Multiple R-squared: 0.4984, Adjusted R-squared: 0.4482
## F-statistic: 9.936 on 13 and 130 DF, p-value: 3.3566e-14Ocena modeli
W modelu regresji wielorakiej widoczne jest najlepsze dopasowanie wykresu funkcji niż w innych modelach, ale wykorzystuje on dość dużo parametrów, a wiele innych modeli przy ich mniejszej ilości uzyskuje podobny wynik (np. temperatura puntu rosy + prędkość wiatru). “Log(ws) + trend + season” ma też bardzo dobre dopasowanie. Model bazujący tylko na sezonowości i trendzie oraz Fourier dobrze odzwierciedlają wzrosty i spadki, ale nie uwidaczniają różnice między kolejnymi latami. Modele z kategorii Arima dość dobrze dopasowują się do rzeczywistych wartości NO2. Model ETS przez większość lat przyjmuje podobny kształt wykresu. Różni się jednak znacząco przy latach granicznych.
bind_rows(fit %>%
augment() %>%
mutate(nazwa = "Liniowy"),
fit_multi %>%
augment() %>%
mutate(nazwa = "Wieloraki"),
fit_season_trend %>%
augment() %>%
mutate(nazwa = "Season + trend"),
fit_k2 %>%
augment() %>%
mutate(nazwa = "Fourier k2"),
fit_dew_point_ws %>%
augment() %>%
mutate(nazwa = "Dew_point + ws"),
fit_wykl_dew_point_ws %>%
augment() %>%
mutate(nazwa = "Log + Dew_point + ws"),
fit_ets %>%
augment() %>%
mutate(nazwa = "ETS"),
fit_arima %>%
augment() %>%
mutate(nazwa = "Arima"),
fit_ar_dew %>%
augment() %>%
mutate(nazwa = "Arima + dew_point"),
fit_ar_fiur %>%
augment() %>%
mutate(nazwa = "Arima + Fourier"),
fit_dodatkowy %>%
augment() %>%
mutate(nazwa = "Log(ws) + trend + season"),
) -> dopasowanie
dopasowanie %>%
ggplot(aes(x = month)) +
geom_line(aes(y = .fitted, color = "model")) +
geom_line(aes(y = NO2, color = "dane")) +
facet_wrap(~nazwa, ncol = 1) +
labs(color = "Reprezentacja")
Dla regresji wielorakiej punkty są trochę lepiej wpasowane do linii w porównaniu do reszty.
dopasowanie %>%
ggplot(aes(x = NO2, y = .fitted, color = nazwa)) +
geom_point() +
geom_abline(intercept = 0, slope = 1) +
facet_wrap(~nazwa, ncol = 1)
Patrząc na poniższą tabelę widzimy, że największe korelacje, czyli dopasowanie do wykresu NO2 mają takie modele jak: Log(ws) + trend + season, wieloraki, Log + Dew_point + ws, Dew_point + ws, Arima + Fourier.
dopasowanie %>%
as.data.frame() %>%
group_by(nazwa) %>%
summarise(kor = cor(NO2, .fitted),
res_mean = mean(.resid),
res_med = median(.resid))Modele, które posiadają rozkład normalny bazując na teście Shapiro-Wilka i wartości p-value to Arima + dew_point, ETS, Fourier k2, Season+ trend.
# Shapiro-Wilk test
normalnosc <- data.frame(model=c("Arima","Arima + dew_point", "Arima + Fourier", "Dew_point + ws", "ETS",
" Fourier k2", "Liniowy"," Log + Dew_point+ ws", "Season+ trend", "wieloraki", "Log(ws) + trend + season" ),
pvalue=c(0.003165, 0.2441, 0.04665, 0.0008828, 0.05943, 0.07144, 0.02056,
0.000003731, 0.09446, 0.0008925, 0.0006556))
datatable(normalnosc)r1 <- ggpubr::ggqqplot(fit %>% residuals() %>% pull(.resid)) +
ggtitle("Liniowy")
r2 <- ggpubr::ggqqplot(fit_multi %>% residuals() %>% pull(.resid)) +
ggtitle("Wieloraki")
r3 <- ggpubr::ggqqplot(fit_season_trend %>% residuals() %>% pull(.resid)) +
ggtitle("Season + trend")
r4 <- ggpubr::ggqqplot(fit_k2 %>% residuals() %>% pull(.resid)) +
ggtitle("Fourier k2")
r5 <- ggpubr::ggqqplot(fit_dew_point_ws %>% residuals() %>% pull(.resid)) +
ggtitle("Dew_point + ws")
r6 <- ggpubr::ggqqplot(fit_wykl_dew_point_ws %>% residuals() %>% pull(.resid)) +
ggtitle("Log + Dew_point + ws")
r7 <- ggpubr::ggqqplot(fit_ets %>% residuals() %>% pull(.resid)) +
ggtitle("ETS")
r8 <- ggpubr::ggqqplot(fit_arima %>% residuals() %>% pull(.resid)) +
ggtitle("Arima")
r9 <- ggpubr::ggqqplot(fit_ar_dew %>% residuals() %>% pull(.resid)) +
ggtitle("Arima + dew_point")
r10 <- ggpubr::ggqqplot(fit_ar_fiur_best %>% residuals() %>% pull(.resid)) +
ggtitle("Arima + Fourier")
r11 <- ggpubr::ggqqplot(fit_dodatkowy %>% residuals() %>% pull(.resid)) +
ggtitle("Log(ws) + trend + season")
grid.arrange(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11, ncol = 4)
t1 <- fit %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t2 <- fit_multi %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t3 <- fit_season_trend %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t4 <- fit_k2 %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t5 <- fit_dew_point_ws %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t6 <- fit_wykl_dew_point_ws %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t7 <- fit_ets %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t8 <- fit_arima %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t9 <- fit_ar_dew %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t10 <- fit_ar_fiur_best %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t11 <- fit_dodatkowy %>%
augment() %>%
features(.resid, ljung_box, lag = 20, dof = 2)
t <- data.frame(rbind(t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11))
datatable(t)Na wykresach reszt widzimy, że spośród najlepiej dopasowanych modeli, najmniejsze średnie reszt ma model Dew_point+ws. Arima + dew_point i ETS są to modele, które posiadają lekko mniejszą korelację niż pozostałe 4 natomiast ich średnie reszty są dużo miejsze.
dopasowanie %>%
ggplot(aes(x = .fitted, y = .resid, color = nazwa)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap(~nazwa, ncol = 1)
Największe R2 ma dew_point + ws i model ze wszystkimi parametrami, mają one również najmniejsze AIC i BIC. Są to zatem najlepsze modele pod względem tych kryteriów.
porownanie <- map(.x = list(fit,fit_multi,fit_season_trend,fit_k2,fit_dew_point_ws),
.f = ~glance(.x) %>%
select(adj_r_squared, CV, AIC, AICc, BIC)
) %>%
do.call(rbind, .)
porownanie <- cbind(c("Liniowy","Wieloraki","Season + Trend","Fourier k","dew_point + ws"),porownanie)
names(porownanie)[1] <- "Model"
datatable(porownanie)Biorąc pod uwagę te kryteria najlepiej wypadały model wieloraki, dew_point + ws (+ jego wersja z logartymem) i model dynamicznej regresji harmonicznej i log(ws)+season+trend. Wypadały one najlepiej w większości modeli.
stats <- bind_rows(
fit_arima %>% accuracy(),
fit_ets %>% accuracy(),
fit %>% accuracy(),
fit_multi %>% accuracy(),
fit_season_trend %>% accuracy(),
fit_k2 %>% accuracy(),
fit_dew_point_ws %>% accuracy(),
fit_wykl_dew_point_ws %>% accuracy(),
fit_ar_dew %>% accuracy(),
fit_ar_fiur_best %>% accuracy(),
fit_dodatkowy %>% accuracy()
)
stats <- stats %>% mutate_at(vars(RMSE, MAE, MPE, MAPE, MASE, RMSSE, ACF1), funs(round(., 3)))
datatable(stats)Wybór modeli
Do dalszej analizy wybrane zostały modele:- Log(ws) + season + trend
- ws + dew_point
- regresja wieloraka
- model dynamicznej regresji harmonicznej
Posiadały one najdokładniejsze dopasowania do rzeczywistego stężenia NO2 oraz najlepsze wyniki w testach.
Prognozy
Najlepiejszą prognozę przeprowadziły modele regresji wielorakiej i dew_point + ws. W tym przypadku różnica pomiędzy nimi jest trudna do zauważenia więc lepszym wyborem jest model wykorzystujący mniejszą liczbę parametrów. Dobrze sprawdził się też model bazujący na prędkości wiatru, sezonowości i trendzie. Lepiej przewidział wzrost na początku 2019r, ale w kolejnych miesiącach tego roku poradził sobie gorzej. Model dynamicznej regresji harmonicznej dobrze przwidział wzrosty i spadki, ale nie słabo dopasował wykres do rzeczywistego stężenia NO2. Wszystkie modele przewidziały kształt wykresu w zakresie błędu 90%.
prog1 <- fit_dodatkowy %>%
forecast(new_data = tsdane_monthk %>%
filter(year(month) >= 2019)) %>%
autoplot(tsdane_monthk %>%
filter(year(month) > 2006),
color = "red", size = 1, level = 90) +
ggtitle("Prognoza z wykorzystaniem modelu log(ws) + season + trend")
prog2 <- fit_dew_point_ws %>%
forecast(new_data = tsdane_monthk %>%
filter(year(month) >= 2019)) %>%
autoplot(tsdane_monthk %>%
filter(year(month) > 2006),
color = "blue", size = 1, level = 90) +
ggtitle("Prognoza z wykorzystaniem modelu dew_point + ws")
prog3 <- fit_multi %>%
forecast(new_data = tsdane_monthk %>%
filter(year(month) >= 2019)) %>%
autoplot(tsdane_monthk %>%
filter(year(month) > 2006),
color = "green", size = 1, level = 90) +
ggtitle("Prognoza z wykorzystaniem modelu regresji wielorakiej")
prog4 <- fit_ar_fiur_best %>%
forecast(new_data = tsdane_monthk %>%
filter(year(month) >= 2019)) %>%
autoplot(tsdane_monthk %>%
filter(year(month) > 2006),
color = "orange", size = 1, level = 90) +
ggtitle("Prognoza z wykorzystaniem modelu dynamicznej regresji harmonicznej")
grid.arrange(prog1, prog2, prog3, prog4, ncol = 1)
Wnioski
Najlepiejszą prognozę przeprowadziły modele regresji wielorakiej i dew_point + ws. W tym przypadku różnica pomiędzy nimi jest trudna do zauważenia więc lepszym wyborem jest model wykorzystujący mniejszą liczbę parametrów. Dobrze sprawdził się też model bazujący na prędkości wiatru, sezonowości i trendzie. Lepiej przewidział wzrost na początku 2019r, ale w kolejnych miesiącach tego roku poradził sobie gorzej. Model dynamicznej regresji harmonicznej dobrze przwidział wzrosty i spadki, ale nie słabo dopasował wykres do rzeczywistego stężenia NO2. Wszystkie modele przewidziały kształt wykresu w zakresie błędu 90%. Dodawanie serii Fouriera do modelu w większości przypadków poprawiała wartości R2 i inne statystyki. Dla naszych danych model Arima sprawdził się lepiej niż ETS, szczególnie przy uwzględnieniu go w dynamicznej regresji harmonicznej, która również wypadła dobrze na tle modeli z innych kategorii. Patrząc na wszystkie modele można zauważyć zależnoźć, że te które uzyskały hipotezę zerową w teście Shapiro-Wilka otrzymywały najgorsze wyniki w testach na dopasowanie modelu do wykresu stężenia NO2. Podsumowując najlepsze są testy oparte na parametrach z ewentualnym dodatkiem informacji o sezonowości bądź trendzie lub Arima z serią Fouriera, ale nie można wykluczać innych modeli, ponieważ wszystko zależy od badanego czynnika.