Wstęp i dane

Dane opisują ilości rowerzystów przejeżdżających przez różne mosty Nowego Jorku w poszczególnych dniach, a także temperaturę i opady występujące w tychże dniach. Obejmują one okres od kwietnia do października 2017 roku, włącznie.

Praca ta będzie polegała na znalezieniu jak najlepszego modelu regresyjnego, pozwalającego przewidzieć łączną ilość osób jeżdżących na rowerze, uwzględniając także miesiąc, dzień tygodnia, jak i obowiązujące w USA święta. Zostaną użyte modele regresyjnego KNN oraz regresji wielorakiej. Pierwszy model jest wymogiem, drugi zaś został wybrany, ponieważ jest prosty w interpretacji i, jak się okaże w części EDA, zmienne są zależne liniowo. Innym modelem który prawdopodobnie dałby dobre wyniki, mógłby być regresyjny las losowy.

Dane nie były najprostsze do wczytania. Były one w kilku plikach (jeden plik na miesiąc) i dodatkowo w formacie excel worksheet. Podczas parsowania plików, przycinania komórek do właściwej tabeli, jednocześnie dane były wstępnie czyszczone:

Na końcu dane zostały złączone w jedną tabelę. Ostatecznie, dane wyglądają następująco:

## # A tibble: 214 × 8
##    Date       HighT  LowT Precipitation Total Month   Day   Dry
##    <date>     <dbl> <dbl>         <dbl> <dbl> <dbl> <int> <dbl>
##  1 2017-04-01  46    37            0     5397     4     6     1
##  2 2017-04-02  62.1  41            0    13033     4     7     1
##  3 2017-04-03  63    50            0.03 16325     4     1     0
##  4 2017-04-04  51.1  46            1.18  6581     4     2     0
##  5 2017-04-05  63    46            0    17991     4     3     1
##  6 2017-04-06  48.9  41            0.73  4896     4     4     0
##  7 2017-04-07  48    43            0    10341     4     5     0
##  8 2017-04-08  55.9  39.9          0    11610     4     6     1
##  9 2017-04-09  66    45            0    14899     4     7     1
## 10 2017-04-10  73.9  55            0    21295     4     1     1
## # ℹ 204 more rows

EDA

Podsumowanie danych:

##       Date                HighT           LowT       Precipitation   
##  Min.   :2017-04-01   Min.   :46.0   Min.   :37.00   Min.   :0.0000  
##  1st Qu.:2017-05-24   1st Qu.:66.9   1st Qu.:55.23   1st Qu.:0.0000  
##  Median :2017-07-16   Median :75.9   Median :64.00   Median :0.0000  
##  Mean   :2017-07-16   Mean   :74.2   Mean   :62.03   Mean   :0.1318  
##  3rd Qu.:2017-09-07   3rd Qu.:82.0   3rd Qu.:70.00   3rd Qu.:0.0375  
##  Max.   :2017-10-31   Max.   :93.9   Max.   :78.10   Max.   :3.0300  
##      Total           Month             Day         Dry        
##  Min.   : 2374   Min.   : 4.000   Min.   :1   Min.   :0.0000  
##  1st Qu.:15705   1st Qu.: 5.000   1st Qu.:2   1st Qu.:0.0000  
##  Median :19367   Median : 7.000   Median :4   Median :1.0000  
##  Mean   :18628   Mean   : 7.009   Mean   :4   Mean   :0.6121  
##  3rd Qu.:23152   3rd Qu.: 9.000   3rd Qu.:6   3rd Qu.:1.0000  
##  Max.   :26969   Max.   :10.000   Max.   :7   Max.   :1.0000

Jak widać dane są bardzo dobrze przygotowane i nie zawierają błędów. Nie ma żadnych wartości nielogicznych (np. ujemnych), ani za bardzo odstających. Wartości maksymalne w zmiennej Precipitation są dość wysokie, ale jest to normalne, że kilka razy w roku opady są wysokie (szczególnie w okresie wiosennym, czy jesiennym).

Histogramy danych są skośne, lecz również wyglądają na poprawne. Jako, że rozważane są dane dotyczące miesięcy cieplejszych, temperatury wyższe są częstsze. Z histogramu opadów można wywnioskować, że większość dni była sucha, a to wszystko tłumaczy też histogram łącznej liczby rowerzystów. Niskie wartości ‘Total’ będą w dniach które są chłodniejsze i mokre, tylko mocni rowerzyści zdecydują się wtedy na taki sposób podróży. A dlatego, że takich dni jest mało, lewa część histogramu jest bardziej płaska.

Bardzo dobrze na boxplocie opadów widać jak większość obserwacji to dni suche, a prawie wszystkie dni z opadami są wartościami odstającymi.

Analiza korelacji

Zmienna ‘Total’ jest skorelowana ze wszystkimi zmiennymi, najmniej jednak ze zmienną ‘Month’. Skorelowane ze sobą są również wysoka i niska temperatura, jak i opady i zmienna ‘Dry’ (oczywiście).

Sprawdzić można również liniowość zależności zmiennej ‘Total’ od innych, wykorzystując różnicę korelacji Spearman’a i Pearson’a:

cp = cor_matrixp[5,]
cs = cor_matrixs[5,]
print(abs(cs) - abs(cp))
##         HighT          LowT Precipitation         Total         Month 
##         -0.04         -0.03          0.03          0.00          0.00 
##           Day           Dry 
##          0.00          0.00

W żadnym wypadku różnica nie wyszła znacząco wyższa od zera, więc można stwierdzić, że zmienne będą zależne liniowo.

Sprawdzić to można na wykresach rozrzutu (dodatkowo zostały przedstawione wykresy boxplot przedstawiające ilość rowerzystów w zależności od miesiąca i dnia tygodnia):

Zmienna ‘Total’ jest zależna liniowo od najwyższej i najniższej temperatury. Jest też zależna od opadów, z wykresu jednak ciężko stwierdzić, czy jest to zależność liniowa.

Na podstawie wykresu ostatniego, można stwierdzić, że ilość rowerzystów zależy też od tego, czy jest to weekend, czy nie. Zostanie więc wprowadzona zmienna o tym informująca.

data$Weekend = as.numeric(data$Day>=6)

Po wykonaniu modeli, czasami występowały predykcje bardzo mocno odstające od rzeczywistych wartości. Po sprawdzeniu w jakie dni występowały one, okazało się, że są to święta w USA. Została więc dodana kolejna zmienna, informująca o tym czy dzień jest świętem (święta w tym okresie zostały sprawdzone i wprowadzone ręcznie):

data$Holiday = (data$Date == '2017-05-29') | (data$Date == '2017-06-19') | (data$Date == '2017-07-04')|
  (data$Date == '2017-09-04') | (data$Date == '2017-10-09')
data$Holiday = as.numeric(data$Holiday)

Podział danych

Z danych zostały wydzielone 4 obserwacje do sprawdzenia modelu, reszta z nich będzie uczyła model. Obserwacje testowe:

sets = sample(nrow(data), 4)
train = data[-sets,]
test = data[sets,]

test
## # A tibble: 4 × 10
##   Date       HighT  LowT Precipitation Total Month   Day   Dry Weekend Holiday
##   <date>     <dbl> <dbl>         <dbl> <dbl> <dbl> <int> <dbl>   <dbl>   <dbl>
## 1 2017-10-11  75    64.9          0.06 20657    10     3     0       0       0
## 2 2017-06-29  81    68            0    24838     6     4     1       0       0
## 3 2017-06-19  87.1  70            1.35 13087     6     1     0       0       1
## 4 2017-06-30  88    73.9          0    20344     6     5     0       0       0

Regresja wieloraka

Model będzie wykorzystywał wszystkie zmienne stworzone dodatkowo i zmienne pogodowe, oprócz ‘LowT’, gdyż jest ona wysoko skorelowana z ‘HighT’.

model = lm(Total~HighT+Precipitation+Dry+Weekend+Holiday, data=train)
summary(model)
## 
## Call:
## lm(formula = Total ~ HighT + Precipitation + Dry + Weekend + 
##     Holiday, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8109.4 -1658.3   203.5  1903.5  6052.2 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4348.70    1502.12   2.895   0.0042 ** 
## HighT           194.83      20.04   9.722  < 2e-16 ***
## Precipitation -4805.97     577.13  -8.327 1.19e-14 ***
## Dry            3357.57     462.02   7.267 7.65e-12 ***
## Weekend       -4995.85     443.55 -11.263  < 2e-16 ***
## Holiday       -8278.74    1478.04  -5.601 6.83e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2913 on 204 degrees of freedom
## Multiple R-squared:  0.732,  Adjusted R-squared:  0.7254 
## F-statistic: 111.4 on 5 and 204 DF,  p-value: < 2.2e-16

Jak widać wszystkie zmienne są statystycznie znaczące.

Należy również sprawdzić residua modelu:

Jak widać są one losowe, nie widać na wykresie zależności. Histogram jest zbliżony do rozkładu normalnego. Wykresy ACF i PACF wskazują lekką autokorelację dla opóźnienia równego 1 i 9, lecz jest ona bardzo niska i być może wynika z wylosowania danych.

Sprawdzenie modelu na danych testowych:

## [1] "R2: 0.731960644160049 MAPE: 0.176574547486 RMSPE: 0.257261019160756"

R^2 wyszło dość wysokie. Wartości MAPE i RMSPE, oznaczające jak dobrze model radzi sobie z przewidywaniami, wyszły trochę gorsze.

Można rownież wizualnie porównać przewidywania z wartościami rzeczywistymi:

## [1] "Różnice:"
##         1         2         3         4 
##  1984.123  1350.193  6535.071 -1150.075

Na wykresie widać, dlaczego wartości MAPE i RMSPE nie wyszły zbyt dobre. Model potrafi się w jakimś stopniu zbliżyć do oczekiwanych wartości, ale różnice nie są małe.

Regresja KNN

Aby skorzystać z modelu KNN należy zawsze pamiętać o standaryzacji zmiennych. Modele te bazują na odległości euklidesowej, jedna zmienna o wysokich wartościach mogłaby kompletnie zepsuć model.

W tym przypadku zostały znormalizowane zmienne temperatury i opadów. Zmienne binarne pozostały bez zmian. W modelu zostały użyte dokładnie te same zmienne, co w przypadku regresji wielorakiej.

## # A tibble: 210 × 5
##      HighT Precipitation   Dry Weekend Holiday
##      <dbl>         <dbl> <dbl>   <dbl>   <dbl>
##  1 -2.70          -0.328     1       1       0
##  2 -1.15          -0.328     1       1       0
##  3 -1.06          -0.251     0       0       0
##  4 -2.21           2.71      0       0       0
##  5 -1.06          -0.328     1       0       0
##  6 -2.42           1.55      0       0       0
##  7 -2.50          -0.328     0       0       0
##  8 -1.74          -0.328     1       1       0
##  9 -0.773         -0.328     1       1       0
## 10 -0.0133        -0.328     1       0       0
## # ℹ 200 more rows

Parametrem, jaki trzeba ustalić w tym modelu, będzie K, czyli liczba sąsiadów na podstawie której model będzie wyliczał wartość przewidywaną. Pomoże nam w tym wykres ilustrujący RMSE, w modelach o różnej wartości K:

## [1] "Najlepsze K:  5"

Najlepszą wartością K jest 5. Można teraz przejść do stworzenia modelu:

model = knnreg(train_std, train_y, k=najlepszeK)

Należy również sprawdzić residua modelu:

Tutaj również residua na wykresie wyglądają na losowe. Histogram jest podobny do rozkładu normalnego. W tym wypadku nie widać autokorelacji na wykresach ACF i PACF.

Sprawdzenie modelu na danych testowych:

## [1] "R2: 0.835904330256657 MAPE: 0.093535530884217 RMSPE: 0.0956277403274971"

W tym przypadku R^2 wyszło jeszcze lepsze, lecz ogromnej poprawie uległy wartości MAPE i RMSPE. Model, więc dużo lepiej nadaje się do wykonywania predykcji.

Również można wizualnie porównać przewidywania z wartościami rzeczywistymi:

## [1] "Różnice:"
## [1]  2598.000  1784.714  1115.400 -1857.200

Poprawę w predykcji widać “gołym okiem”. Model ten jest w stanie dużo dokładniej przewidywać wartości, szczególnie niską wartość trzecią, która we wcześniejszym modelu bardzo odstawała od rzeczywistości.

Trzecia wartość jest dniem 19 czerwca 2017, jest więc to Dzień Wyzwolenia - narodowe święto w USA. Wniosek: model drugi lepiej uwzględnił więc zmienną “Holiday”.

Można porównać modele wykorzystując kryteria Akaike i Bayesa (AIC oraz BIC):

## [1] "Model 1 AIC: 3440.7746550418 BIC: 3457.51019269538"
## [1] "Model 2 AIC: 3183.55397157418 BIC: 3200.28950922777"

Również z tych statystyk wnika, iż model regresji KNN lepiej nadaje się do wykonywania predykcji.

Podsumowanie i wnioski

Dane były bardzo dobrze zebrane, nie było w nich braków, ani wartości niepoprawnych. Zmienna mówiąca o najniższej temperaturze w ciągu dnia była silnie skorelowana z tą z najwyższą temperaturą, a więc nie była konieczna do uwzględnienia w modelach. Dodatkowo wartości najniższych temperatur prawdopodobnie odczytywane były w nocy, a większość rowerzystów preferuje jazdę w dzień, więc raczej zwracają uwagę na tę wyższą temperaturę.

O ile model regresji wielorakiej nie był zbyt dobry (pierwszą tego przesłanką były autokorelacje residuów), tak model regresji KNN poradził sobie bardzo dobrze. Oba modele poprawnie uwzględniają temperaturę, opady i dni weekendowe, lecz model pierwszy miał problem z uwzględnieniem święta. Dobrze byłoby też uwzględnić w modelach fakt, że całe weekendy, po których (lub przed którymi) występuje święto, mają obniżoną ilość rowerzystów (być może wyjeżdżają wtedy z miasta na długi weekend).

Wprowadzone dodatkowe zmienne świąt, czy wskazanie dni weekendowych poprawiły działanie modelu, lecz możliwe, że brakuje jeszcze paru zmiennych, aby modele mogły lepiej przewidywać ilości rowerzystów (np. wymieniona wcześniej zmienna długiego weekendu - nie tylko dnia święta).

Prawdopodobnie dobrym rozwiązaniem byłoby zastosowanie dwóch osobnych modeli, dla dni zwykłych i długich weekendów. Potrzeba byłoby jednak większej ilości danych dla tego drugiego.