library(ggplot2)
library(tidyverse)
summary(sykler$nrides) Min. 1st Qu. Median Mean 3rd Qu. Max.
22 842 1709 2041 2854 11715 sykler %>%
ggplot() +
geom_histogram(aes(x = nrides))
MET4 - HEKS - H23
Oppgave 1
Gi en introduksjon til besvarelsen med deskriptiv statistikk som er relevant for oppgavene under.
Løsning
Det kan være lurt å starte med å undersøke hovedvariabelen nrides:
Vi ser at det er stor variasjon i antall turer. De minste observasjonene knytter seg nok til oppstarten da det kanskje ikke var så mange sykler tilgjengelig (vi har ikke data på antall sykler tilgjengelig). Fordelingen til nrideshar en relativt tung høyre hale, med flere dager med ekstremt mange turer.
Vi ser så på utviklingen i tid og variasjon gjennom året. Vi starter med å plotte nrides som en tidsrekke og se på autokorrelasjonen:
sykler %>%
ggplot() +
geom_line(aes(x = date, y = nrides))
acf(sykler$nrides)
Det helt tydelig en avhengighetsstruktur i denne tidsrekken, med noen litt rare sesongkomponenter gjennom året: Det kan se ut til at vi har noen topper vår og høst, mens det i sommerferien(?) er en liten dump igjen. I tillegg ser vi at i oppstartsfasen var det betydelig mindre turer, muligens p.g.a. få sykler. I denne perioden er det også et litt utypisk mønster.
Vi kan bruke et boxplott for å se nærmere på mønsteret gjennom året:
boxplot(nrides ~ month, data = sykler)
Her ser vi at det er minst turer i vintermånedene, og vi får også bekreftet at det sykles mindre i Juli.
Vi kan ta en tilsvarende titt på den årlige utviklingen:
boxplot(nrides ~ year, data = sykler)
Her skiller spesielt 2018 seg ut med færre turer.
Vi kan så se litt på korrelasjonen mellom antall turer og værvariablene:
sykler %>%
select(-c(date, isweekend, month, year)) %>%
cor() nrides precip precipIntensity temp snow
nrides 1.0000000 -0.21257048 -0.14035871 0.43556435 -0.23117023
precip -0.2125705 1.00000000 0.56909989 0.03283886 -0.07311786
precipIntensity -0.1403587 0.56909989 1.00000000 0.05295014 -0.06907122
temp 0.4355644 0.03283886 0.05295014 1.00000000 -0.45627821
snow -0.2311702 -0.07311786 -0.06907122 -0.45627821 1.00000000Dette ser i utgangspunktet lovende ut, med tanke på å forklare antall turer med værvariablene. Videre ser vi at spesielt variablene precip og precipIntensity er sterkt korrelerert.
Men er assosiasjonen mellom nrides og værvariablene lineær?
pairs(nrides~., data = sykler %>% select(-c(date, isweekend, month, year)))
Fra den øverste raden her ser vi at med unntak av kanskje temp, så er avhengighetstrukturen ikke lineær. En kan da eksperimentere litt med log() på de ulike variablene og se hvilke varianter som gir tilnærmet en lineær sammenheng. Det vi da må være obs på er at værvariablene kan ta verdien 0. En løsning er å legge til en verdi, f.eks 1, i log-transformasjonen:
pairs(log(nrides) ~ log(precip + 1) + log(precipIntensity + 1) + temp + log(snow + 1), data = sykler)
og vi ser at vi nærmer oss noe. Dette er nyttig informasjon når vi skal lage regresjonsmodellen i neste oppgave! Merk at selv om vi har brukt pairs() i MET4, så er det også helt greit å se på ett og ett spredningsplott.
Til slutt tar vi en visuell sjekk om det er forskjell antall turer på helgedager og ukedager:
boxplot(nrides ~ isweekend, data = sykler)
Det er flere turer i ukedager sammenlignet med helgedager.
Oppgave 2
Sett opp en eller flere regresjonsmodeller for å undersøke hvordan antall sykkelturer henger sammen med de andre variablene. Tolk modellen og fokuser spesielt på værvariablene. Valg som gjøres skal begrunnes. Diagnostiser minst én modell.
Løsning
Fra oppgave 1 har vi sett at det er en slags sesongvariasjon gjennom året og at spesielt 2018 skiller seg ut. Dette kan vi til dels fange opp med å legge til de kategoriske variablene month og year som forklaringsvariabler. Videre så vi at log() transformasjonen var nyttige til å forbedre funksjonsformen mellom nrides og flere av værvariablene. Vi setter også opp en modell der vi ikke gjør noen transformasjoner.
modela <- lm(log(nrides) ~ log(precip + 1) + log(precipIntensity + 1) + temp +
log(snow + 1) + isweekend + year + month, data = sykler)
modelb <- lm(nrides ~ precip + precipIntensity + temp +
snow + isweekend + year + month, data = sykler)
library(stargazer)
stargazer(modela, modelb, type = "text", omit = c("month", "year"))
============================================================
Dependent variable:
----------------------------
log(nrides) nrides
(1) (2)
------------------------------------------------------------
log(precip + 1) -0.135***
(0.014)
log(precipIntensity + 1) -0.041**
(0.020)
precip -28.894***
(2.649)
precipIntensity -15.635***
(4.950)
temp 0.035*** 86.586***
(0.004) (7.797)
log(snow + 1) -0.077***
(0.030)
snow 6.867
(15.593)
isweekend -0.553*** -906.940***
(0.025) (50.469)
Constant 5.678*** -405.338***
(0.063) (120.570)
------------------------------------------------------------
Observations 1,613 1,613
R2 0.721 0.658
Adjusted R2 0.718 0.654
Residual Std. Error (df = 1592) 0.460 914.353
F Statistic (df = 20; 1592) 205.874*** 153.019***
============================================================
Note: *p<0.1; **p<0.05; ***p<0.01Det er ikke mulig å sammenligne forklaringsgraden mellom disse to modellene da den ene beskriver forklart variasjon i nrides, mens den andre beskriver forklart variasjon i log(nrides). Diagnostikken som følger er mer tilfredstillende for modellen hvor vi gjort transformasjoner. Vi fokuserer derfor på tolkningen av denne.
For helgevariabelen isweekend har vi en log-lin relasjon og vi ser at verdien av effekten er så stor at en direkte prosenttolkning vil være dårlig. Helg henger sammen at antall turer blir multiplisert med exp(-0.553) altså en (exp(-0.553) - 1)100% = -42.4% reduksjon, gitt alt annet likt. Det er betydelig flere turer i ukedagene.
Samtlige værvariabler ser ut til å ha en signifikant effekt på sykkelturer. En grads økning i temperatur henger sammen med en 0.034*100% = 3.4% (log-lin) økning i i antall turer, gitt alt annet likt. Siden vi har log(x + 1) på de andre værvariablene må vi se for oss en ny forklarings variabel \(x^* = x + 1\) og endringer i denne når vi tolker. Tolkningen er da analog til log-log tolkningen: 1% økning i precip + 1 henger sammen med en -0.135% reduksjon i antall sykkelturer. De tilsvarende reduksjonene for 1% økning i precipIntensity + 1 og snow + 1 er hhv. -0.041% og -0.077%. Med andre ord: syklister ser ut til å foretrekke varmt vær og minst mulig nedbør og snødekke.
Her viser vi ikke månedseffektene med stargazer (omit = c("month", "year")), men de reflekterer i stor grad boksplottetene vi så i oppgave 1 for begge modellene.
Vi diagnostiserer så residualene for modela og begynner med et residualplott:
residualer <- data.frame(residualer = modela$residuals,
predikert = modela$fitted.values)
ggplot(residualer, aes(x = predikert, y = residualer)) +
geom_point()
Tendenser til heteroskedastisitet med litt større spredning rundt 0 for små prediksjoner sammenlignet med store, men ingen mønster. Alt i alt ikke så verst, iallefall dersom man sammenligner med en modell uten noen log() transformasjoner.
Hva med så med normalantagelsen?
ggplot(residualer, aes(x = residualer)) +
geom_histogram()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(residualer, aes(sample = residualer)) +
stat_qq() +
stat_qq_line()
Her ser vi et betydelig avvik i halene så vi skal være forsiktige med å stole på dekningsgraden til prediksjons- og konfidensintervaller.
Har vi kvittet oss med tidsavhengigheten?
acf(modela$residuals)
Ennå struktur igjen, men i forhold til autokorrelasjonen til nrides i oppgave 1, er det mindre autokorrelasjon og vi har iallefall delvis fanget opp sesongkomponentene siden vi ikke lenger ser disse periodiske svingningene i like stor grad. I praksis betyr det at det er rom for forbedring dersom vi vil bruke modellen til prediksjon frem i tid. Mer om det i oppgave 4!
Vi kan også ta en titt på hvordan de tilpassede verdiene stemmer med de faktiske observasjonene i et tidsrekkeplott. Da må vi huske at prediksjonene er på log-skala så vi bruke exp() for transformere tilbake til originalskalaen:
sykler %>%
ggplot() +
geom_line(aes(x = date, y = nrides)) +
geom_line(aes(x = date, y = exp(modela$fitted)), color = "blue")
Se forøvrig INFO kommentar i oppgave 4b).
Oppgave 3
Undersøk om været har forskjellig effekt på antall turer i helger sammenlignet med ukedager.
Løsning
Fra boxplottet i oppgave 1 og regresjonsmodellen i oppgave 2 så vi at det var betydelig færre turer på helgedager sammenlignet med ukedager. Men dette kan skyldes andre ting enn at syklister reagerer forskjellig på været på de to dagtypene. Det kan f.eks. være p.g.a. et større behov for sykler i ukedager.
Man kunne tenke seg at man delte datasettet i to: ett datasett for ukedager og ett for helg og så sammenlignet regresjonsmodellen vi brukte i oppgave 2 for de to datasettene. Men dersom vi formelt skal teste om været har forskjellig effekt på turer i helger og ukedager må vi bruke interaksjonsledd i èn enkelt regresjon:
modelinter <- lm(log(nrides) ~ + isweekend + year + month +
log(precip + 1) + log(precipIntensity + 1) + temp + log(snow + 1)
+ isweekend*log(precip + 1)
+ isweekend*log(precipIntensity + 1) +
+ isweekend*temp +
+ isweekend*log(snow + 1),
data = sykler)
stargazer(modelinter, type = "text", omit = c("month", "year"))
==================================================================
Dependent variable:
---------------------------
log(nrides)
------------------------------------------------------------------
isweekend -0.494***
(0.064)
log(precip + 1) -0.118***
(0.016)
log(precipIntensity + 1) -0.033
(0.023)
temp 0.033***
(0.004)
log(snow + 1) -0.062*
(0.034)
isweekendTRUE:log(precip + 1) -0.054*
(0.030)
isweekendTRUE:log(precipIntensity + 1) -0.037
(0.043)
isweekendTRUE:temp 0.008
(0.005)
isweekendTRUE:log(snow + 1) -0.040
(0.054)
Constant 5.663***
(0.066)
------------------------------------------------------------------
Observations 1,613
R2 0.724
Adjusted R2 0.720
Residual Std. Error 0.458 (df = 1588)
F Statistic 173.972*** (df = 24; 1588)
==================================================================
Note: *p<0.1; **p<0.05; ***p<0.01Retningen på effekten (fortegnet) av alle interaksjonsledd er i henhold til at syklister foretrekker det enda varmere og vil ha mindre nedbør og snødekke sammenlignet med syklister i ukedager. De tilhørende standardavvikene er relativt store og det er f.eks. kun effekten av interaksjonen mellom helgedag og nedbør (log(precip + 1)) som er signifikant forskjellig fra 0 med et 10% signifikansnivå. En tolkning blir da at en 1% økning i precip + 1 er assosiert med en -0.118% reduksjon i antall turer på ukedager, mens den tilsvarende reduksjonen på helgedager er (-0.118 - 0.054)% = -0.172 %.
Oppgave 4a
Bruk modellen(e) fra oppgave 2 til å gi en generell (og kort) anbefaling om når Bergen Bysykkel bør vedlikeholde sine sykler.
Løsning
Ut fra modellen i oppgave 2 bør Bergen Bysykkel utføre vedlikeholdet en helg i desember eller januar. Ideelt sett er det meldt om kalde temperaturer og nedbør og det kan gjerne ligge litt snødekke.
OBS: Det er meningsløst å si at de bør utføre vedlikeholdet i 2018, med mindre du er skråsikker på at Bergen Bysykkel er i besittelse av en tidsmaskin.
Oppgave 4b
Bruk modellen(e) fra oppgave 2 til å gi en anbefaling om hvilken av de neste 5 dagene Bergen Bysykkel bør utføre vedlikehold.
Løsning
Vi bruker regresjonsmodellen og dataene for de neste 5 dagene til å predikere antall turer (med prediksjonsintervall). Vi må huske at vi da får prediksjoner for log(nrides). For å få prediksjonene (og prediksjonsintervallene) på originalskala kan vi bruke exp() funksjonen:
prediksjoner <- predict(modela,
newdata = neste5,
interval = "prediction",
level = 0.95)
exp(prediksjoner) fit lwr upr
1 876.6275 353.4955 2173.934
2 798.8366 322.0216 1981.668
3 432.9979 174.3798 1075.166
4 441.2860 177.6586 1096.110
5 779.2322 314.1623 1932.768Ut fra denne prediksjonen ser det ut til at dag 3 (lørdag 03.12.2022) er en ideell dag for vedlikehold. Merk at på kort sikt, hvis ikke det er meldt om store væromslag så er det “helgeeffekten” som er viktigst.
INFO: Prediksjonsintervallet over vil være riktig, men selve punktprediksjonen av nrides er generelt underestimert hvis vi gjør det på denne måten. Modellen impliserer at nrides er log-normalfordelt og det kan vises at en bedre prediksjon er gitt ved \(\exp(\textit{prediksjoner\$fit} + \sigma^2/2)\) der \(\sigma = 0.460\) (Residual std.error fra utskriften).
Vi kan se hvordan dette ser ut i en figur der vi tar med antall turer de siste 20 dagene sammen med prediksjonsintervallene for de fem neste dagene. Figuren under er ikke ventet i besvarelsen:
tpred <- as_tibble(exp(prediksjoner))
tpred$date <- neste5$date
tail(sykler, 40) %>%
ggplot() +
geom_line(aes(x = date, y = nrides)) +
geom_line(data = tpred, aes(x = date, y = fit), color = "blue") +
geom_ribbon(data = tpred, aes(x = date, ymin = lwr, ymax = upr), fill = "blue", alpha = 0.2) 
Oppgave 4c
Anta at du kun har tilgjengelig variabelen i datasettet . Finn den tidsrekkemodellen som best predikerer antall sykkelturer, og bruk denne modellen til å gi en anbefaling om hvilken av de neste 5 dagene Bergen Bysykkel bør utføre vedlikehold.
Løsning
Skal vi sammenligne hvor godt noen tidsrekkemodeller klarer å predikere må vi først lage et trenings og et testsett:
trening <- head(sykler$nrides, length(sykler$nrides) - 10)
test <- tail(sykler$nrides, 10)Det er mulig å gjøre en log-transformasjon av tidsrekken først, men da må vi i såfall forholde oss til dette og transformere prediksjonene tilbake til originalskala til slutt. Vi velger for enkelthetsskyld å jobbe på originalskalaen. Man kan også vurdere å droppe de første observasjonene da de kanskje ikke er representative for resten av tidsrekken.
Vi forsøker å predikere testsettet med en arima modell og eksponensiell glatting:
library(forecast)
arimamod <- auto.arima(trening)
arimapred <- forecast(arimamod, h = 10)
expmod <- HoltWinters(trening, gamma = FALSE)
exppred <- forecast(expmod, h = 10)
accuracy(arimapred, test) ME RMSE MAE MPE MAPE MASE
Training set 1.010313 746.7697 511.5097 -12.7818893 33.99807 0.9193044
Test set 145.511125 342.5443 275.0605 -0.5357157 26.25911 0.4943490
ACF1
Training set -0.05459547
Test set NAaccuracy(exppred, test) ME RMSE MAE MPE MAPE MASE ACF1
Training set -26.1311 830.5513 567.3608 -15.69980 37.63904 1.019682 0.0657854
Test set 626.4911 743.9689 701.5265 33.34008 50.01463 1.260810 NAVi ser at den estimerte arimamodellen er i en klasse for seg selv. Ved hjelp av summary(arimamod) ser vi at vi har estimert en arima(2,1,2) modell.
Skal vi bruke denne modellen til å predikere de neste fem dagene er det viktig at vi reestimerer modellen på hele datasettet først. Da slipper vi litt styr for å få predikert de riktige dagene (vi kan bare sette h = 5 i forecast() funksjonen), og vi bruker all tilgjengelig informasjon til estimering av modellen. Det siste er spesielt viktig hvis testdatasettet er stort.
tsmodel <- auto.arima(sykler$nrides)
tsprediksjoner <- forecast(tsmodel, h = 5)
tsprediksjoner Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1614 1550.047 593.84182 2506.252 87.65746 3012.436
1615 1330.543 129.99423 2531.092 -505.53784 3166.624
1616 1256.517 -17.55742 2530.592 -692.01185 3205.047
1617 1327.262 28.97331 2625.552 -658.29945 3312.824
1618 1459.547 135.48942 2783.605 -565.42452 3484.519Også her predikerer vi at dag 3 er den dagen med færrest sykkelturer og følgelig den beste dagen å vedlikeholde på.
Vi kan også lett lage en figur av denne prediksjonen:
autoplot(tsprediksjoner, include = 40)
INFO: Dekomponering ved hjelp av stl() er ikke rett frem fordi vi må konvertere tidsrekken om til et såkalt ts() objekt, som kan være litt vrient. Det fører derfor ikke til noe trekk i besvarelsen dersom dette ikke har blitt prøvd ut. Men for de som er interessert kunne dette ha blitt gjort som følger:
library(forecast)
trening_ts <- ts(trening, start = c(2018, 181), frequency = 365)
stlmod <- stl(trening_ts, s.window = "periodic", )
stlforecast <- forecast(stlmod, h = 10)
accuracy(stlforecast, test) ME RMSE MAE MPE MAPE MASE
Training set 1.475952 734.8036 534.8124 -10.52611 40.21459 0.9611849
Test set 785.299798 1043.2606 934.2778 38.19572 66.18934 1.6791190
ACF1
Training set 0.1031973
Test set NADette er også en ganske god modell som kunne blitt forbedret med å tilpasse en arimamodell til residualtidsrekken.
Oppgave 4d
Vis at forklaringsgraden til modellen du brukte i oppgave 4a/b kan forbedres ved å ta hensyn til tidsavhengigheten
Løsning
Fra autokorrelasjonensplottet i oppgave 2 så vi at de inkluderte forklaringsvariablene ikke helt klarte å redegjøre for tidsavhengigheten i antall turer. Generelt hvis vi legger til \(Y_{t-1}\) som en forklaringsvariabel i en regresjonsmodell så får vi en slags AR(1) modell med forklaringsvariabler.
Vi begynner med å legge til lag-1 av nrides til datasettet sykler:
library(dplyr)
sykler <- sykler %>%
mutate(lag1_nrides = lag(nrides, 1))Siden vi modellerer log(nrides) i modellen fra oppgave 2, er det fornuftig å bruke log(lag1_nrides) som den ekstra forklaringsvariabelen i modellen:
modelny <-lm(log(nrides) ~ log(lag1_nrides) + log(precip + 1) + log(precipIntensity + 1) + temp +
log(snow + 1) + isweekend + year + month, data = sykler)
stargazer(modelny, type = "text", omit = c("month", "year"))
====================================================
Dependent variable:
---------------------------
log(nrides)
----------------------------------------------------
log(lag1_nrides) 0.591***
(0.015)
log(precip + 1) -0.127***
(0.010)
log(precipIntensity + 1) 0.003
(0.014)
temp 0.015***
(0.003)
log(snow + 1) -0.036*
(0.021)
isweekend -0.484***
(0.018)
Constant 2.526***
(0.091)
----------------------------------------------------
Observations 1,612
R2 0.861
Adjusted R2 0.859
Residual Std. Error 0.324 (df = 1590)
F Statistic 469.093*** (df = 21; 1590)
====================================================
Note: *p<0.1; **p<0.05; ***p<0.01Vi ser at dette øker forklaringsgraden betraktelig.
Vi kan ta en titt på autokorrelasjonen til residualene:
acf(modelny$residuals)
og vi ser også her at vi har fått en kraftig forbedring. Toppene som dukker opp med periode 7 dager kan tenkes er at det er spesielt lite aktivitet på søndager.
Ekstra: Vi har ikke sett på hvordan vi kan evaluere prediksjonsevnen på et testsett for slike modeller, men denne modellen er ganske mye bedre enn både tidsrekkemodellen og modellen fra oppgave 2a. Det finnes også et egent argument i auto.arima() kalt xreg som lar deg legge til forklaringsvariabler direkte i arimamodellen. På denne måten kunne vi ha spisset modellen enda mer.