Cel projektu

Celem projektu było przeprowadzenie modelowania danych o stężeniu PM10 na podstawie danych meteorologicznych oraz wybór najbardziej optymalnego modelu.

Pozyskanie danych meteorologicznych

Do uzyskania danych meteorologicznych zastosowano funkcje importNOAA z pakietu `worldmet. Do wykonania analizy wybrano dane z jednej ze stacji w Warszawie - Okęcie, z lat 2004-2019. Poniżej widoczne jest wczytanie danych oraz zapisanie ich do pliku.

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

Pozyskanie danych o jakości powietrza

Dane o jakości powietrza zaimportowano Z Głównego Inspektoratu Ochrony środowiska za pomocą pakietu giosimport. Wybrano stację Warszawa-Tołstoja znajdującą się w północnej części Warszawy. Do wykonania prognozy wybrano substancję PM10 (pomiary 1-godzinne). Poniżej pokazane jest, jak zostały wyselekcjonowane dane oraz zapisane do pliku.

inp_pm10 <- map(.x = pliki_all,
                .f = ~.[str_detect(.,pattern = "PM10_1g")])

inp_pm10 <- flatten_chr(inp_pm10)

PM10 <- map_df(.x = inp_pm10,
               .f = gios_read,
               czas_mu = "1g",
               path = kat_dost)

pm10_kody <- gios_kody(data = PM10, meta=meta)


pm10_kody %>%
  filter(kod == "MzWarTolstoj") %>%
  select(kod,date,obs) -> pm10_war

colnames(pm10_war) <- c("kod", "data", "PM10")
save(pm10_war, file = "pm10_war.Rdata")

Przygotowanie danych

Łączenie danych

Wczytano zapisane dane, a następnie d danych meteorologicznych dodano jedną godzinę, aby obie zmienne były w tej samej strefie czasowej. Następnie za pomocą funckji full_join połączono oba zbiory danych po dacie. Poniżej widoczny jest fragment powstałej tabeli.Poza prognozowaną substancją, widoczne są również predyktory:

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

#setwd("C:/PDS/projekt_Rzeszut/rzeszut_projekt")
load("meteo_war.Rdata")
load("pm10_war.Rdata")
colnames(pm10_war)[2] <- "date"
meteo_war <- meta
attr(meteo_war$date, "tzone") <- "Etc/GMT-1"
dane <- full_join(meteo_war,pm10_war,by="date")
head(dane)

## # A tibble: 6 x 10
##   date                   ws    wd air_temp atmos_pres visibility dew_point    RH
##   <dttm>              <dbl> <dbl>    <dbl>      <dbl>      <dbl>     <dbl> <dbl>
## 1 2004-01-01 01:00:00  3.05  350.    -0.05      1019        2250     -0.8   94.8
## 2 2004-01-01 02:00:00  3.05  355.    -0.2       1019        5000     -0.95  94.8
## 3 2004-01-01 03:00:00  3.55  340     -0.75      1019.       4250     -1.15  97.2
## 4 2004-01-01 04:00:00  3.3   345.    -0.95      1018.       3500     -1.7   94.7
## 5 2004-01-01 05:00:00  3.8   345.    -1.05      1018        3750     -1.8   94.7
## 6 2004-01-01 06:00:00  3.3   355.    -1.15      1018.       2750     -1.9   94.7
## # … with 2 more variables: kod <chr>, PM10 <dbl>

Uśrednianie danych i sprawdzenie kompletności

Za pomocą funkcji timeAverage utworzono zmienne zawierające informacje o danych średnio-miesięcznych oraz średnio-dobowych.

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

Powyżej widoczny jest wykres szeregu czasowego dla danych miesięcznych.

timePlot(mydata = dane_day, pollutant = c("ws", "wd", "air_temp", "atmos_pres", "visibility", "dew_point", "RH", "PM10"), y.relation = "free")

Powyżej widoczny jest wykres szeregu czasowego dla danych dobowych. Jak widać, żadna substancja nie posiada dużych braków danych, więc możliwe było uwzględnienie ich w dalszej analizie.

Konwersja na obiekt tsibble

Za pomocą funkcji as_tsibble zamieniono połączone oraz uśrednione dane na obiekt tsibble, który przechowuje szeregi czasowe.

tsdane <- dane %>% 
  as_tsibble(index = date)

tsdane_day <- dane_day %>% 
  mutate(date = ymd(date)) %>% 
  as_tsibble(index = date)

tsdane_mies <- dane_mies %>% 
  mutate(date = yearmonth(date)) %>% 
  as_tsibble(index = date)

ts_gather <- gather(tsdane_mies, key = "obs", value = "value", c(ws, wd,air_temp, atmos_pres, visibility, dew_point, RH, PM10))

Analiza danych

Wykres czasowy dla danych uśrednionych miesięcznie

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

Na powyższym wykresie widoczna jest sezonowość,szczególnie w drugiej połowie badanego okresu. Można też zauważyć niewielki trend malejący.

Wykresy sezonowe

Ponizej widoczne sa wykresy sezonowosci dla danych miesiecznych. Mozna zauwazyc duze zroznicowanie dla poszczeoglnych lat, ale wszystkie lata charakteryzuja sie zmniejszeniem ilosc PM10 w powietrzu w okresie letnim.

gg_season(tsdane_mies, PM10, labels = "both") +
  xlab("miesiąc") +
  ylab("obs") +
  ggtitle("Wykres sezonowy dla prognozowanej zmiennej")

gg_season(tsdane_mies, PM10, polar = T) +
  xlab("miesiąc") +
  ylab("obs") +
  ggtitle("Wykres seoznowy dla prognozowanej zmiennej")

Wykres podserii

Aby dokladniej przyjrzec sie sezonowosci wykonano wykres podserii dla prognozowanej zmiennej, ktory rowniez obrazuje nizsze stezenie w okresie letnim.

gg_subseries(tsdane_mies, PM10) +
  xlab("rok") +
  ylab("PM10") +
  ggtitle("Wykres podserii dla PM10")

Wykresy rozrzutu

W celu wstepnego zbadania korelacji miedzy prognozowana zmienna a przykladowymi predyktorami utworzono wykresy rozrzutu. Na ponizszym wykresie widac dosc wysoka korelacje ujemna dla PM10 i temperatury powietrza (-0.51). W przypadku RH korelacja jest znacznie mniejsza (0.18).

ggplot(data = tsdane_mies %>% 
         as.data.frame(), 
       aes(PM10, air_temp)) +
  ylab("temperatura powietrza") +
  geom_point() +
  geom_smooth() +
  ggtitle("Wykres rozrzutu dla zmiennej prognozowanej i temperatury powietrza")

cor(tsdane_mies[,"air_temp"], tsdane_mies[,"PM10"], method = "pearson")

##                PM10
## air_temp -0.5063934

ggplot(data = tsdane_mies %>% 
         as.data.frame(), 
       aes(PM10, RH)) +
  ylab("wilgotność względna") +
  geom_point() +
  geom_smooth() +
  ggtitle("Wykres rozrzutu dla zmiennej prognozowanej i wilgotności względnej")

cor(tsdane_mies[,"RH"], tsdane_mies[,"PM10"], method = "pearson")

##         PM10
## RH 0.1798173

Zbadano rowniez sezonowosc predyktorow. Przedstawia ja ponizszy wykres. Najbardziej wyrazna sezonowosc widoczna jest dla temperatury, temperatury punktu rosy, wilgotnosci wzglednej oraz widocznosci.

ts_gather %>%
  ggplot(aes(x = date, y = value)) +
  geom_line() +
  facet_grid(vars(obs), scales = "free_y")

wykresy korelacji

Przed budowa modeli zbadano korelacje wszystkich zmiennych. Miedzy zmiennymi air_temp i dew_point wystepuje zbyt duza korelacja (0,97), wiec postanowiono wykluczyc temperature punktu rosy z dalszej analizy. Predyktory najbardziej skorelowane ze zmienna prognozowana to: temperatura powietrza, widocznosc oraz cisnienie atmosferyczne. Przeprowadzono rowniez transformacje stezenia PM10 (wykonano pierwiastkowanie i logarytmowanie). Po tych operacjach korelacja zmienila sie jednak nieznacznie.

tsdane_mies %>%
  GGally::ggpairs(columns = 2:9)

tsdane_mies %>% mutate(PM10 = log10(PM10)) %>% 
  GGally::ggpairs(columns = 2:9)

tsdane_mies %>% mutate(PM10 = sqrt(PM10)) %>% 
  GGally::ggpairs(columns = 2:9)

Wykres opoznien

Utworzono wykres opoznien dla danych srednio-miesiecznych. Zaleznosc jest najsilniejsza dla opoznienia lag12, co odzwierciedla silna sezonowosc.

tsdane_mies %>%
  gg_lag(PM10, geom = "point", lags = c(2,4,6,8,10,12)) +
  labs(x = "PM10")

Wykresy autokorelacji

Utworzono wykresy autokorelacji w celu zidentyfikowania, w ktorych zmiennych wystepuje bialy szum. Sa to: wd i atmos_pres. Pozostale zmienne sa autoskorelowane.

w1 <- tsdane_mies %>%
  ACF(ws, lag_max = 36) %>%
  autoplot() + labs(title="ws")

w2 <- tsdane_mies %>%
  ACF(wd, lag_max = 36) %>%
  autoplot() + labs(title="wd")

w3 <- tsdane_mies %>%
  ACF(air_temp, lag_max = 36) %>%
  autoplot() + labs(title="air_temp")

w4 <- tsdane_mies %>%
  ACF(atmos_pres, lag_max = 36) %>%
  autoplot() + labs(title="atmos_pres")

w5 <- tsdane_mies %>%
  ACF(visibility, lag_max = 36) %>%
  autoplot() + labs(title="visibility")

w6 <- tsdane_mies %>%
  ACF(dew_point, lag_max = 36) %>%
  autoplot() + labs(title="dew_point")

w7 <- tsdane_mies %>%
  ACF(RH, lag_max = 36) %>%
  autoplot() + labs(title="RH")

ggarrange(w1,w2,w3,w4,w5,w7, 
          ncol = 2, nrow = 3)

Modelowanie

Zestawy treningowe i testowe

Ze wszystkich dostepnych danych wydzielono zestaw testowy (ostatnie pol roku) i zestaw treningowy (pozostale dane).

testowe <- tsdane_mies[187:192,]
treningowe <- tsdane_mies[-c(187:192),]

Regresja liniowa prosta

Przetestowano 13 przykladowych modeli regresji liniowej prostej. Dokonano oceny modeli na podstawie parametrów:

adj_r_squared - wspolczynnik determinacji (im blizszy wartosci 1 tym lepszy model)
CV - kroswalidacja (im nizsza wartosc tym lepszy model)
AIC - kryterium informacyjne Akaike dla duzej proby (im nizsza wartosc tym lepszy model)
AICc - kryterium informacyjne Akaike dla malej proby (im nizsza wartosc tym lepszy model)
BIC - kryterium informacyjne Schwarza-Bayesa

list(m1 = PM10 ~ air_temp,
     m2 = PM10 ~ atmos_pres,
     m3 = PM10 ~ visibility,
     m4 = PM10 ~ RH,
     m5 = log10(PM10) ~ air_temp,
     m6 = log10(PM10) ~ atmos_pres,
     m7 = log10(PM10) ~ visibility,
     m8 = log10(PM10) ~ RH,
     m9 = sqrt(PM10) ~ air_temp,
     m10 = sqrt(PM10) ~ atmos_pres,
     m11 = sqrt(PM10) ~ visibility,
     m12 = sqrt(PM10) ~ RH,
     m13 = log10(PM10) ~ sqrt(visibility)) -> formy

map(formy, as.formula)  -> formy  

out <- map(.x = formy, 
           .f = ~model(treningowe, 
                       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
log10(PM10) ~ sqrt(visibility)	0.36	0.01	-930.52	-930.39	-920.84
log10(PM10) ~ visibility	0.35	0.01	-927.71	-927.57	-918.03
log10(PM10) ~ air_temp	0.26	0.01	-902.57	-902.44	-892.90
log10(PM10) ~ atmos_pres	0.07	0.01	-861.79	-861.66	-852.11
log10(PM10) ~ RH	0.03	0.01	-853.14	-853.01	-843.47
sqrt(PM10) ~ visibility	0.34	0.42	-160.48	-160.35	-150.80
sqrt(PM10) ~ air_temp	0.26	0.47	-139.05	-138.92	-129.37
sqrt(PM10) ~ atmos_pres	0.08	0.58	-98.83	-98.70	-89.15
sqrt(PM10) ~ RH	0.03	0.61	-89.25	-89.11	-79.57
PM10 ~ visibility	0.32	80.47	818.52	818.65	828.20
PM10 ~ air_temp	0.26	88.63	836.07	836.20	845.75
PM10 ~ atmos_pres	0.08	109.99	875.30	875.44	884.98
PM10 ~ RH	0.03	115.71	885.76	885.89	895.44

Na podstawie analizy powyzszych parametrow mozna stwierdzic, ze sama regresja liniowa nie jest wystarczajaca, zeby modelowac dane. Wybrano jednak najlepszy z modeli (m13) i jego zgodnosc z zalozeniami teoretycznymi.

Badanie dopasowania modelu

Na wykresie Time plot widac, ze ekstrema znajduja sie w tych samych przedzialach czasowych, ale model znaczaco ich nie doszacowal. Na wykresie rozrzutu (PM10 vs fitted) mozna zauwazyc slabe wpasowanie danych w model. Na podstawie wykresu reszt mozna wstepnie zauwazyc, ze reszty sa autoskorelowane,a ich rozklad jest zblizony do normalnego, co potwierdza test Shapiro-Wilka. Wykres rozrzutu (fitted vs residuals) pokazuje, że nie da sie porawic tego modelu, gdyz prosta na powyzszym wykresie jest rownolegla do osi x.

ml13 <- treningowe %>% 
  model(TSLM(formula = log10(PM10) ~ sqrt(visibility)))

bind_rows(ml13 %>% 
            augment() %>% 
            mutate(nazwa = "model")) -> dopas

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

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

#Wykres reszt
ml13 %>% gg_tsresiduals() +
  ggtitle("Wykres reszt")

#autokorelacja test
ml13 %>% 
  augment() %>% 
  features(.resid, ljung_box, lag = 20, dof = 2)

## # A tibble: 1 x 3
##   .model                                         lb_stat lb_pvalue
##   <chr>                                            <dbl>     <dbl>
## 1 TSLM(formula = log10(PM10) ~ sqrt(visibility))    86.6  5.81e-11

#Normalność testy (p < 0.05 brak normalnosci)
shapiro.test(ml13 %>% residuals() %>% pull(.resid))

## 
##  Shapiro-Wilk normality test
## 
## data:  ml13 %>% residuals() %>% pull(.resid)
## W = 0.99306, p-value = 0.5258

#wykres rozrzutu  
dopas %>% 
  ggplot(aes(x = .fitted, y = .resid, color = nazwa)) +
  ggtitle("Wykres rozrzutu (fitted vs rediduals)") +
  geom_point() +
  geom_smooth(method = "lm", color = "blue") +
  facet_wrap(~nazwa, ncol = 1)

Model mozna rowniez oceniac poprzez zbadanie korelacji oraz sredniej i mediany reszt. Mozna zauwazyc, ze dla wybranego modelu wystepuje korelacja, ale nie jest zbyt wysoka, a srednia i mediana reszt maja wartosci bliskie 0.

dopas %>% 
  as.data.frame() %>% 
  summarise(kor = cor(PM10, .fitted), 
            res_mean = mean(.resid),
            res_med = median(.resid)) %>% 
  knitr::kable()

kor	res_mean	res_med
0.5918274	0.8046324	-0.3882818

Regresja liniowa z trendem i sezonowoscia

Przetestowano 15 przykladowych modeli regresji liniowej prostej dodajac do nich trend i sezonowosc. Na podstawie regresji liniowej prostej zauwazono, ze najlepsze modele proste tworzy sie na podstawie predyktorow visibility oraz air_temp, dlatego wykorzystano je do tworzenia tych modeli. Dokonano oceny modeli na podstawie parametrów (adj_r_squared, CV, AIC, AICc, BIC).

list(m1 = PM10 ~ visibility + trend(),
     m2 = PM10 ~ visibility + season(),
     m3 = PM10 ~ visibility + season() + trend(),
     m4 = PM10 ~ air_temp + trend(),
     m5 = PM10 ~ air_temp + season(),
     m6 = PM10 ~ air_temp + season() + trend(),
     m7 = log10(PM10) ~ visibility + trend(),
     m8 = log10(PM10) ~ visibility + season(),
     m9 = log10(PM10) ~ visibility + season() + trend(),
     m10 = log10(PM10) ~ air_temp + trend(),
     m11 = log10(PM10) ~ air_temp + season(),
     m12 = log10(PM10) ~ air_temp + season() + trend(),
     m13 = log10(PM10) ~ sqrt(visibility) + trend(),
     m14 = log10(PM10) ~ sqrt(visibility) + season(),
     m15 = log10(PM10) ~ sqrt(visibility) + season() + trend()) -> formy1

map(formy1, as.formula)  -> formy1  

out1 <- map(.x = formy1, 
            .f = ~model(treningowe, 
                       TSLM(formula = .x)) %>% 
             glance()) %>% 
  do.call(rbind, .) %>%  
  select(.model, adj_r_squared, CV, AIC, AICc, BIC)  

out1 %>% 
  mutate(.model = formy1) %>% 
  arrange(AIC) %>% 
  knitr::kable(digits = 2)

.model	adj_r_squared	CV	AIC	AICc	BIC
log10(PM10) ~ sqrt(visibility) + season() + trend()	0.50	0.01	-965.90	-963.08	-917.51
log10(PM10) ~ visibility + season() + trend()	0.49	0.01	-962.44	-959.62	-914.06
log10(PM10) ~ sqrt(visibility) + season()	0.49	0.01	-960.87	-958.42	-915.71
log10(PM10) ~ visibility + season()	0.47	0.01	-954.50	-952.05	-909.34
log10(PM10) ~ air_temp + season() + trend()	0.46	0.01	-949.65	-946.82	-901.26
log10(PM10) ~ air_temp + trend()	0.40	0.01	-940.94	-940.72	-928.04
log10(PM10) ~ sqrt(visibility) + trend()	0.37	0.01	-931.16	-930.93	-918.25
log10(PM10) ~ visibility + trend()	0.36	0.01	-928.78	-928.56	-915.88
log10(PM10) ~ air_temp + season()	0.31	0.01	-905.08	-902.62	-859.92
PM10 ~ visibility + season() + trend()	0.47	67.43	784.79	787.61	833.17
PM10 ~ visibility + season()	0.44	70.13	792.17	794.62	837.33
PM10 ~ air_temp + season() + trend()	0.44	71.63	793.33	796.15	841.72
PM10 ~ air_temp + trend()	0.39	73.58	801.41	801.63	814.32
PM10 ~ visibility + trend()	0.33	80.13	817.82	818.04	830.72
PM10 ~ air_temp + season()	0.31	87.81	832.10	834.56	877.26

Na podstawie analizy powyzszych parametrow wybrano najlepsze modele - m15, m9, m14.

Badanie dopasowania modeli

Na wykresach Time plot widac, ze ekstrema znajduja sie w tych samych przedzialach czasowych i sa lepiej dopasowane niz w przypadku samej regresji liniowej. Na wykresach rozrzutu (PM10 vs fitted) mozna zauwazyc niezbyt dobre wpasowanie danych w modele.

m15st <- treningowe %>% 
  model(TSLM(formula = log10(PM10) ~ sqrt(visibility) + season() + trend()))
m9st <- treningowe %>% 
  model(TSLM(formula = log10(PM10) ~ visibility + season() + trend()))
m14st <- treningowe %>% 
  model(TSLM(formula = log10(PM10) ~ sqrt(visibility) + season()))


bind_rows(m15st %>% 
            augment() %>% 
            mutate(nazwa = "Spierwiastkowana widocznosc z trendem i sezonowoscia"),
          m9st %>% 
            augment() %>% 
            mutate(nazwa = "Widocznosc z trendem i sezonowoscia"),
          m14st %>% 
            augment() %>% 
            mutate(nazwa = "Spierwiastkowana widocznosc z sezonowoscia")
          ) -> dop_sez_trend


#time_plot = fitted vs orginal data (czy dobrze wpasowane)
dop_sez_trend %>% 
  ggplot(aes(x = date)) +
  geom_line(aes(y = .fitted, color = "model")) +
  geom_line(aes(y =PM10, color = "dane")) +
  labs(color = "Reprezentacja", x = "Rok", y = "fitted", 
       title = "Timpe plot") +
  facet_wrap(~nazwa, ncol = 1)

#wykres rozrzutu = fitted vs orginal data (czy dobrze wpasowane)
dop_sez_trend %>% 
  ggplot(aes(x = PM10, y = .fitted, color = nazwa)) +
  ggtitle("Wykres rozrzutu (PM10 vs fitted)") +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  facet_wrap(~nazwa, ncol = 1)

Na podstawie ponizszej tabeli mozna zauwazyc, ze wspolczynniki korelacji sa bardzo zblizone, ale srednia reszt jest najbardziej zblizona do 0 dla modelu ze spierwiastkowana widocznoscia, trendem i sezonowoscia.

#Korelacja (orginal vs fitted) oraz średnia i mediana (residuals)
dop_sez_trend %>% 
  as.data.frame() %>% 
  group_by(nazwa) %>% 
  summarise(kor = cor(PM10, .fitted), 
            res_mean = mean(.resid),
            res_med = median(.resid)) %>% 
  knitr::kable()

nazwa	kor	res_mean	res_med
Spierwiastkowana widocznosc z sezonowoscia	0.7170452	0.6146240	0.6683194
Spierwiastkowana widocznosc z trendem i sezonowoscia	0.7217515	0.5882843	0.7186077
Widocznosc z trendem i sezonowoscia	0.7126984	0.6002936	0.9237461

Na podstawie ponizszego wykresu reszt mozna wstepnie zauwazyc, ze reszty sa autoskorelowane, co potwierdzil test autokorelacji,a ich rozklad nie jest normalny, co potwierdzil test Shapiro-Wilka (p_value < 0,05).

m15st %>% gg_tsresiduals() + ggtitle("Wykres reszt dla modelu ze spierwiastkowana widocznoscia, trendem i sezonowoscia")

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

## 
##  Shapiro-Wilk normality test
## 
## data:  m15st %>% residuals() %>% pull(.resid)
## W = 0.98412, p-value = 0.03339

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

## # A tibble: 1 x 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 TSLM(formula = log10(PM10) ~ sqrt(visibility) + season() + …    49.5 0.0000907

m9st %>% gg_tsresiduals() + ggtitle("Wykres reszt dla modelu z widocznoscia, trendem i sezonowoscia")

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

## 
##  Shapiro-Wilk normality test
## 
## data:  m9st %>% residuals() %>% pull(.resid)
## W = 0.98403, p-value = 0.03244

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

## # A tibble: 1 x 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 TSLM(formula = log10(PM10) ~ visibility + season() + trend(…    46.9  0.000222

Na podstawie ponizszego wykresu reszt mozna wstepnie zauwazyc, ze reszty sa autoskorelowane, co potwierdzil test autokorelacji. W tym przypadku rozklad reszt jest normalny (p_value = 0.08).

m14st %>% gg_tsresiduals() + ggtitle("Wykres reszt dla modelu z widocznoscia i sezonowoscia")

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

## 
##  Shapiro-Wilk normality test
## 
## data:  m14st %>% residuals() %>% pull(.resid)
## W = 0.98689, p-value = 0.08176

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

## # A tibble: 1 x 3
##   .model                                                    lb_stat   lb_pvalue
##   <chr>                                                       <dbl>       <dbl>
## 1 TSLM(formula = log10(PM10) ~ sqrt(visibility) + season())    65.9 0.000000221

Wykresy rozrzutu (fitted vs residuals) pokazuja, że nie da sie poprawic tych modeli.

dop_sez_trend %>% 
  ggplot(aes(x = .fitted, y = .resid, color = nazwa)) +
  ggtitle("Wykres rozrzutu (fitted vs residuals)") +
  geom_point() +
  geom_smooth(method = "lm", color = "purple") +
  facet_wrap(~nazwa, ncol = 1)

Regresja wieloraka

Przetestowano 11 przykladowych modeli regresji liniowej wielorakiej dodajac do nich trend i sezonowosc. Dokonano oceny modeli na podstawie parametrów (adj_r_squared, CV, AIC, AICc, BIC).

list(m1 = log10(PM10) ~ sqrt(visibility) + air_temp + season() + trend(),
     m2 = log10(PM10) ~ visibility + air_temp + season() + trend(),
     m3 = log10(PM10) ~ visibility + atmos_pres + season(),
     m4 = log10(PM10) ~ air_temp + atmos_pres + season() + trend(),
     m5 = log10(PM10) ~ air_temp + atmos_pres + season(), 
     m6 = log10(PM10) ~ air_temp + atmos_pres +trend(),
     m7 = log10(PM10) ~ visibility + air_temp + atmos_pres,
     m8 = log10(PM10) ~ visibility + air_temp + season(),
     m9 = log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH,
     m10 = log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH + season() + trend(),
     m11 = log10(PM10) ~ visibility + ws + atmos_pres + RH + season() + trend()) -> formy2

test9 <- treningowe %>% 
  model(TSLM(log10(PM10) ~ visibility + air_temp + ws  + RH)) %>% 
  report()

## Series: PM10 
## Model: TSLM 
## Transformation: log10(PM10) 
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.2772904 -0.0399311 -0.0005332  0.0458016  0.1829817 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.527e+00  8.368e-02  30.193  < 2e-16 ***
## visibility  -1.125e-05  1.253e-06  -8.983 3.40e-16 ***
## air_temp    -8.636e-03  1.028e-03  -8.403 1.25e-14 ***
## ws          -5.072e-02  1.065e-02  -4.765 3.87e-06 ***
## RH          -5.970e-03  7.450e-04  -8.013 1.33e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0671 on 181 degrees of freedom
## Multiple R-squared: 0.571,   Adjusted R-squared: 0.5615
## F-statistic: 60.22 on 4 and 181 DF, p-value: < 2.22e-16

map(formy2, as.formula)  -> formy2  

out2 <- map(.x = formy2, 
           .f = ~model(treningowe, 
                       TSLM(formula = .x)) %>% 
             glance() 
           ) %>% 
  do.call(rbind, .) %>%  
  select(.model, adj_r_squared, CV, AIC, AICc, BIC)  

out2 %>% 
  mutate(.model = formy2) %>% 
  arrange(AIC) %>% 
  knitr::kable(digits = 2)

.model	adj_r_squared	CV	AIC	AICc	BIC
log10(PM10) ~ visibility + ws + atmos_pres + RH + season() + , trend()	0.65	0.00	-1026.22	-1022.12	-968.15
log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + , RH + season() + trend()	0.64	0.00	-1023.10	-1018.01	-958.58
log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + , RH	0.58	0.00	-1003.37	-1002.56	-977.57
log10(PM10) ~ visibility + atmos_pres + season()	0.53	0.01	-977.09	-974.26	-928.70
log10(PM10) ~ air_temp + atmos_pres + season() + trend()	0.51	0.01	-969.57	-966.35	-917.96
log10(PM10) ~ sqrt(visibility) + air_temp + season() + trend()	0.50	0.01	-964.07	-960.85	-912.46
log10(PM10) ~ visibility + air_temp + season() + trend()	0.49	0.01	-960.96	-957.75	-909.35
log10(PM10) ~ air_temp + atmos_pres + trend()	0.45	0.01	-955.42	-955.09	-939.29
log10(PM10) ~ visibility + air_temp + season()	0.47	0.01	-953.38	-950.56	-905.00
log10(PM10) ~ visibility + air_temp + atmos_pres	0.42	0.01	-947.61	-947.28	-931.49
log10(PM10) ~ air_temp + atmos_pres + season()	0.35	0.01	-916.47	-913.64	-868.08

Na podstawie analizy powyzszych parametrow wybrano najlepsze modele - m11, m10, m9

Badanie dopasowania modeli

Na wykresach Time plot widac, ze ekstrema znajduja sie w tych samych przedzialach czasowych i sa dopasowane do oryginalnych danych. Na wykresach rozrzutu (PM10 vs fitted) rowniez mozna zauwazyc dobre wpasowanie danych do modeli.

m11w <- treningowe %>% 
  model(TSLM(formula = log10(PM10) ~ visibility + ws + atmos_pres + RH + season() + trend()))
m10w <- treningowe %>% 
  model(TSLM(formula = log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH + season() + trend()))
m9w <- treningowe %>% 
  model(TSLM(formula = log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH))


bind_rows(m11w %>% 
            augment() %>% 
            mutate(nazwa = "m11"),
          m10w %>% 
            augment() %>% 
            mutate(nazwa = "m10"),
          m9w %>% 
            augment() %>% 
            mutate(nazwa = "m9")
          ) -> dop_wiel

dop_wiel %>% 
  ggplot(aes(x = date)) +
  geom_line(aes(y = .fitted, color = "model")) +
  geom_line(aes(y =PM10, color = "dane")) +
  labs(color = "Reprezentacja", x = "Rok", y = "fitted", 
       title = "Timpe plot") +
  facet_wrap(~nazwa, ncol = 1)

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

Na podstawie ponizszej tabeli mozna zauwazyc, ze najlepsze wspolczynniki korelacji maja modele m10 i m11, a ich srednie oraz mediany reszt sa do siebie bardzo podobne.

#Korelacja (orginal vs fitted) oraz średnia i mediana (residuals)
dop_wiel %>% 
  as.data.frame() %>% 
  group_by(nazwa) %>% 
  summarise(kor = cor(PM10, .fitted), 
            res_mean = mean(.resid),
            res_med = median(.resid)) %>% 
  knitr::kable()

nazwa	kor	res_mean	res_med
m10	0.8336962	0.4058311	0.7898891
m11	0.8325610	0.4088464	0.8640550
m9	0.7746714	0.5115516	-0.1954516

m11w %>% gg_tsresiduals() + ggtitle("Wykres reszt dla m11")

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

## 
##  Shapiro-Wilk normality test
## 
## data:  m11w %>% residuals() %>% pull(.resid)
## W = 0.95937, p-value = 3.359e-05

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

## # A tibble: 1 x 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 "TSLM(formula = log10(PM10) ~ visibility + ws + atmos_pres …    92.6  4.90e-12

m10w %>% gg_tsresiduals() + ggtitle("Wykres reszt dla m10")

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

## 
##  Shapiro-Wilk normality test
## 
## data:  m10w %>% residuals() %>% pull(.resid)
## W = 0.95916, p-value = 3.196e-05

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

## # A tibble: 1 x 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 "TSLM(formula = log10(PM10) ~ visibility + air_temp + ws + …    96.9  8.02e-13

Na podstawie ponizszego wykresu reszt mozna wstepnie zauwazyc, ze reszty sa autoskorelowane, co potwierdzil test autokorelacji,a ich rozklad jest zblizony do normalnego, jednak odrzucono ta hipoteze na podstawie testu Shapiro-Wilka (p_value = 0,03).

m9w %>% gg_tsresiduals() + ggtitle("Wykres reszt dla m9")

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

## 
##  Shapiro-Wilk normality test
## 
## data:  m9w %>% residuals() %>% pull(.resid)
## W = 0.98372, p-value = 0.02938

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

## # A tibble: 1 x 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 "TSLM(formula = log10(PM10) ~ visibility + air_temp + ws + …    92.9  4.27e-12

Wykresy rozrzutu (fitted vs residuals) pokazuja, że nie da sie poprawic tych modeli.

#wykres rozrzutu = fitted vs residuals data (brak korelacji)
dop_wiel %>% 
  ggplot(aes(x = .fitted, y = .resid, color = nazwa)) +
  ggtitle("Wykres rzorzutu (fitted vs residuals)") +
  geom_point() +
  geom_smooth(method = "lm", color = "purple") +
  facet_wrap(~nazwa, ncol = 1)

Model ETS i Arima

Przetestowano rowniez model wykladniczy ETS oraz model Arima. Mozna zauwazyc, ze przy obu modelach automatycznych wartosci AIC oraz AICc wychodza najwyzsze ze wszystkich dotychczas analizowanych modeli.Jednak dla Arimy te wartosci sa nizsze.

ets <- treningowe %>% 
  model(ETS(PM10 ~ error("M") + trend("N") + season("M"))) %>% 
  report()

## Series: PM10 
## Model: ETS(M,N,M) 
##   Smoothing parameters:
##     alpha = 0.1431043 
##     gamma = 0.0001000106 
## 
##   Initial states:
##        l       s1      s2       s3        s4        s5        s6        s7
##  52.2486 1.010341 1.01407 1.005153 0.9721095 0.8195656 0.7952248 0.8263654
##         s8      s9      s10     s11      s12
##  0.8663066 1.09019 1.229157 1.19842 1.173099
## 
##   sigma^2:  0.0268
## 
##      AIC     AICc      BIC 
## 1738.202 1741.026 1786.588

arima <-treningowe %>%
  model(arima = ARIMA(PM10, stepwise = F, approximation = F)) %>%  
  report()

## Series: PM10 
## Model: ARIMA(2,1,1)(2,0,0)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1    sar1    sar2
##       0.1518  0.2396  -0.9723  0.1997  0.2240
## s.e.  0.0806  0.0754   0.0188  0.0749  0.0754
## 
## sigma^2 estimated as 81.66:  log likelihood=-669
## AIC=1350.01   AICc=1350.48   BIC=1369.33

Badanie dopasowania modeli

Na wykresach Time plot widac, ze ekstrema znajduja sie w tych samych przedzialach czasowych. Jednak mozna stwierdzic, ze model ETS jest lepiej wpasowany do danych. Na wykresach rozrzutu (PM10 vs fitted) rowniez mozna zauwazyc lepsze wpasowanie ETS.

bind_rows(arima %>% 
            augment() %>% 
            mutate(nazwa = "arima"), 
            ets %>% 
            augment() %>% 
            mutate(nazwa = "ets")) -> ar_ets_dop

#time_plot = fitted vs orginal data (czy dobrze wpasowane)
ar_ets_dop %>% 
  ggplot(aes(x = date)) +
  geom_line(aes(y = .fitted, color = "model")) +
  geom_line(aes(y =PM10, color = "dane")) +
  labs(color = "Reprezentacja", x = "Rok", y = "fitted", 
       title = "Timpe plot") +
  facet_wrap(~nazwa, ncol = 1)

#wykres rozrzutu = fitted vs orginal data (czy dobrze wpasowane)
ar_ets_dop %>% 
  ggplot(aes(x = PM10, y = .fitted, color = nazwa)) +
  ggtitle("Wykres rozrzutu (PM10 vs fitted)") +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  facet_wrap(~nazwa, ncol = 1)

Na podstawie ponizszej tabeli mozna zauwazyc, ze model ETS ma lepszy wspolczynnik korelacji, srednia oraz mediane reszt.

#Korelacja (orginal vs fitted) oraz średnia i mediana (residuals)
ar_ets_dop %>% 
  as.data.frame() %>% 
  group_by(nazwa)%>% 
  summarise(kor = cor(PM10, .fitted), 
            res_mean = mean(.resid),
            res_med = median(.resid)) %>% 
  knitr::kable()

nazwa	kor	res_mean	res_med
arima	0.5745709	-0.6840190	-2.4278113
ets	0.6895522	-0.5250333	-0.2306371

arima %>% gg_tsresiduals() + ggtitle("Wykres reszt dla modelu Arima")

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

## 
##  Shapiro-Wilk normality test
## 
## data:  arima %>% residuals() %>% pull(.resid)
## W = 0.97428, p-value = 0.00165

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

## # A tibble: 1 x 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 arima     31.9    0.0224

Na podstawie ponizszego wykresu reszt mozna wstepnie zauwazyc, ze reszty sa autoskorelowane, co potwierdzil test autokorelacji,a ich rozklad jest normalny, co potwierdzil test Shapiro-Wilka (p_value = 0,91).

ets %>% gg_tsresiduals() + ggtitle("Wykres reszt dal modelu ETS")

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

## 
##  Shapiro-Wilk normality test
## 
## data:  ets %>% residuals() %>% pull(.resid)
## W = 0.99602, p-value = 0.91

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

## # A tibble: 1 x 3
##   .model                                                    lb_stat lb_pvalue
##   <chr>                                                       <dbl>     <dbl>
## 1 "ETS(PM10 ~ error(\"M\") + trend(\"N\") + season(\"M\"))"    36.0   0.00697

Regresja dynamiczna

Przetestowano 10 przykladowych modeli regresji dynamicznej (Arima + regresja). Dokonano oceny modeli na podstawie parametrów (adj_r_squared, CV, AIC, AICc, BIC).

list(m1 = log10(PM10) ~ sqrt(visibility) + air_temp,
     m2 = log10(PM10) ~ visibility + air_temp,
     m3 = log10(PM10) ~ visibility + atmos_pres,
     m4 = log10(PM10) ~ air_temp + atmos_pres,
     m5 = log10(PM10) ~ air_temp + atmos_pres + I(atmos_pres^2), 
     m6 = log10(PM10) ~ air_temp + atmos_pres + I(air_temp^2) ,
     m7 = log10(PM10) ~ visibility + air_temp + atmos_pres,
     m8 = log10(PM10) ~ visibility + air_temp ,
     m9 = log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH,
     m10 = log10(PM10) ~ visibility + ws + atmos_pres + RH) -> formy3

map(formy3, as.formula)  -> formy3  

out3 <- map(.x = formy3, 
           .f = ~model(treningowe, 
                       ARIMA(formula = .x)) %>% 
             glance()) %>% 
  do.call(rbind, .) %>%  
  select(.model,AIC, AICc, BIC)  

out3 %>% 
  mutate(.model = formy3) %>% 
  arrange(AIC) %>% 
  knitr::kable(digits = 2)

.model	AIC	AICc	BIC
log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + , RH	-519.18	-518.15	-490.20
log10(PM10) ~ visibility + ws + atmos_pres + RH	-490.11	-489.29	-464.30
log10(PM10) ~ visibility + atmos_pres	-443.29	-442.03	-411.03
log10(PM10) ~ visibility + air_temp + atmos_pres	-442.54	-441.90	-419.99
log10(PM10) ~ air_temp + atmos_pres	-438.69	-438.22	-419.37
log10(PM10) ~ air_temp + atmos_pres + I(air_temp^2)	-437.13	-436.50	-414.59
log10(PM10) ~ air_temp + atmos_pres + I(atmos_pres^2)	-436.86	-436.23	-414.32
log10(PM10) ~ visibility + air_temp	-420.16	-419.69	-400.84
log10(PM10) ~ visibility + air_temp	-420.16	-419.69	-400.84
log10(PM10) ~ sqrt(visibility) + air_temp	-417.53	-416.89	-394.98

Na podstawie analizy powyzszych parametrow wybrano najlepsze modele - m9 i m10.

Badanie dopasowania modeli

Na wykresach Time plot widac, ze ekstrema znajduja sie w tych samych przedzialach czasowych i sa bardzo dobrze dopasowane do oryginalnych danych. Na wykresach rozrzutu (PM10 vs fitted) rowniez mozna zauwazyc dobre wpasowanie danych do modeli.

m9_dyn <- treningowe %>%
  model(ARIMA(log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH)) %>% 
  report()

## Series: PM10 
## Model: LM w/ ARIMA(1,1,1) errors 
## Transformation: log10(PM10) 
## 
## Coefficients:
##          ar1      ma1  visibility  air_temp       ws     wd  atmos_pres
##       0.3716  -0.9214           0   -0.0098  -0.0721  1e-04      0.0035
## s.e.  0.0942   0.0455           0    0.0012   0.0088  1e-04      0.0002
##            RH
##       -0.0057
## s.e.   0.0007
## 
## sigma^2 estimated as 0.003333:  log likelihood=268.59
## AIC=-519.18   AICc=-518.15   BIC=-490.2

m10_dyn <- treningowe %>%
  model(ARIMA(log10(PM10) ~ visibility + ws + atmos_pres + RH)) %>% 
  report()

## Series: PM10 
## Model: LM w/ ARIMA(1,0,0)(2,0,0)[12] errors 
## Transformation: log10(PM10) 
## 
## Coefficients:
##          ar1    sar1    sar2  visibility       ws  atmos_pres       RH
##       0.5172  0.2218  0.2749           0  -0.0727      0.0024  -0.0045
## s.e.  0.0698  0.0655  0.0677         NaN   0.0074      0.0001   0.0009
## 
## sigma^2 estimated as 0.003932:  log likelihood=253.05
## AIC=-490.11   AICc=-489.29   BIC=-464.3

bind_rows(m9_dyn %>% 
            augment() %>% 
            mutate(nazwa = "m9"),
          m10_dyn %>% 
            augment() %>% 
            mutate(nazwa = "m10")) -> dynam_dop

dynam_dop %>% 
  ggplot(aes(x = date)) +
  geom_line(aes(y = .fitted, color = "model")) +
  geom_line(aes(y =PM10, color = "dane")) +
  labs(color = "Reprezentacja", x = "Rok", y = "fitted", 
       title = "Timpe plot") +
  facet_wrap(~nazwa, ncol = 1)

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

Na podstawie ponizszej tabeli mozna zauwazyc, ze oba modele maja podobne wspolczynniki korelacji oraz srednie reszt, jednak m9 wypada nieco korzystniej.

#Korelacja (orginal vs fitted) oraz średnia i mediana (residuals)
dynam_dop %>% 
  as.data.frame() %>% 
  group_by(nazwa) %>% 
  summarise(kor = cor(PM10, .fitted), 
            res_mean = mean(.resid),
            res_med = median(.resid)) %>% 
  knitr::kable()

nazwa	kor	res_mean	res_med
m10	0.7948263	0.3121170	0.5408346
m9	0.8311435	0.2160248	0.0821071

Na podstawie ponizszego wykresu reszt mozna zauwazyc, ze reszty nie sa autoskorelowane, co potwierdzil test autokorelacji,a ich rozklad jest normalny, co potwierdzil test Shapiro-Wilka (p_value = 0.95).

m9_dyn %>% gg_tsresiduals() + ggtitle("Wykres reszt dla modelu m9")

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

## 
##  Shapiro-Wilk normality test
## 
## data:  m9_dyn %>% residuals() %>% pull(.resid)
## W = 0.99646, p-value = 0.9462

m9_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(log10(PM10) ~ visibility + air_temp + ws + wd + atmo…    27.1    0.0776

Na podstawie ponizszego wykresu reszt mozna zauwazyc, ze reszty sa autoskorelowane, co potwierdzil test autokorelacji,a ich rozklad jest normalny, co potwierdzil test Shapiro-Wilka (p_value = 0.36).

m10_dyn %>% gg_tsresiduals() + ggtitle("Wykres reszt dla modelu m10")

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

## 
##  Shapiro-Wilk normality test
## 
## data:  m10_dyn %>% residuals() %>% pull(.resid)
## W = 0.99169, p-value = 0.3635

m10_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(log10(PM10) ~ visibility + ws + atmos_pres + RH)    50.0 0.0000748

Wykresy rozrzutu (fitted vs residuals) pokazuja, że nie da sie poprawic tych modeli.

#wykres rozrzutu = fitted vs residuals data (brak korelacji)
dynam_dop %>% 
  ggplot(aes(x = .fitted, y = .resid, color = nazwa)) +
  geom_point() +
  geom_smooth(method = "lm", color = "blue") +
  facet_wrap(~nazwa, ncol = 1)

Prognoza

Na podstawie powyzszych ocen modeli wybrano 6 tych, ktore uznano za najlepsze:

dynamiczne
- log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH (m9)
- log10(PM10) ~ visibility + ws + atmos_pres + RH (m10)
model ETS
regresji wielorakiej
- log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH + season() + trend() (m10)
- log10(PM10) ~ visibility + ws + atmos_pres + RH + season() + trend() (m11)
regresji liniowej z trendem i sezonowoscia
- log10(PM10) ~ sqrt(visibility) + season() + trend() (m15)

Przeprowadzono prognozy dla powyzszych modeli na danych testowych, ktore nie zostaly wykorzystane do budowy modelu. Ponizej przedstawione sa wykresy przedstawiajace prognozy dla poszczegolnych modeli.

forecast1 <- m15st %>% forecast(testowe)
forecast1 %>% 
  autoplot(tsdane_mies, level = 95) + xlab("Year") +
  ggtitle("Prognoza dla modelu regresji liniowej z trendem i sezonowoscia")

forecast2 <- m11w %>% forecast(testowe)
forecast2 %>% 
  autoplot(tsdane_mies, level = 95) + xlab("Year") +
  ggtitle("Prognoza dla modelu regresji wielorakiej - 
          log10(PM10) ~ visibility + ws + atmos_pres + RH + season() + trend()")

forecast3 <- m10w %>% forecast(testowe)
forecast3 %>% 
  autoplot(tsdane_mies, level = 95) + xlab("Year") +
  ggtitle("Prognoza dla modelu regresji wielorakiej -
          log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH + season() + trend()")

forecast4 <- ets %>% forecast(testowe)
forecast4 %>% 
  autoplot(tsdane_mies, level = 95) + xlab("Year") + 
  ggtitle("Prognoza dla modelu ETS")

forecast5 <- m9_dyn %>% forecast(testowe)
forecast5 %>% 
  autoplot(tsdane_mies, level = 95) + xlab("Year") +
  ggtitle("Prognoza dla modelu regresji dynamicznej -
          log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH")

forecast6 <- m10_dyn %>% forecast(testowe)
forecast6 %>% 
  autoplot(tsdane_mies, level = 95) + xlab("Year") +
  ggtitle("Prognoza dla modelu regresji dynamicznej - 
          log10(PM10) ~ visibility + ws + atmos_pres + RH")

dop_15st <- accuracy(m15st)
dop_11w <- accuracy(m11w)
dop_10w <- accuracy(m10w)
dop_ets <- accuracy(ets)
dop_9dyn <- accuracy(m9_dyn)
dop_10dyn <- accuracy(m10_dyn)

porownanie <- bind_rows(dop_15st ,
                        forecast1 %>% accuracy(testowe),
                        dop_11w,
                        forecast2 %>% accuracy(testowe),
                        dop_10w,
                        forecast3 %>% accuracy(testowe),
                        dop_ets,
                        forecast4 %>% accuracy(testowe),
                        dop_9dyn,
                        forecast5 %>% accuracy(testowe),
                        dop_10dyn,
                        forecast6 %>% accuracy(testowe))

porownanie

## # A tibble: 12 x 10
##    .model        .type     ME  RMSE   MAE     MPE  MAPE    MASE   RMSSE     ACF1
##    <chr>         <chr>  <dbl> <dbl> <dbl>   <dbl> <dbl>   <dbl>   <dbl>    <dbl>
##  1 "TSLM(formul… Trai…  0.588  7.52  5.79  -1.29   12.6   0.677   0.645  1.44e-1
##  2 "TSLM(formul… Test  -5.53   6.21  5.53 -18.0    18.0 NaN     NaN      1.47e-1
##  3 "TSLM(formul… Trai…  0.409  6.02  4.67  -0.918  10.5   0.547   0.516  3.83e-1
##  4 "TSLM(formul… Test  -5.00   5.55  5.00 -16.0    16.0 NaN     NaN     -1.57e-1
##  5 "TSLM(formul… Trai…  0.406  6.00  4.64  -0.914  10.4   0.543   0.515  3.95e-1
##  6 "TSLM(formul… Test  -5.06   5.53  5.06 -16.2    16.2 NaN     NaN     -1.67e-1
##  7 "ETS(PM10 ~ … Trai… -0.525  7.97  5.98  -3.57   13.4   0.699   0.684  6.67e-4
##  8 "ETS(PM10 ~ … Test  -7.12   8.11  7.12 -23.1    23.1 NaN     NaN      2.14e-1
##  9 "ARIMA(log10… Trai…  0.216  6.03  4.62  -1.13   10.2   0.541   0.517 -4.32e-2
## 10 "ARIMA(log10… Test  -7.37   7.83  7.37 -23.7    23.7 NaN     NaN     -3.23e-1
## 11 "ARIMA(log10… Trai…  0.312  6.58  5.18  -1.31   11.4   0.606   0.564 -8.73e-2
## 12 "ARIMA(log10… Test  -5.84   6.58  5.84 -18.6    18.6 NaN     NaN     -7.43e-1

Wnioski

Po analizie powyzszych wykresow oraz tabeli stwierdzono, ze najlepsze modele to obydwa modele regresji wielorakiej: log10(PM10) ~ visibility + air_temp + ws + wd + atmos_pres + RH + season() + trend() log10(PM10) ~ visibility + ws + atmos_pres + RH + season() + trend() Mimo, ze po analizie teoretycznej modele te wydawaly sie jednymi ze slabszych sposrod wybranych, to po wykonaniu prognoz okazaly sie najlepsze w praktyce. Mialy one najnizsze bledy RMSE oraz MAE.

Prognozowanie

Przetwarzanie Danych Środowiskowych

Patrycja Majdak

07 czerwiec, 2021

Cel projektu

Pozyskanie danych meteorologicznych

Pozyskanie danych o jakości powietrza

Przygotowanie danych

Łączenie danych

Uśrednianie danych i sprawdzenie kompletności

Konwersja na obiekt tsibble

Analiza danych

Wykres czasowy dla danych uśrednionych miesięcznie

Wykresy sezonowe

Wykres podserii

Wykresy rozrzutu

wykresy korelacji

Wykres opoznien

Wykresy autokorelacji

Modelowanie

Zestawy treningowe i testowe

Regresja liniowa prosta

Badanie dopasowania modelu

Regresja liniowa z trendem i sezonowoscia

Badanie dopasowania modeli

Regresja wieloraka

Badanie dopasowania modeli

Model ETS i Arima

Badanie dopasowania modeli

Regresja dynamiczna

Badanie dopasowania modeli

Prognoza

Wnioski