Jak zdefiniować spóźnienie?

Nic trudnego, zdawałoby się. Jeżeli tramwaj miał przyjechać o 8:02, a przyjechał o 8:06, to zapisujemy 4 minuty spóźnienia. Problem pojawia się jeżeli tramwaj nie przyjechał. Jeżeli po 10 minutach przyjechał następny, to być może powinniśmy zapisać 10 minut spóźnienia. Co jednak zrobić jeśli ostatecznie przyjedzie on na przystanek, ale po 12 minutach. Taka sytuacja może się zdarzyć jeśli spóźnienie wynika z nadzwyczajnej zmiany trasy. Może więc 12 minut? Czekających pasażerów, można by odpowiedzieć, niespecjalnie interesuje jednak co się stało już po odjeździe kolejnego tramwaju. Ale z drugiej strony, jeżeli taka sytuacja zdarzy się w miejscu, gdzie wszystkie linie jadące w danym kierunku mają tą samą trasę (np. niedaleko pętli), to ci pasażerowie mogli wsiąść w tramwaj innej linii i czekać mniej niż 10 minut, więc to też nie jest najlepsze kryterium.

Z pewną pomocą przychodzi TTSS, w którym kursy anulowane znikają (choć nie zawsze) z tablic oczekiwanych przyjazdów. Na przykład w środę 25.07 odbywał się remont przy tunelu pod Dworcem Głównym. Jadące tam linie miały zmienione trasy i informacje na temat ich odjazdów pojawiały się tylko na tych ich częściach, które zmieniane nie były. Żeby zachować spójność z systemem z którego korzystam przyjąłem następującą procedurę określania spóźnienia:

1. Jeżeli tramwaj przestał być śledzony przez system, to kurs jest traktowany tak jakby w ogóle nie był planowany.
2. Jeżeli tramwaj zostanie odnotowany na przystanku, to spóźnienie jest liczone jest liczone jako różnica między rzeczywistym, a planowanym czasem przyjazdu. Przy czym jako rzeczywisty czas przyjazdu będę rozumiał albo czas zapytania w którym system pokazał, że tramwaj zatrzymuje się na przystanku, albo czas ostatniego zapytania, gdy przyjazd był jeszcze oczekiwany. Na informacjach o tym, że tramwaj się zatrzymuje nie można do końca polegać.
3. Jeżeli tak wyliczone spóźnienie jest duże, to sprawdzam w komunikatach publikowanych przez MPK, czy da się je powiązać z jakąś awarią. Jeżeli okaże się, że takie spóźnienie wynika ze zbyt późnego „anulowania kursu” (co miało miejsce na przykład przy tymczasowym skróceniu linii 6 i 24 do Prokocimia), to usuwam takie obserwacje ręcznie.

Czas badania

Dane zbierałem od wtorku 24.07 do wtorku 31.07, z pominięciem weekendu. Był to pierwszy tydzień obwiązywania poprawionego rozkładu jazdy, po zakończeniu remontu na Bieńczyckiej.

W weekend odbyła się naprawa szyn w okolicach skrzyżowania ulic Nowosądeckiej i Spółdzielców. Po naprawie średnie opóźnienia obniżyły się, a istotne różnice w rozkładach spóźnień występowały na przystankach Piaski Nowe i Nowosądecka, dla tramwajów jadących do Kurdwanowa, oraz od Witosa aż do Podgórze SKA, dla tych jadących w drugą stronę. Przynajmniej część z tych różnic trudno powiązać bezpośrednio z remontem, bo pojawiają się na kilka przystanków przed naprawianym odcinkiem. Tak czy owak, wyniki dla okolic pętli Kurdwanów P+R oraz dla początkowych odcinków linii 6 i 24 kursujących z tej pętli, mogą być z powodu remontu niemiarodajne.

repair_test <- function(x){
print(x)

#select departures to Kurdwanów
to <- report %>% filter(stopName == x & direction == 'Kurdwanów P+R') %>% select(delay, plannedTime) %>% mutate(repair = day(plannedTime)>=30, delay = as.factor(delay)) %>% select(-plannedTime)

#determine which levels of delays should be merged to meet chi square test assumptions
drop <- to %>% group_by(repair) %>% count(delay) %>% filter(n >= 5) %>% arrange(n) %>% summarize(threshold = first(delay)) %>% `$`(threshold) %>% as.character() %>% as.numeric() %>% min()

#merge levels
to$delay <- fct_lump(to$delay, n = drop, other_level = 'bigger')

#test
table(to) %>% chisq.test() %>% print()

from <- report %>% filter(number %in% c(6, 24) & stopName == x & direction != 'Kurdwanów P+R') %>% select(delay, plannedTime) %>% mutate(repair = day(plannedTime)>=30, delay = as.factor(delay)) %>% select(-plannedTime)

drop <- from %>% group_by(repair) %>% count(delay) %>% filter(n >= 5) %>% arrange(n) %>% summarize(threshold = first(delay)) %>% `$`(threshold) %>% as.character() %>% as.numeric() %>% min()

from$delay <- fct_lump(from$delay, n = drop, other_level = 'bigger')

table(from) %>% chisq.test() %>% print()
}

purrr::walk(c("Witosa", "Nowosądecka", "Piaski Nowe", "Dauna", "Bieżanowska", "Kabel", "Dworcowa", "Cmentarz Podgórski", "Podgórze SKA", "Limanowskiego"), repair_test)

## [1] "Witosa"
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 12.435, df = 6, p-value = 0.05295
## 
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 6.9135, df = 2, p-value = 0.03153
## 
## [1] "Nowosądecka"
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 12.803, df = 5, p-value = 0.02529
## 
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 10.589, df = 3, p-value = 0.01417
## 
## [1] "Piaski Nowe"
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 12.371, df = 6, p-value = 0.05419
## 
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 9.1451, df = 3, p-value = 0.02742
## 
## [1] "Dauna"
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 9.2774, df = 6, p-value = 0.1586
## 
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 10.592, df = 3, p-value = 0.01415
## 
## [1] "Bieżanowska"
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 10.293, df = 5, p-value = 0.06734
## 
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 15.393, df = 3, p-value = 0.00151
## 
## [1] "Kabel"
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 11.917, df = 7, p-value = 0.1033
## 
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 8.8538, df = 4, p-value = 0.06486
## 
## [1] "Dworcowa"
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 13.434, df = 5, p-value = 0.01963
## 
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 14.169, df = 5, p-value = 0.01457
## 
## [1] "Cmentarz Podgórski"
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 11.296, df = 8, p-value = 0.1855
## 
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 17.803, df = 4, p-value = 0.001349
## 
## [1] "Podgórze SKA"
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 11.187, df = 5, p-value = 0.0478
## 
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 12.338, df = 5, p-value = 0.03043
## 
## [1] "Limanowskiego"
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 0.84645, df = 3, p-value = 0.8383
## 
## 
##  Pearson's Chi-squared test
## 
## data:  .
## X-squared = 2.0559, df = 4, p-value = 0.7255

Poza tym, wyników analizy nie należy stosować poza wakacyjny rozkład jazdy.

Czyszczenie i braki danych

Z pobranych danych trzeba było wyrzucić te błędne i niemiarodajne:

Rozkład spóźnień

Oto, jak przedstawia się rozkład spóźnień:

Pesymista powiedziałby, że ponad połowa tramwajów przyjeżdża opóźniona. Optymista, że zdecydowaną większość stanowią małe spóźnienia, do 2 minut. Tak czy owak średnie spóźnienie w próbie wyniosło minutę i 0.3 sekundy, a odchylenie standardowe minutę i 23 sekundy (1.39 minuty). Ten rozkład nie daje się łatwo utożsamić z żadnym z rozkładów dobrze znanych z podręczników i to niezależnie czy potraktujemy spóźnienia jako zmienną dyskretną, czy jako ciągłą lecz obarczoną błędami zaokrągleń.

Spóźnienia w czasie

Zobaczmy jak układa się średnia opóźnienia w zależności od godziny:

Prawdę mówiąc spodziewałem się raczej, że będzie wyraźnie widać poranny i popołudniowy szczyt komunikacyjny. Tymczasem rano przyrost spóźnień nie jest tak dynamiczny jak po południu i kumuluje się później niż się myślałem.

Krzywa została dopasowana przy użyciu uogólnionego modelu liniowego, co jest metodologicznie śliskie, zważywszy na problemy z identyfikacją rozkładu. Lepsze byłyby metody nieparametryczne. Próbowałem tutaj użyć metody loess, ale brakło mi pamięci (loess wymaga O(n^2), pamięci, obserwacji w zbiorze jest ok. 400 tysięcy). Niżej będę korzystał z regresji kwantylowej.

Rankingi spóźnień

Przystanki z największym średnim opóźnieniem:

Wysokie miejsca zajmują jak widać przysatnki, na których wpływ na średnie spóźnienie miał remont torowiska 28 i 29 lipca. Gdyby ranking zrobić tylko dla danych z dni po remoncie przystanki Kabel i Nowosądecka wciąż byłyby w pierwszej piątce, a zamiast Dauna i Piasów Nowych pojawiłyby się Zabłocie i Muzeum Lotnictwa.

Przystanki z najmniejszym średnim opóźnieniem (z pominięciem pętli):

Co ciekawe, przystani z obydwu rankingów znajdują się niedaleko pętli. Można się domyślać (potwierdzą to mapy poniżej), że tramwaje które dojeżdżają do pętli mają spore opóźnienia, a te które dopiero co wyjechały - małe. W rankingach znalazły się te przystanki w których jedna strona tego układu przeważa w średniej.

Szczególnym przypadkiem jest przystanek Łagiewniki ZUS. Jest on położony niedaleko pętli w Borku Fałęckim i obsługuje tylko jeden tor: do Borku właśnie. Stąd nie dziwi jego pozycja na czele niechlubnego rankingu.

Linie z największym średnim opóźnieniem:

Linie z najmniejszym średnim opóźnieniem:

Z takiego prostego zestawienia niewiele jeszcze wynika. Przyjrzyjmy się zatem wszystkim średnim spóźnieniom na wszystkich liniach.

W rankingu najmniej spóźniających się linii wysokie miejsca zajęły krótkie linie 2 i 16. Zobaczmy, czy istnieje związek między długością linii, a średnim spóźnieniem.

O ile najkrótsze linie mają niewielkie średnie opóźnienia, to nie widać ogólnej, silnej zależności.

Wnioski z map

• Mapy potwierdzają postawioną wyżej hipotezę o spóźnieniach w okolicach pętli: na samych pętlach spóźnienie jest małe i tramwaje, które dopiero co wyjechały również mają niewielkie spóźnienia. Z kolei tramwaje kończące trasę przyjeżdżają z większym opóźnieniem. W zależności od skali dysproporcji średnie opóźnienia w okolicach pętli należą do jednych z niższych na mapie (Wzgórza Krzesławickie, Os. Piastów, Mistrzejowice, Cmentarz Rakowicki, Cichy Kącik, Czerwone Maki, Nowy Bieżanów), przeciętnych (Kopiec Wandy, Krowodrza Górka, Bronowice Małe, Salwator), lub wyższych (Walcownia, Borek Fałęcki, Kurdwanów, Mały Płaszów). Część tej zależności widzieliśmy w rankingach przystanków.

• Sąsiadujące przystanki mają podobne średnie spóźnienia, co świadczy o sporej korelacji pomiędzy nimi. Przyjrzyjmy się korelacji bliżej.

  Ten wykres powstał następująco: wziąłem wszystkie fragmenty tras o długości 6 przystanków i wyliczyłem średnie spóźnienia na wszystkich przystanków na ścieżce. Następnie wyznaczyłem korelacje między średnimi spóźnieniami na pierwszym i drugim przystanku ścieżki, pierwszym i trzecim itd. Przy liczeniu korelacji każdemu fragmentowi nadawałem wagę równą liczbie zarejestrowanych na niej przejazdów.

• Tramwaje nadrabiają spóźnienia! Na mapce można zauważyć, że po najmocniej czerwonych odcinkach, oznaczających wzrost spóźnienia, często występują mocno niebieskie, oznaczające spadek. Zobaczmy tę zależność na wykresie.

  Jak czytać ten wykres? Przypuśćmy, że przy Filharmonii tramwaj nie miał opóźnienia. Na następnym przystanku, Jubilat, dojechał jednak z dwuminutowym opóźnieniem. Potem udało mu się nadrobić część spóźnienia i na Komorowskiego przyjechał spóźniony minutę. Dla przystanku Jubilat przyrost opóźnienia przed przystankiem wyniósł 2 minuty (od 0 do 2), a po przystanku -1 minutę (od 2 do 1). Taka historia zostałaby w powyższym układzie współrzędnych zaznaczona w punkcie (2, -1). Wykres powstał przez wzięcie średnich z takich sekwencji dla każdego przystanku. Jeżeli na jakimś przystanku dochodziło do rozwidlenia, to wszystkie ścieżki były liczone osobno.

  W prawej części wykresu średnia spada poniżej 0, czyli po odcinkach o dużym przyroście spóźnienia następują te z ujemnym przyrostem (spadkiem). Niestety rzadko udaje się nadrobić całe spóźnienie: większość punktów znajduje się powyżej przekątnej (przekątna odpowiada sytuacjom w których tramwaj nadrabia tyle ile wcześniej stracił)

• Oprócz fragmentów z dużymi wahaniami warto zwrócić uwagę na fragmenty, gdzie tramwaje przez dłuższy odcinek łapią nieduże spóźnienia, ale mają ich gdzie nadrobić. Najwięcej takich miejsc leży w obszarze pomiędzy przystankami Stradom, Św. Gertrudy, Dworzec Główny, Rondo Mogilskie, Rondo Grzegórzeckie. Choć niektóre przyrosty są bardzo małe, to w całym tym obszarze tylko na dwóch odcinkach opóźnienie spada. Podobne odcinki to: wyjazd z Borku Fałęckiego, od Sanktuarium Bożego Miłosierdzia do Smolki, przejazd Aleją Pokoju od Ronda Czyżyńskiego do Grzegórzeckiego, połączenie Dworca Towarowego ze Starym Kleparzem, czy odcinek od Placu Wolnica do Dworca Płaszów.

Regresja kwantylowa

Dopasujmy teraz modele regresji kwantylowej dla obydwu kierunków, uwzględniające także planowany czas odjazdu i interakcję między tymi zmiennymi. Ponieważ szczegółowe potraktowanie przystanków dobrze się sprawdziło wymusimy szczegółowość przez duży wymiar bazy przestrzeni z której będą dobierane splajny.

df52b <- dplyr::filter(report, number == 52, direction == 'Czerwone Maki P+R')
df52b$delay <- jitter(as.numeric(df52b$delay), amount = 0.5)
df52b$delay[df52b$delay < 0] <- 0

qgam_model_A <- mqgam(delay ~ ti(seq_num, k = 30) + ti(minutefs, k = 8) + ti(seq_num, minutefs), data = df52, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_A, fun = summary)

qgam_model_B <- mqgam(delay ~ ti(seq_num, k = 30) + ti(minutefs, k = 8) + ti(seq_num, minutefs), data = df52b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_B, fun = summary)

Nie wklejam wyników podsumowań modelu, żeby dodatkowo nie wydłużać wpisu, ale można sprawdzić, że wszystkie współczynniki wyszły statystycznie istotne.

Przyjrzyjmy się spóźnieniom w godzinach, w których są one największe.

W obydwu kierunkach początkowe spóźnienia są bardzo małe, z kilkoma odstającymi przystankami. Potem następuje rozłożony na kilka odcinków wzrost opóźnień i ich rozpiętości. Co ciekawe, w obydwu kierunkach odbywa się to na podobnym odcinku - między Stradomiem, a Cystersów/Lubicz - który wcześniej wskazaliśmy jako jedno z miejsc gdzie występują małe, ale następujące bez przerwy przyrosty spóźnień.

W kilku miejscach na trasie następują krótkie skoki spóźnień. Dobrze byłoby sprawdzić ich źródła w terenie. Mogą to być np. niekorzystne usytuowanie przystanku względem sygnalizacji świetlnej (podwójne czekanie: najpierw na sygnalizacji, potem na przystanku; dotyczy przystanków Chmieleniec, Norymberska, Kleeberga), trudności ze sprawnym przejazdem przez węzeł komunikacyjny (Centrum Kongresowe ICE, DH Wanda), czy zły stan torowiska (Muzeum Lotnictwa).

Przyjrzyjmy się ja wyglądały spóźnienia na Stradomiu i na Cystersów w zależności od godziny.

▼ Stradom ▼

▼ Cystersów ▼

Ciekawe, że w stronę Czerwonych Maków szczyt występuje rano, a w stronę Osiedla Piastów po południu. Przynajmniej na tych dwóch przystankach, bo na całej trasie zależność ta jest jeszcze bardziej interesująca.

Ponieważ różnice w przeciętnych spóźnieniach pomiędzy przystankami są znacznie większe niż pomiędzy różnymi porami dnia, mediana została przeskalowana do wartości [0,1]. Inaczej mówiąc, wykresy powyżej powstały jakby przez przeklejenie samych tylko kolorów z wykresów takich jak ten poniżej, zrobionych dla każdego przystanku osobno:

Najbardziej uderza oczywiście, że szczyt popołudniowy zaczyna się na trasie do Czerwonych Maków gdzieś między Rondem Mogilskim, a Dworcem Głównym. Na poprzednich przystankach spóźnienia ok. 17:00 nie odstają tak mocno. W drugą stronę taka sytuacja ma z kolei miejsce rano, choć kontrast nie jest tak wyraźny. Interesujące byłoby porównanie tych wyników z danymi o przepływie pasażerów.

Na koniec jedna obserwacja na temat porannego szczytu spóźnień. Przyszło mi do głowy następujace wyjaśnienie faktu, że występuje on później niż zwczajowo kojarzony poranny szczyt komunikacyjny: ok. 10 tramwaje mają bardziej wyśrubowane rozkłady, ale ruch nie zamiera aż tak szybko. Bardzo łatwo było tą hipotezę obalić: planowane czasy przejazdów pomiędzy przystankami generalnie nie różnią się w ciągu dnia (ewentualnie na niektórych odcinkach, późno w nocy).

fit_plot <- function(mqgamObject, type){
  tau <- as.numeric(names(mqgamObject$fit))#used quantiles

  resids <- qdo(mqgamObject, fun = residuals, type = 'response') %>%
    as.data.frame()#model residuals
  
  if(type == 'seq_num'){
    explan <- mqgamObject$data$seq_num 
    subtitle <- 'Dopasowanie względem numeru przystanku'
    xlab <- 'Numer przystanku na trasie'
  }
  if(type == 'minutefs'){
    explan <- hour(round_date(mqgamObject$data$plannedTime, unit = 'hours'))
    subtitle <- 'Dopasowanie względem numeru godziny'
    xlab <- 'Godzina'
  }
  
  est <- as.data.frame(resids <= 0)#proportion of values below fitted quantile. Ideally will equal quantile value
  names(est) <- tau
  
  est <- cbind(explan, est) %>%
    group_by(explan) %>% 
    summarise_all(funs(mean)) %>% #proportion of residuals under 0
    gather(key = 'quantile', value = 'est', -explan)

  upper <- as.data.frame(resids <= 1/6) # 10secs upper tolerance
  names(upper) <- tau
  upper <- cbind(explan, upper) %>%
    group_by(explan) %>% summarise_all(funs(mean)) %>%
    gather(key = 'quantile', value = 'upper', -explan)

  lower <- as.data.frame(resids <= -1/6)
  names(lower) <- tau
  lower <- cbind(explan, lower) %>%
    group_by(explan) %>% summarise_all(funs(mean)) %>%
    gather(key = 'quantile', value = 'lower', -explan)

  data <- left_join(est, upper, by = c('explan','quantile')) %>%
    left_join(lower, by = c('explan', 'quantile'))

  ggplot(data, aes(x = explan, y = est, ymin = lower, ymax = upper)) +
    geom_line(aes(color = quantile)) + geom_ribbon(alpha = 0.2) +
    geom_hline(aes(yintercept = as.numeric(quantile)), color = swatch()[1], linetype = 2) +
    facet_grid(quantile ~ ., scales = 'free_y') + 
    scale_x_discrete(limits = sort(explan)) + 
    theme(legend.position = 'none') + 
    labs(title = paste('Linia:', mqgamObject$data$number[1], '| Kierunek:', mqgamObject$data$direction[1]), subtitle = subtitle, x = xlab, y = 'Empiryczne kwantyle') 
}

rho <- function(u, tau) u*(tau - (u < 0))

R1 <- function(mqgamObject){

  tau <- as.numeric(names(mqgamObject$fit))
  output <- numeric(length(tau))
  names(output) <- names(mqgamObject$fit)
  
  sample_quantiles <- quantile(mqgamObject$model$delay, probs = tau) %>% t() %>% as.data.frame()
  null_residuals <- sapply(sample_quantiles, function(x) mqgamObject$model$delay - x) %>% as.data.frame()
  
  #extracted <- qdo(mqgamObject, fun = predict, type = 'response')
  
  for(i in seq_along(tau)){
  Vhat <- sum(rho(mqgamObject$model$delay - mqgamObject$fit[[i]]$fitted.values, tau[i]))
  V <- sum(rho(null_residuals[[i]], tau[i]))
  output[i] <- 1-Vhat/V
    }
  output
}

Pozostałe linie

Rozwiń ▼

W dalszej części pracy będę umieszczał same tylko wykresy, bez komentarza. Wszystkie wykresy spóźnienia względem przystanków są ustawione na 16:30.

Wydzielmy najpierw zbiory do modelowania poszczególnych linii.

lines <- unique(line_ranking$number)

lin_dir_asgn <- function(x){
  dirs <- filter(line_ranking, number == x) %>% `$`(direction)
  
  a <- dplyr::filter(report, number == x, direction == dirs[[1]])
  a$delay <- jitter(as.numeric(a$delay), amount = 0.5)
  a$delay[a$delay < 0] <- 0
  
  b <- dplyr::filter(report, number == x, direction == dirs[[2]])
  b$delay <- jitter(as.numeric(b$delay), amount = 0.5)
  b$delay[b$delay < 0] <- 0
  
  assign(paste0('df', x, 'a'), a, pos = 1)
  assign(paste0('df', x, 'b'), b, pos = 1)
  
  return(NULL)
}

purrr::walk(lines, lin_dir_asgn)

W przypadku linii 52 wymiar bazy spaljnów był dobrany na podstawie cross-validacji modelu z przystankami jako zmienną jakościową. Dla następnych linii przyjąłem inną strategię. Ponieważ przeprowadzenie tak wielu testów trawłoby zbyt długo dla każdej linii podiz zbiór na część testową i trenigową, a następnie sprawdziłem wartość funkcji kosztu dla różnych wymiarów bazy. W przypadkach spornych (różne wyniki dla róznych kwantyli) decydowałem na podstawie wykresów dopasowań brzegowych. Wyniki modelowanie poniżej:

• Linia 1
• Linia 2
• Linia 3
• Linia 4
• Linia 5
• Linia 6
• Linia 8
• Linia 9
• Linia 10
• Linia 11
• Linia 13
• Linia 14
• Linia 16
• Linia 18
• Linia 19
• Linia 20
• Linia 21
• Linia 22
• Linia 24
• Linia 44
• Linia 50

Linia 1

qgam_model_1A <- mqgam(delay ~ ti(seq_num, k = 27) + ti(minutefs, k = 7) + ti(seq_num, minutefs), data = df1a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_1A, fun = summary) #hidden because of lenght

qgam_model_1B <- mqgam(delay ~ ti(seq_num, k = 27) + ti(minutefs, k = 8) + ti(seq_num, minutefs), data = df1b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_1B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0	0.06	0.09	0.12	0.15	0.16	0.18	0.22

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.01	0.05	0.08	0.1	0.12	0.13	0.15	0.16

Linia 2

qgam_model_2A <- mqgam(delay ~ ti(seq_num, k = 7) + ti(minutefs, k = 7) + ti(seq_num, minutefs), data = df2a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_2A, fun = summary) #hidden because of lenght

qgam_model_2B <- mqgam(delay ~ ti(seq_num, k = 6) + ti(minutefs, k = 8) + ti(seq_num, minutefs), data = df2b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_2B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
-0.03	-0.02	0.01	0.03	0.04	0.07	0.12	0.17

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0	0.08	0.14	0.18	0.23	0.29	0.33	0.38

Linia 3

qgam_model_3A <- mqgam(delay ~ ti(seq_num, k = 18) + ti(minutefs, k = 8) + ti(seq_num, minutefs), data = df3a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_3A, fun = summary) #hidden because of lenght

qgam_model_3B <- mqgam(delay ~ ti(seq_num, k = 22) + ti(minutefs, k = 8) + ti(seq_num, minutefs), data = df3b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_3B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0	0.05	0.09	0.12	0.15	0.17	0.19	0.23

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.06	0.12	0.16	0.19	0.21	0.23	0.26	0.29

Linia 4

qgam_model_4A <- mqgam(delay ~ ti(seq_num, k = 28) + ti(minutefs, k = 13) + ti(seq_num, minutefs), data = df4a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_4A, fun = summary) #hidden because of lenght

qgam_model_4B <- mqgam(delay ~ ti(seq_num, k = 28) + ti(minutefs, k = 13) + ti(seq_num, minutefs), data = df4b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_4B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
-0.03	0.02	0.06	0.09	0.12	0.15	0.17	0.21

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.05	0.11	0.14	0.15	0.16	0.16	0.16	0.16

Linia 5

qgam_model_5A <- mqgam(delay ~ ti(seq_num, k = 22, bs = 'cs') + ti(minutefs, k = 14) + ti(seq_num, minutefs), data = df5a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_5A, fun = summary) #hidden because of lenght

qgam_model_5B <- mqgam(delay ~ ti(seq_num, k = 19) + ti(minutefs, k = 14) + ti(seq_num, minutefs), data = df5b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_5B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
-0.01	0.01	0.04	0.06	0.08	0.1	0.08	0.09

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.07	0.13	0.16	0.17	0.18	0.18	0.19	0.19

Linia 6

qgam_model_6A <- mqgam(delay ~ ti(seq_num, k = 19) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df6a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_6A, fun = summary) #hidden because of lenght

qgam_model_6B <- mqgam(delay ~ ti(seq_num, k = 19) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df6b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_6B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.05	0.13	0.17	0.21	0.25	0.27	0.3	0.29

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.07	0.11	0.13	0.14	0.15	0.16	0.17	0.18

Linia 8

qgam_model_8A <- mqgam(delay ~ ti(seq_num, k = 22) + ti(minutefs, k = 14) + ti(seq_num, minutefs), data = df8a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_1A, fun = summary) #hidden because of lenght

qgam_model_8B <- mqgam(delay ~ ti(seq_num, k = 22) + ti(minutefs, k = 14) + ti(seq_num, minutefs), data = df8b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_1B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.01	0.07	0.12	0.17	0.21	0.24	0.27	0.3

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.02	0.1	0.14	0.19	0.23	0.25	0.28	0.3

Linia 9

Pomijam z powodu zbyt dużego odsetka niezarejestrowanych odjazdów.

Linia 10

qgam_model_10A <- mqgam(delay ~ ti(seq_num, k = 24) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df10a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_10A, fun = summary) #hidden because of lenght

qgam_model_10B <- mqgam(delay ~ ti(seq_num, k = 25) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df10b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_10B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.09	0.15	0.19	0.22	0.24	0.25	0.26	0.24

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0	0.06	0.1	0.13	0.16	0.17	0.18	0.18

Linia 11

qgam_model_11A <- mqgam(delay ~ ti(seq_num, k = 15) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df11a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_11A, fun = summary) #hidden because of lenght

qgam_model_11B <- mqgam(delay ~ ti(seq_num, k = 17) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df11b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_11B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
-0.04	-0.01	0.02	0.05	0.08	0.11	0.15	0.14

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.09	0.18	0.24	0.29	0.34	0.37	0.4	0.4

Linia 13

qgam_model_13A <- mqgam(delay ~ ti(seq_num, k = 25) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df13a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_13A, fun = summary) #hidden because of lenght

qgam_model_13B <- mqgam(delay ~ ti(seq_num, k = 25) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df13b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_13B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.04	0.1	0.14	0.17	0.2	0.21	0.24	0.27

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.06	0.11	0.15	0.18	0.2	0.22	0.25	0.28

Linia 14

qgam_model_14A <- mqgam(delay ~ ti(seq_num, k = 31) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df14a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_14A, fun = summary) #hidden because of lenght

qgam_model_14B <- mqgam(delay ~ ti(seq_num, k = 33) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df14b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_14B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.02	0.1	0.14	0.18	0.21	0.24	0.26	0.28

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.03	0.09	0.12	0.14	0.15	0.16	0.17	0.19

Linia 16

Pomijam z powodu zbyt dużego odsetka niezarejestrowanych odjazdów.

Linia 18

qgam_model_18A <- mqgam(delay ~ ti(seq_num, k = 23) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df18a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_18A, fun = summary) #hidden because of lenght

qgam_model_18B <- mqgam(delay ~ ti(seq_num, k = 23) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df18b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_18B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
-0.01	0.04	0.07	0.09	0.12	0.12	0.13	0.15

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
-0.01	0.04	0.08	0.11	0.14	0.15	0.15	0.17

Linia 19

qgam_model_19A <- mqgam(delay ~ ti(seq_num, k = 17) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df19a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_19A, fun = summary) #hidden because of lenght

qgam_model_19B <- mqgam(delay ~ ti(seq_num, k = 16) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df19b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_19B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.14	0.17	0.19	0.2	0.22	0.23	0.24	0.27

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.04	0.12	0.17	0.21	0.25	0.28	0.3	0.32

Linia 20

qgam_model_20A <- mqgam(delay ~ ti(seq_num, k = 16) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df20a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_20A, fun = summary) #hidden because of lenght

qgam_model_20B <- mqgam(delay ~ ti(seq_num, k = 16) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df20b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_20B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.02	0.08	0.11	0.13	0.14	0.12	0.1	0.1

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.1	0.16	0.2	0.23	0.25	0.28	0.3	0.33

Linia 21

Pomijam ze względu na zbyt duży odsetek niezarejestrowanych odjazdów.

Linia 22

qgam_model_22A <- mqgam(delay ~ ti(seq_num, k = 34) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df22a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_22A, fun = summary) #hidden because of lenght

qgam_model_22B <- mqgam(delay ~ ti(seq_num, k = 34) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df22b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_22B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.03	0.09	0.13	0.16	0.18	0.2	0.21	0.22

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.07	0.12	0.15	0.18	0.19	0.2	0.21	0.23

Linia 24

qgam_model_24A <- mqgam(delay ~ ti(seq_num, k = 23) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df24a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_24A, fun = summary) #hidden because of lenght

qgam_model_24B <- mqgam(delay ~ ti(seq_num, k = 25) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df24b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_24B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.07	0.1	0.11	0.11	0.12	0.12	0.13	0.17

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.09	0.14	0.17	0.2	0.22	0.25	0.28	0.31

Linia 44

qgam_model_44A <- mqgam(delay ~ ti(seq_num, k = 23) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df44a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_44A, fun = summary) #hidden because of lenght

qgam_model_44B <- mqgam(delay ~ ti(seq_num, k = 23) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df44b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_44B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
-0.01	0.04	0.07	0.09	0.12	0.13	0.14	0.19

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.04	0.13	0.17	0.2	0.23	0.26	0.28	0.31

Linia 50

qgam_model_50A <- mqgam(delay ~ ti(seq_num, k = 15) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df50a, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_50A, fun = summary) #hidden because of lenght

qgam_model_50B <- mqgam(delay ~ ti(seq_num, k = 15) + ti(minutefs, k = 15) + ti(seq_num, minutefs), data = df50b, qu = seq(0.2, 0.9, 0.1))
qdo(qgam_model_50B, fun = summary) #hidden because of lenght

Diagnostyka

Rozwiń ▼

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.05	0.11	0.13	0.15	0.16	0.16	0.16	0.16

0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
0.18	0.22	0.25	0.27	0.28	0.29	0.29	0.27