Przetwarzanie Danych Środowiskowych
Projekt regresja liniowa
Jakub Rudek, Geoinformacja, II rok
21 sierpień, 2020
1 Wprowadzenie
Celem projektu jest utworzenie modelu liniowego na przykładowych danych, diagnostyka i optymalizacja modelu.
Dane pobrano z UCI - repozytorium danych do uczenia maszynowego: https://archive.ics.uci.edu/ml/datasets/Seoul+Bike+Sharing+Demand
Dane dotyczą ilości wypożyczonych rowerów w Seulu i warunków meteorologicznych.
Plik w formacie csv zawiera ilość wypożyczonych rowerów w ciągu kolejnych godzin każdego dnia od grudnia 2017 roku do listopada 2018 roku.
2 Przygotowanie danych
Sprawdziłem, że największa ilość rowerów jest wypożyczana o godzinie 18. Do dalszej analizy zostały wybrane wiersze z danymi z tej godziny.
Po odrzuceniu kolumn z niepotrzebnymi do modelu regresji danymi, otrzymano następującą tabelę:
3 Sprawdzenie danych
3.1 Wykresy rozrzutu dla zmiennych
Większe wykresy rozrzutu zrobione “ręcznie”:
3.2 Testy czy mamy rozkłady normalne zmiennych
## [1] "liczba"
##
## Shapiro-Wilk normality test
##
## data: dane[, i]
## W = 0.93335, p-value = 1.044e-11
##
## [1] "temp"
##
## Shapiro-Wilk normality test
##
## data: dane[, i]
## W = 0.97071, p-value = 1.009e-06
##
## [1] "wilg"
##
## Shapiro-Wilk normality test
##
## data: dane[, i]
## W = 0.9514, p-value = 1.331e-09
##
## [1] "wiatr"
##
## Shapiro-Wilk normality test
##
## data: dane[, i]
## W = 0.97932, p-value = 4.296e-05
##
## [1] "wid"
##
## Shapiro-Wilk normality test
##
## data: dane[, i]
## W = 0.76676, p-value < 2.2e-16
##
## [1] "rosa"
##
## Shapiro-Wilk normality test
##
## data: dane[, i]
## W = 0.96362, p-value = 7.01e-08
##
## [1] "slonce"
##
## Shapiro-Wilk normality test
##
## data: dane[, i]
## W = 0.84975, p-value < 2.2e-16Sprawdzenie założenia 1 - rozkład w grupach jest zbliżony do normalnego. Dla każdej zmiennej p-value jest mniejsze od 0.05 czyli niestety nie mamy rozkładów normalnych zmiennych. Musimy się z tym pogodzić i wiedząc to robimy dalej regresję.
3.3 Test zależności między zmiennymi
ggplot(data = dane, aes(x = temp, y = liczba)) +
geom_point() +
geom_smooth(method = 'lm', formula = y~x-1) +
labs(title = "Zależność między temperaturą a liczbą rowerów",
caption = "Dane dla godziny 18:00 w poszczególnych dniach") +
xlab("temperatura [st.C]") +
scale_y_continuous(limits = c(0, 3600)) +
theme_bw()
ggplot(data = dane, aes(x = wilg, y = liczba)) +
geom_point() +
geom_smooth(method = 'lm', formula = y~x-1) +
labs(title = "Zależność między wilgotnością a liczbą rowerów",
caption = "Dane dla godziny 18:00 w poszczególnych dniach") +
xlab("wilgotnoć [%]") +
scale_y_continuous(limits = c(0, 3600)) +
theme_bw()
Sprawdzenie założenia 2 zależność liniowa między predyktorami a zmienną objśnianą. Istnieje zależność między zmienną objaśnianą - liczbą pożyczonych rowerów o godzinie 18, a zmiennymi objaśniającymi - warunkami meteorologicznymi - prawidłowo. Sprawdzenie założenia 3 - brak korelacji między danymi (zmiennymi objaśniającymi). Niestety występują też wzajemne korelacje pomiędzy zmiennymi objaśniającymi.
4 Model regresji wielorakiej
4.1 Utworzenie modelu liniowego regresji
##
## Call:
## lm(formula = liczba ~ ., data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2203.2 -392.2 -31.5 546.7 1548.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4543.51093 795.89733 5.709 2.39e-08 ***
## temp -84.88064 28.93074 -2.934 0.00356 **
## wilg -50.60599 8.92131 -5.672 2.91e-08 ***
## wiatr -31.95718 39.84549 -0.802 0.42307
## wid 0.08475 0.08168 1.038 0.30013
## rosa 144.76159 30.60393 4.730 3.23e-06 ***
## slonce 366.93582 144.22322 2.544 0.01137 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 699.5 on 358 degrees of freedom
## Multiple R-squared: 0.5458, Adjusted R-squared: 0.5382
## F-statistic: 71.7 on 6 and 358 DF, p-value: < 2.2e-16Wartość R2 informuje nas o dopasowaniu modelu. Tutaj model wyjaśnia około 54% danych.
4.2 Wykresy diagnostyczne dla modelu
Żeby sprawdzić, czy spełnione są założenia modelu analizujemy, jak kształtują sie reszty.
Odczytujemy to z czterech wykresów utworzonych funkcją plot(regresja).
Na pierwszym wykresie widać, że reszty mają średnią koło zera i niewiele zależą od wartości dopasowanych modelu – jest to dobra właściwość. Na wykresie kwantylów reszt w stosunku do rozkładu normalnego reszt punkty niezbyt układają się liniowo. Nie jest to zbyt dobra właściwość. Na trzecim wykresie na początku wariancja reszt nie jest jednorodna, co jest niekorzystne. Na czwartym wykresie „dźwigni” też na początku widać, że występują wartości nietypowe, „dźwigające” współczynniki modelu, więc nie jest to najlepiej. Na tej podstawie można stwierdzić, że ten model jest nezbyt dobry.
4.3 Model regresji wielorakiej z interakcjami
##
## Call:
## lm(formula = liczba ~ wiatr + wid + slonce * temp + rosa * wilg,
## data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2342.74 -331.84 -15.37 511.46 1505.37
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3509.36307 835.63344 4.200 3.38e-05 ***
## wiatr -31.87424 39.08211 -0.816 0.415291
## wid 0.16331 0.08057 2.027 0.043402 *
## slonce 1271.79725 292.41472 4.349 1.79e-05 ***
## temp -51.13960 29.90291 -1.710 0.088102 .
## rosa 146.95293 29.97451 4.903 1.44e-06 ***
## wilg -38.10145 9.60260 -3.968 8.77e-05 ***
## slonce:temp -40.82650 10.37554 -3.935 0.000100 ***
## rosa:wilg -0.63161 0.16576 -3.810 0.000163 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 675.7 on 356 degrees of freedom
## Multiple R-squared: 0.5785, Adjusted R-squared: 0.569
## F-statistic: 61.07 on 8 and 356 DF, p-value: < 2.2e-165 Wybór najlepszych parametrów dla modelu
Współczynnik AIC im jest niższy, tym lepiej. Według tego kryterium oceniamy, który model jest lepszy.
## df AIC
## regresja 8 5826.511
## regresja2 10 5803.234Wyszło nieco lepiej dla drugiego modelu, ale wtedy jest więcej parametrów modelu (liczba stopni swobody, oznaczana df). Model jest wtedy bardziej obciążony.
5.1 Regresja krokowa z AIC
Algorytm stepAIC korzysta z modelu regresji zrobionego wcześniej. Tworzy kolejne modele, usuwając po kolei zmienne objaśniające.
## Start: AIC=4788.69
## liczba ~ temp + wilg + wiatr + wid + rosa + slonce
##
## Df Sum of Sq RSS AIC
## - wiatr 1 314730 175478073 4787.3
## - wid 1 526819 175690162 4787.8
## <none> 175163342 4788.7
## - slonce 1 3167160 178330502 4793.2
## - temp 1 4211710 179375052 4795.4
## - rosa 1 10947435 186110777 4808.8
## - wilg 1 15743689 190907031 4818.1
##
## Step: AIC=4787.34
## liczba ~ temp + wilg + wid + rosa + slonce
##
## Df Sum of Sq RSS AIC
## - wid 1 542524 176020597 4786.5
## <none> 175478073 4787.3
## - slonce 1 2852991 178331064 4791.2
## - temp 1 4067281 179545354 4793.7
## - rosa 1 10833130 186311203 4807.2
## - wilg 1 15550091 191028164 4816.3
##
## Step: AIC=4786.47
## liczba ~ temp + wilg + rosa + slonce
##
## Df Sum of Sq RSS AIC
## <none> 176020597 4786.5
## - slonce 1 2697051 178717648 4790.0
## - temp 1 4515271 180535868 4793.7
## - rosa 1 11922172 187942769 4808.4
## - wilg 1 19539639 195560236 4822.9##
## Call:
## lm(formula = liczba ~ temp + wilg + rosa + slonce, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2158.96 -398.08 -28.07 533.63 1501.30
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4770.472 720.996 6.617 1.33e-10 ***
## temp -86.972 28.620 -3.039 0.00255 **
## wilg -53.260 8.425 -6.322 7.66e-10 ***
## rosa 149.051 30.185 4.938 1.21e-06 ***
## slonce 321.976 137.091 2.349 0.01938 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 699.2 on 360 degrees of freedom
## Multiple R-squared: 0.5436, Adjusted R-squared: 0.5385
## F-statistic: 107.2 on 4 and 360 DF, p-value: < 2.2e-16## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## liczba ~ temp + wilg + wiatr + wid + rosa + slonce
##
## Final Model:
## liczba ~ temp + wilg + rosa + slonce
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 358 175163342 4788.686
## 2 - wiatr 1 314730.4 359 175478073 4787.341
## 3 - wid 1 542524.1 360 176020597 4786.468Algorytm stepAIC w finalnym modelu uznał parametry wiatr i wid za zbędne.
5.2 Algorytm LEAPS
ALgorytm LEAPS wybiera najlepszy podzestaw zmiennych x dla regresji liniowej. LEAPS sam sobie tworzy model regresji jeszcze raz.
Ustawiamy liczbę podzbiorów każdego rozmiaru do zarejestrowania nbest = 6. Jako statystykę ustawiamy tym razem nie R2, ale adjR2. Ze sposobu wyliczania R2 wynika, że rośnie on zawsze ze wzrostem liczby zmiennych. Statystyka adjusted R2 uwzględnia poprawkę zależną od liczby zmiennych.
## Subset selection object
## Call: regsubsets.formula(liczba ~ ., data = dane, nbest = 6)
## 6 Variables (and intercept)
## Forced in Forced out
## temp FALSE FALSE
## wilg FALSE FALSE
## wiatr FALSE FALSE
## wid FALSE FALSE
## rosa FALSE FALSE
## slonce FALSE FALSE
## 6 subsets of each size up to 6
## Selection Algorithm: exhaustive
## temp wilg wiatr wid rosa slonce
## 1 ( 1 ) "*" " " " " " " " " " "
## 1 ( 2 ) " " " " " " " " " " "*"
## 1 ( 3 ) " " " " " " " " "*" " "
## 1 ( 4 ) " " " " " " "*" " " " "
## 1 ( 5 ) " " "*" " " " " " " " "
## 1 ( 6 ) " " " " "*" " " " " " "
## 2 ( 1 ) " " "*" " " " " "*" " "
## 2 ( 2 ) "*" "*" " " " " " " " "
## 2 ( 3 ) "*" " " " " " " "*" " "
## 2 ( 4 ) "*" " " " " "*" " " " "
## 2 ( 5 ) "*" " " " " " " " " "*"
## 2 ( 6 ) "*" " " "*" " " " " " "
## 3 ( 1 ) "*" "*" " " " " "*" " "
## 3 ( 2 ) " " "*" " " " " "*" "*"
## 3 ( 3 ) " " "*" " " "*" "*" " "
## 3 ( 4 ) " " "*" "*" " " "*" " "
## 3 ( 5 ) "*" "*" " " " " " " "*"
## 3 ( 6 ) "*" "*" " " "*" " " " "
## 4 ( 1 ) "*" "*" " " " " "*" "*"
## 4 ( 2 ) "*" "*" " " "*" "*" " "
## 4 ( 3 ) "*" "*" "*" " " "*" " "
## 4 ( 4 ) " " "*" " " "*" "*" "*"
## 4 ( 5 ) " " "*" "*" " " "*" "*"
## 4 ( 6 ) " " "*" "*" "*" "*" " "
## 5 ( 1 ) "*" "*" " " "*" "*" "*"
## 5 ( 2 ) "*" "*" "*" " " "*" "*"
## 5 ( 3 ) "*" "*" "*" "*" "*" " "
## 5 ( 4 ) " " "*" "*" "*" "*" "*"
## 5 ( 5 ) "*" "*" "*" "*" " " "*"
## 5 ( 6 ) "*" " " "*" "*" "*" "*"
## 6 ( 1 ) "*" "*" "*" "*" "*" "*"## Subset selection object
## Call: regsubsets.formula(liczba ~ ., data = dane, nbest = 6)
## 6 Variables (and intercept)
## Forced in Forced out
## temp FALSE FALSE
## wilg FALSE FALSE
## wiatr FALSE FALSE
## wid FALSE FALSE
## rosa FALSE FALSE
## slonce FALSE FALSE
## 6 subsets of each size up to 6
## Selection Algorithm: exhaustive
Analizując wykres, kierujemy się tym, że jeżeli usuniemy czynnik, który ma dużo „białych” pól, to niewiele straci model. W naszym przypadku możemy usunąć wiatr i wid i adjR2 niewiele się zmniejszy (podobnie zalecał algorytm stepAIC).
5.3 Algorytm Earth
Algorytm Earth też tworzy własny model.Używa własnego typu regresji. Następnie stosujemy funkcję evimp() do oszacowania istotności zmiennych w obiekcie earth.
Algorytm earth uznał za najbardziej istotne dla modelu zmienne (czarna linia nsubsets) temperaturę i wilgotność. Podobnie jak poprzednie algorytmy też uznał widoczność i wiatr(nie uwzględniony na wykresie) za zbędne.
5.4 Algorytm boruta
Algorytm boruta stosuje z kolei algorytm lasu losowego, który uruchamia, aby wybrać najlepszą funkcję.
## [1] "temp" "wilg" "wiatr" "wid" "rosa" "slonce"
Największy wkład do modelu według algorytmu boruta mają parametry: temperatura, wilgotność, punkt rosy, nasłonecznienie. Jest to zgodny obraz z poprzednimi algorytmami.
5.5 Dodatkowo inna metoda - sprawdzenie wspóliniowosci zmiennych objasniajacych
VIF - Variance Inflation Factor
## t H R
## 86.00907 19.41772 115.60011## t H W V R S
## 91.720605 21.783168 1.170080 1.640349 119.697590 2.652569Wspolczynnik VIF dla R - punktu rosy ma bardzo duza wartosc. Z tego wynika, ze punkt rosy nalezy odrzucic, bo jest bardziej zalezny od innych predyktorów niż od odpowiedzi modelu.
Jest to zgodne z zaleznosciami fizycznymi dla punktu rosy.
5.6 Walidacja krzyżowa
Walidacja krzyżowa polega na tworzeniu podmodeli danych i potem następuje łączenie wyników. Najpierw tworzymy skrypt kontrolny, gdzie ustawiamy parametry: metoda-kroswalidacja, liczbę zestawów 10. Potem uruchamiamy trenowanie, określając metrykę według której następuje sprawdzanie dopasowania modelu do rzeczywistości, oraz ustawiamy długość zestawów danych z których tworzony jest model (tuneLength = 30).
Control <- trainControl(method = "cv", number = 10, verboseIter = F)
regresja3 <- caret::train(liczba~.,
data=dane,
method = "lm",
metric = 'Rsquared',
trControl = Control,
na.action = na.omit,
tuneLength = 30)5.6.1 Sprawdzenie predykcji modelu utworzonego przez walidację krzyżową - porównanie z obserwacjami.
liczba_pred3=predict(regresja3)
ggplot(data = dane, aes(x = liczba, y = liczba_pred3)) +
geom_point() +
geom_smooth(method = 'lm', formula = y~x-1) +
labs(title = "Obserwacje kontra model z walidacji krzyżowej",
caption = "Dane dla godziny 18:00 w poszczególnych dniach") +
ylab("liczba model3") +
theme_bw()
Z wykresu widać, że ten model liniowy też jest średniej jakości. Punkty rzeczywiste i prognozowane niezbyt układają się w linię prostą. Ponadto model zaniża wartości wyjściowe – liczbę wypożyczanych rowerów.
5.6.2 Sprawdzenie dopasowania modelu za pomocą funkcji ggof (graphical goodness of fit)
## [,1]
## ME 0.00
## MAE 534.79
## MSE 479899.57
## RMSE 692.75
## NRMSE % 67.30
## PBIAS % 0.00
## RSR 0.67
## rSD 0.74
## NSE 0.55
## mNSE 0.42
## rNSE -Inf
## d 0.84
## md 0.66
## rd -Inf
## cp 0.54
## r 0.74
## R2 0.55
## bR2 0.47
## KGE 0.63
## VE 0.64
Wykres parametrów modelu z pakietu hydroGOF potwierdza słabej jakości model utworzony przez walidację krzyżową. Nie prognozuje on dobrze nagłych dużych skoków wartości wypożyczanych rowerów.
Poszczególne parametry na tym wykresie oznaczają:
- br2 – współczynnik determinacji pomnożony przez nachylenie linii regresji. Model, który systematycznie przewiduje zawyżenie lub zaniżenie przez cały czas, nadal będzie skutkował „dobrym” r2 (bliskim 1), nawet jeśli wszystkie przewidywania były błędne. Współczynnik br2 umożliwia rozliczanie dla rozbieżności w wielkości dwóch sygnałów (przedstawionych przez „b”), jak również ich dynamiki (przedstawiony przez r2)
- cp - współczynnik stałości między predykcją i obserwacjami. Używany do porównywania wydajności modelu w stosunku do prostego modelu z wykorzystaniem wartości obserwowanej z poprzedniego dnia jako prognozy dla dnia bieżącego.
- d - Ta funkcja oblicza indeks zgodności między predykcją i obserwacjami, z uwzględnieniem braków wartości.
- KGE- Efektywność Klinga-Gupty między symulacją a obserwacjami , z uwzględnieniem brakujących wartości. Ta miara dopasowania została opracowana, aby zapewnić diagnostykę rozkładu wydajności Nasha-Sutcliffe’a (a tym samym MSE), który ułatwia analiza względnego znaczenia różnych jego składników (korelacja, odchylenie i zmienność).
- mae - Mean Absolute Error
- md - Ta funkcja oblicza zmodyfikowany indeks zgodności między simulacją i obserwacjami, z uwzględnieniem brakujących wartości
- me - Mean Error
- mNSE - Modified Nash-Sutcliffe efficiency
- mse - Mean Squared Error
- nrmse - Normalized Root Mean Square Error
- NSE - Nash-Sutcliffe Efficiency
- pbias - Percent Bias
- pbiasfdc - współczynnik w hydrologii
- pfactor - współczynnik P to procent obserwacji, które mieszczą się w określonych granicach niepewności.
- rd - Relative Index of Agreement – względny wskaźnik zgodności między symulacją a obserwacjami
- rfactor - współczynnik R reprezentuje średnią szerokość podanych granic niepewności podzieloną przez odchylenie standardowe obserwacji.
- rmse - Root Mean Square Error. RMSE podaje odchylenie standardowe błędu prognozowania modelu. Mniejsza wartość oznacza lepszą wydajność modelu.
- rNSE - Relative Nash-Sutcliffe efficiency - rSD - Ratio of Standard Deviations. Stosunek odchyleń standardowych między symulacją i obserwacjami, z uwzględnieniem brakujących wartości.
- rsr - Ratio of RMSE to the standard deviation of the observations
- ve - współczynnik w hydrologii
5.6.3 Sprawdzenie modelu na podstawie analizy reszt
Najlepiej jest, gdy koniec na powyższym wykresie „podnosi się” jak najpóźniej.
5.6.4 Histogram reszt
Histogram reszt potwierdza to, że model nie jest najlepszy. Występuje pewna ilość reszt o dużej wartości.
6 Utworzenie modelu z podziałem danych na treningowe i testowe.
6.1 Podział danych na treningowe i testowe.
Dzielimy zbiór danych ta zbiór treningowy (uczący) i testowy (walidacyjny) mniej więcej po połowie. Tworzymy model na danych treningowych prz użyciu wszystkich predyktorów.
trening = dane[1:180, ] # zbiór treningowy (uczący)
test = dane[181:365,] # zbior testowy (walidacyjny)
regresja4=lm(liczba ~ ., data = trening)
summary(regresja4)##
## Call:
## lm(formula = liczba ~ ., data = trening)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1982.1 -257.0 42.1 293.5 982.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2845.91065 643.19631 4.425 1.70e-05 ***
## temp -33.23855 22.64165 -1.468 0.143913
## wilg -33.46476 7.26650 -4.605 7.95e-06 ***
## wiatr -0.32401 32.97812 -0.010 0.992172
## wid 0.16656 0.07535 2.210 0.028389 *
## rosa 99.80990 24.23162 4.119 5.89e-05 ***
## slonce 587.94119 171.86584 3.421 0.000779 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 447.3 on 173 degrees of freedom
## Multiple R-squared: 0.7535, Adjusted R-squared: 0.7449
## F-statistic: 88.11 on 6 and 173 DF, p-value: < 2.2e-166.2 Dobór najistotniejszych parametrów nowego modelu
6.2.1 Algorytm LEAPS
## Subset selection object
## Call: regsubsets.formula(liczba ~ ., data = trening, nbest = 5)
## 6 Variables (and intercept)
## Forced in Forced out
## temp FALSE FALSE
## wilg FALSE FALSE
## wiatr FALSE FALSE
## wid FALSE FALSE
## rosa FALSE FALSE
## slonce FALSE FALSE
## 5 subsets of each size up to 6
## Selection Algorithm: exhaustive
## temp wilg wiatr wid rosa slonce
## 1 ( 1 ) "*" " " " " " " " " " "
## 1 ( 2 ) " " " " " " " " " " "*"
## 1 ( 3 ) " " " " " " " " "*" " "
## 1 ( 4 ) " " "*" " " " " " " " "
## 1 ( 5 ) " " " " "*" " " " " " "
## 2 ( 1 ) " " "*" " " " " "*" " "
## 2 ( 2 ) "*" "*" " " " " " " " "
## 2 ( 3 ) "*" " " " " "*" " " " "
## 2 ( 4 ) "*" " " " " " " " " "*"
## 2 ( 5 ) "*" " " " " " " "*" " "
## 3 ( 1 ) " " "*" " " " " "*" "*"
## 3 ( 2 ) " " "*" " " "*" "*" " "
## 3 ( 3 ) "*" "*" " " " " "*" " "
## 3 ( 4 ) " " "*" "*" " " "*" " "
## 3 ( 5 ) "*" " " " " "*" " " "*"
## 4 ( 1 ) " " "*" " " "*" "*" "*"
## 4 ( 2 ) "*" "*" " " " " "*" "*"
## 4 ( 3 ) " " "*" "*" " " "*" "*"
## 4 ( 4 ) "*" "*" " " "*" "*" " "
## 4 ( 5 ) " " "*" "*" "*" "*" " "
## 5 ( 1 ) "*" "*" " " "*" "*" "*"
## 5 ( 2 ) " " "*" "*" "*" "*" "*"
## 5 ( 3 ) "*" "*" "*" " " "*" "*"
## 5 ( 4 ) "*" "*" "*" "*" "*" " "
## 5 ( 5 ) "*" "*" "*" "*" " " "*"
## 6 ( 1 ) "*" "*" "*" "*" "*" "*"## Subset selection object
## Call: regsubsets.formula(liczba ~ ., data = trening, nbest = 5)
## 6 Variables (and intercept)
## Forced in Forced out
## temp FALSE FALSE
## wilg FALSE FALSE
## wiatr FALSE FALSE
## wid FALSE FALSE
## rosa FALSE FALSE
## slonce FALSE FALSE
## 5 subsets of each size up to 6
## Selection Algorithm: exhaustive
6.2.2 “Nowa” metoda doboru parametrów modelu (których predyktorów użyć, a których nie uwzględniać) - algorytm bestGLM (best subset GLM)
library(bestglm)
dane.for.bestglm <- within(dane, {
y <- liczba # liczba do y
liczba <- NULL # Delete bwt
})
res.bestglm <-
bestglm(Xy = dane.for.bestglm,
family = gaussian,
IC = "AIC", # Information criteria for
method = "exhaustive")
res.bestglm$BestModels## temp wilg wiatr wid rosa slonce Criterion
## 1 TRUE TRUE FALSE FALSE TRUE TRUE 4784.468
## 2 TRUE TRUE FALSE TRUE TRUE TRUE 4785.341
## 3 TRUE TRUE TRUE FALSE TRUE TRUE 4785.782
## 4 TRUE TRUE TRUE TRUE TRUE TRUE 4786.686
## 5 TRUE TRUE FALSE FALSE TRUE FALSE 4788.018##
## Call:
## lm(formula = y ~ ., data = data.frame(Xy[, c(bestset[-1], FALSE),
## drop = FALSE], y = y))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2158.96 -398.08 -28.07 533.63 1501.30
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4770.472 720.996 6.617 1.33e-10 ***
## temp -86.972 28.620 -3.039 0.00255 **
## wilg -53.260 8.425 -6.322 7.66e-10 ***
## rosa 149.051 30.185 4.938 1.21e-06 ***
## slonce 321.976 137.091 2.349 0.01938 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 699.2 on 360 degrees of freedom
## Multiple R-squared: 0.5436, Adjusted R-squared: 0.5385
## F-statistic: 107.2 on 4 and 360 DF, p-value: < 2.2e-16Algorytm wyznaczył najlepsze predyktory: temp, wilg, rosa.
6.2.3 Kolejna “nowa metoda” - sprawdzenie wspóliniowosci zmiennych objasniajacych VIF - Variance Inflation Factor
## t H R
## 86.00907 19.41772 115.60011## t H W V R S
## 91.720605 21.783168 1.170080 1.640349 119.697590 2.652569Współczynnik VIF dla R - punktu rosy ma bardzo dużą wartość. Z tego wynika, że punkt rosy należy odrzucic, bo jest bardziej zależny od innych predyktorów niż od odpowiedzi modelu. Jest to zgodne z zależnosciami fizycznymi dla punktu rosy.
7 Model finalny
7.1 Utworzenie modelu finalnego tylko dla najistotniejszych predyktorów
##
## Call:
## lm(formula = liczba ~ temp + wilg, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2106.80 -397.88 -4.99 556.41 1625.76
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1363.408 112.733 12.094 < 2e-16 ***
## temp 59.056 3.159 18.697 < 2e-16 ***
## wilg -15.025 1.999 -7.518 4.43e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 724 on 362 degrees of freedom
## Multiple R-squared: 0.508, Adjusted R-squared: 0.5053
## F-statistic: 186.9 on 2 and 362 DF, p-value: < 2.2e-16
Równanie finalnego modelu regresji wielokrotnej: \[liczba = 1363.408+59.056*temp-15.026*wilg\]
Dopasowanie tego modelu jest najlepsze w porównaniu z poprzednimi modelami. R2 = 0.695
7.2 Diagnostyka modelu finalnego
##
## Shapiro-Wilk normality test
##
## data: reszty
## W = 0.97251, p-value = 2.095e-06Algorytm shapiro-wilk podaje, że p_value < 0.05 czyli nie ma rozkładu normalnego reszt.
##
## Title:
## Jarque - Bera Normalality Test
##
## Test Results:
## STATISTIC:
## X-squared: 14.2266
## P VALUE:
## Asymptotic p Value: 0.0008142
##
## Description:
## Fri Aug 21 19:34:36 2020 by user: RyszardW teście jarquebera wskaźnik X wynosi około 14 czyli jest większy od wartości krytycznej 6. Oznacza to, że reszty nie maja rozkładu normalnego.
Durbin-Watson test, czyli testowanie wartości resztowych pod względem autokorelacji. Występowanie autokorelacji wskazuje na to, ze modelowi brakuje jakiegoś predyktora (zmiennej objaśniającej) lub należy uwzględnić składową szeregu czasowego, np. trend lub wskaźnik sezonowości.
##
## Durbin-Watson test
##
## data: regresja5
## DW = 1.5646, p-value = 1.155e-05
## alternative hypothesis: true autocorrelation is greater than 0ACF - AutoCorrelation Function 
To potwierdza, ze reszty tu są ze sobą skorelowane - wartości na wykresie przekraczają niebieską linię.
Test Breuscha-Pagana jednorodności wariancji
##
## studentized Breusch-Pagan test
##
## data: regresja5
## BP = 37.578, df = 2, p-value = 6.921e-09Mała wartość p oznacza występowanie istotnej zależności pomiędzy wariancją reszt, a liniowymi predyktorami modelu.
Wykresy diagnostyczne dla modelu 
7.3 Predykcja na zbiorze testowym
7.3.1 Rzeczywistość (obserwacje) kontra model (prognoza)
ggplot(data = test, aes(x = liczba, y = liczba_pred5)) +
geom_point() +
geom_smooth(method = 'lm', formula = y~x-1) +
ylab("liczba rowerów (model5)") +
labs(title = "Predykcja finalnego modelu na danych testowych",
subtitle = "Obserwacje kontra model",
caption = "Dane dla godziny 18:00 w poszczególnych dniach") +
theme_bw()
7.3.2 Rzeczywistość kontra model - przebieg w czasie
test0 = dane0[181:365, ] # dane0 z kolumna data
test0 = cbind(test0, liczba_pred5)
typ1 = as.data.frame(rep(1, 185))
colnames(typ1) <- "rodzaj_danych"
typ2 = as.data.frame(rep(2, 185))
colnames(typ2) <- "rodzaj_danych"
typ = rbind(typ1, typ2)
typ = factor(typ$rodzaj_danych)
typ = as.data.frame(typ)
colnames(typ) <- "rodzaj_danych"
rm(typ1, typ2)
dane_typ1 = test0[ ,1:2]
dane_typ2 = test0[ ,c(1,9)] # lub c(1, 10) jezeli jest kolumna nr_wiersza
colnames(dane_typ2) <- c("data", "liczba")
dane_typ = rbind(dane_typ1, dane_typ2)
test00 = cbind(typ, dane_typ)
rm(dane_typ1, dane_typ2)
ggplot(data = test00) +
geom_path(aes(x = data , y = liczba, color = rodzaj_danych)) +
labs(title = "Predykcja finalnego modelu na danych testowych",
subtitle = "Porównanie zmian danych w czasie",
caption = "Dane dla godziny 18:00 w poszczególnych dniach") +
xlab("Data (2018 rok)") +
ylab("Liczba wypożyczonych rowerów") +
scale_color_discrete(name = "Rodzaj danych",
labels = c("rzeczywiste", "predykcja")) +
theme_bw()
7.3.3 Rzeczywistosc kontra model - test funkcją ggof
## [,1]
## ME -60.81
## MAE 727.66
## MSE 775163.59
## RMSE 880.43
## NRMSE % 92.90
## PBIAS % -3.10
## RSR 0.93
## rSD 0.61
## NSE 0.13
## mNSE 0.04
## rNSE -Inf
## d 0.60
## md 0.42
## rd -Inf
## cp 0.44
## r 0.42
## R2 0.17
## bR2 0.15
## KGE 0.30
## VE 0.63
7.3.4 Walidacja założeń modelu
##
## Call:
## lm(formula = liczba ~ temp + wilg, data = dane)
##
## Coefficients:
## (Intercept) temp wilg
## 1363.41 59.06 -15.02
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = regresja5)
##
## Value p-value Decision
## Global Stat 61.60509 1.334e-12 Assumptions NOT satisfied!
## Skewness 9.27731 2.320e-03 Assumptions NOT satisfied!
## Kurtosis 0.08074 7.763e-01 Assumptions acceptable.
## Link Function 9.95288 1.606e-03 Assumptions NOT satisfied!
## Heteroscedasticity 42.29415 7.853e-11 Assumptions NOT satisfied!7.3.5 Wykrywanie najbardziej znaczących obserwacji
## Potentially influential observations of
## lm(formula = liczba ~ temp + wilg, data = dane) :
##
## dfb.1_ dfb.temp dfb.wilg dffit cov.r cook.d hat
## 6 -0.01 -0.01 0.02 0.02 1.03_* 0.00 0.02
## 18 -0.02 -0.03 0.04 0.05 1.03_* 0.00 0.02
## 53 -0.02 -0.02 0.03 0.04 1.03_* 0.00 0.02
## 55 0.02 -0.02 -0.01 0.03 1.03_* 0.00 0.02
## 56 0.02 -0.02 -0.01 0.02 1.03_* 0.00 0.02
## 57 0.03 -0.03 -0.01 0.04 1.03_* 0.00 0.02
## 61 -0.02 -0.03 0.04 0.05 1.03_* 0.00 0.02
## 90 0.02 0.02 -0.03 -0.04 1.03_* 0.00 0.02
## 105 -0.02 -0.01 0.03 0.03 1.03_* 0.00 0.02
## 126 0.03 0.02 -0.05 -0.06 1.03_* 0.00 0.02
## 129 0.01 0.02 -0.03 -0.03 1.03_* 0.00 0.02
## 132 -0.24 -0.03 0.23 -0.27 0.96_* 0.02 0.01
## 144 0.05 0.03 -0.08 -0.09 1.03_* 0.00 0.02
## 200 -0.09 -0.17 0.13 -0.24 0.96_* 0.02 0.01
## 201 -0.02 0.08 0.03 0.15 0.97_* 0.01 0.00
## 234 -0.07 -0.24 0.14 -0.30_* 0.97_* 0.03 0.01
## 241 -0.07 -0.24 0.14 -0.30_* 0.97_* 0.03 0.01
## 292 -0.03 -0.13 0.04 -0.20 0.94_* 0.01 0.00
## 293 0.01 -0.10 -0.02 -0.18 0.95_* 0.01 0.00
## 302 0.00 -0.05 -0.03 -0.14 0.96_* 0.01 0.00
## 304 -0.14 -0.05 0.13 -0.20 0.95_* 0.01 0.01
## 306 -0.05 -0.04 0.03 -0.14 0.96_* 0.01 0.00
## 308 -0.04 -0.10 0.04 -0.18 0.95_* 0.01 0.00
## 313 -0.08 -0.03 0.06 -0.14 0.96_* 0.01 0.00
Tam gdzie jest gwiazdka, to punkt jest nadmiernie znaczący, czyli ta obsewacja ma najwiekszy wpływ na model regresji, czyli zaburza caly model.
Tutaj wyświetlamy tylko podsumowanie, bo influence to cały duży obiekt. W pierwszej kolumnie jest numer wiersza z odstającą obserwacją.
7.3.6 Wnioski
Model dobrze się dopasował do danych treningowych, ale predykcja na danych testowych wychodzi słabo. Jest to spowodowane niezbyt dużym zbiorem danych (365 wierszy) oraz nie spenieniem założeń dla danych do regresji liniowej.
8 Porównanie modeli
## df AIC
## regresja4 8 2716.855
## regresja5 4 5847.688## df BIC
## regresja4 8 2742.398
## regresja5 4 5863.287
Wykres pokazuje wpływ zmiennych na poszczególne modele regresji. Największy wpływ (po wyrazie wolnym) ma temperatura i skorelowane z nią nasłonecznienie.
9 Regresja metodą najmniejszych kwadratów
##
## Call:
## lm(formula = N_test ~ t_test)
##
## Coefficients:
## (Intercept) t_test
## 1253.82 31.02##
## Call:
## lm(formula = N_test ~ t_test)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2023.1 -475.9 159.5 701.3 1554.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1253.818 190.753 6.573 4.98e-10 ***
## t_test 31.018 7.793 3.980 9.92e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 911.7 on 183 degrees of freedom
## Multiple R-squared: 0.07967, Adjusted R-squared: 0.07464
## F-statistic: 15.84 on 1 and 183 DF, p-value: 9.923e-05dane6 = data.frame(t, N)
ggplot(dane6, aes(x = t, y = N)) +
geom_point() +
stat_smooth(method = "lm", col = "blue") +
xlab("temperatura [st.C]") +
ylab("liczba rowerów") +
labs(title = "Regresja metodą najmniejszych kwadratów dla temperatury",
subtitle = "Zbiór danych treningowy",
caption = "Dane dla godziny 18:00 w poszczególnych dniach") +
scale_y_continuous(limits = c(0, 3600)) +
theme_bw()
Przedzial ufności dla współczynników a i b.
## 2.5 % 97.5 %
## (Intercept) 877.46022 1630.17668
## t_test 15.64188 46.39505Górna granica ufności parametru b w tym modelu (temperatura) wynosi 46.3905