Help us understand the problem. What is going on with this article?

kaggleのkernels学習ノート~NYC Taxi EDA - Update: The fast & the curious~

More than 1 year has passed since last update.

題名の通りkaggleのkernelsを和訳したものです。自分の勉強の備忘録をまとめていきます。そして、データ分析に携わる誰かのお役に立てれば幸いです。和訳とか得意じゃないので間違っていたらごめんなさい。
今回はHitesh palamadaHeadさんのNYC Taxi EDA - Update: The fast & the curiousを和訳していきます。途中から睡魔に襲われて、訳がグチャグチャですが、少しづづ手直しします。また、コードブロックを見れば著者の行いたいことはわかるかと思います…

1 はじめに

これは、tidy Rggplot2を使って、New York City Taxi Trip Duration competitionの探索的データ分析です。
この遊び場の目的は、ひと走り時の座標やピックアップの日時などの特徴量に基づいてニューヨークのタクシー移動時間を予測することです。データは150万回の訓練データ(../input/train.csv)と630kのテストデータ(../input/test.csv)です。各行には、1回のタクシー移動情報が含まれています。
このノートブックでは、最初に元のデータを調べて可視化し、新しい特徴量を作り、潜在的な異常値を調べます。その後、NYCの天気と理論的に最速のルートを含む2つの外部データセットを追加します。これらのデータセット内の新しい特徴量と、ターゲットであるtrip_durationへの影響を視覚化して分析します。最後に、この課題を分類問題と考え、XGBoostを使って基本的な予測を行い、このノートブックを完成させます。
このノートブックが、この挑戦を始める際にあなたの手助けになることを願っています。新しい特徴量や視覚化のためのアイデアを考える余地はたくさんあります。特に、分類アプローチと最終モデルは、パフォーマンスを向上させるために改善し、最適化するための基本的な出発点にすぎません。いつものように、フィードバック、質問、建設的な批判をお待ちしております。
もともと、この分析には、当時のカーネル環境のメモリ/実行時の制限のために、いくつかの隠された値が含まれていました。これは、いくつかの補助的な視覚化を目隠しをしなければなりませんでした。分析の流れにあまり影響を与えないプロットだけを削除しようとしました。対応するコードは、すべてこのノートブックに含まれており、これらの隠しプロットを簡単に再表示できます。さて、新しく改良されたカーネルでは、それらのプロットは "公式"バージョンに戻っています。
最後に、皆さんにこのカーネルを見たり、アップしたり、コメントしたりすることに感謝したいと思います!私はまだKaggleの初心者ですし、あなたのサポートは私にとっては実りのあるものです。また、NockedDownのおかげで、タイトルが新しく更新されました。さぁ、始めましょう。

1.1 ライブラリのロード

R
library('ggplot2') # visualisation
library('scales') # visualisation
library('grid') # visualisation
library('RColorBrewer') # visualisation
library('corrplot') # visualisation
library('alluvial') # visualisation
library('dplyr') # data manipulation
library('readr') # input/output
library('data.table') # data manipulation
library('tibble') # data wrangling
library('tidyr') # data wrangling
library('stringr') # string manipulation
library('forcats') # factor manipulation
library('lubridate') # date and time
library('geosphere') # geospatial locations
library('leaflet') # maps
library('leaflet.extras') # maps
library('maps') # maps
library('xgboost') # modelling
library('caret') # modelling

R Cookbooksのmultiplot functionを使用して、マルチパネルプロットを作成しました。

R
# Define multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

1.2 データの読み込み

R
train <- as.tibble(fread('../input/nyc-taxi-trip-duration/train.csv'))
test <- as.tibble(fread('../input/nyc-taxi-trip-duration/test.csv'))
sample_submit <- as.tibble(fread('../input/nyc-taxi-trip-duration/sample_submission.csv'))

data.tableのfreadを使用してデータの読み込みを高速化します。

1.3 ファイル構造と内容

summaryglimpseを使って、まずは訓練データを見てみましょう。

R
summary(train)
##       id              vendor_id     pickup_datetime    dropoff_datetime  
##  Length:1458644     Min.   :1.000   Length:1458644     Length:1458644    
##  Class :character   1st Qu.:1.000   Class :character   Class :character  
##  Mode  :character   Median :2.000   Mode  :character   Mode  :character  
##                     Mean   :1.535                                        
##                     3rd Qu.:2.000                                        
##                     Max.   :2.000                                        
##  passenger_count pickup_longitude  pickup_latitude dropoff_longitude
##  Min.   :0.000   Min.   :-121.93   Min.   :34.36   Min.   :-121.93  
##  1st Qu.:1.000   1st Qu.: -73.99   1st Qu.:40.74   1st Qu.: -73.99  
##  Median :1.000   Median : -73.98   Median :40.75   Median : -73.98  
##  Mean   :1.665   Mean   : -73.97   Mean   :40.75   Mean   : -73.97  
##  3rd Qu.:2.000   3rd Qu.: -73.97   3rd Qu.:40.77   3rd Qu.: -73.96  
##  Max.   :9.000   Max.   : -61.34   Max.   :51.88   Max.   : -61.34  
##  dropoff_latitude store_and_fwd_flag trip_duration    
##  Min.   :32.18    Length:1458644     Min.   :      1  
##  1st Qu.:40.74    Class :character   1st Qu.:    397  
##  Median :40.75    Mode  :character   Median :    662  
##  Mean   :40.75                       Mean   :    959  
##  3rd Qu.:40.77                       3rd Qu.:   1075  
##  Max.   :43.92                       Max.   :3526282

glimpse(train)
## Observations: 1,458,644
## Variables: 11
## $ id                 <chr> "id2875421", "id2377394", "id3858529", "id3...
## $ vendor_id          <int> 2, 1, 2, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2, 2...
## $ pickup_datetime    <chr> "2016-03-14 17:24:55", "2016-06-12 00:43:35...
## $ dropoff_datetime   <chr> "2016-03-14 17:32:30", "2016-06-12 00:54:38...
## $ passenger_count    <int> 1, 1, 1, 1, 1, 6, 4, 1, 1, 1, 1, 4, 2, 1, 1...
## $ pickup_longitude   <dbl> -73.98215, -73.98042, -73.97903, -74.01004,...
## $ pickup_latitude    <dbl> 40.76794, 40.73856, 40.76394, 40.71997, 40....
## $ dropoff_longitude  <dbl> -73.96463, -73.99948, -74.00533, -74.01227,...
## $ dropoff_latitude   <dbl> 40.76560, 40.73115, 40.71009, 40.70672, 40....
## $ store_and_fwd_flag <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N"...
## $ trip_duration      <int> 455, 663, 2124, 429, 435, 443, 341, 1551, 2...

次はテストデータです。

R
summary(test)
##       id              vendor_id     pickup_datetime    passenger_count
##  Length:625134      Min.   :1.000   Length:625134      Min.   :0.000  
##  Class :character   1st Qu.:1.000   Class :character   1st Qu.:1.000  
##  Mode  :character   Median :2.000   Mode  :character   Median :1.000  
##                     Mean   :1.535                      Mean   :1.662  
##                     3rd Qu.:2.000                      3rd Qu.:2.000  
##                     Max.   :2.000                      Max.   :9.000  
##  pickup_longitude  pickup_latitude dropoff_longitude dropoff_latitude
##  Min.   :-121.93   Min.   :37.39   Min.   :-121.93   Min.   :36.60   
##  1st Qu.: -73.99   1st Qu.:40.74   1st Qu.: -73.99   1st Qu.:40.74   
##  Median : -73.98   Median :40.75   Median : -73.98   Median :40.75   
##  Mean   : -73.97   Mean   :40.75   Mean   : -73.97   Mean   :40.75   
##  3rd Qu.: -73.97   3rd Qu.:40.77   3rd Qu.: -73.96   3rd Qu.:40.77   
##  Max.   : -69.25   Max.   :42.81   Max.   : -67.50   Max.   :48.86   
##  store_and_fwd_flag
##  Length:625134     
##  Class :character  
##  Mode  :character  

glimpse(test)
## Observations: 625,134
## Variables: 9
## $ id                 <chr> "id3004672", "id3505355", "id1217141", "id2...
## $ vendor_id          <int> 1, 1, 1, 2, 1, 1, 1, 1, 2, 2, 1, 2, 1, 2, 1...
## $ pickup_datetime    <chr> "2016-06-30 23:59:58", "2016-06-30 23:59:53...
## $ passenger_count    <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 4, 1, 1, 1, 1...
## $ pickup_longitude   <dbl> -73.98813, -73.96420, -73.99744, -73.95607,...
## $ pickup_latitude    <dbl> 40.73203, 40.67999, 40.73758, 40.77190, 40....
## $ dropoff_longitude  <dbl> -73.99017, -73.95981, -73.98616, -73.98643,...
## $ dropoff_latitude   <dbl> 40.75668, 40.65540, 40.72952, 40.73047, 40....
## $ store_and_fwd_flag <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N"...

わかったことは下記の通りです。

  • vendor_idは1か2の値しかないので、おそらく2つのタクシー会社を区別する
  • pickup_datetimedropoff_datetimeは、分析可能な形に再加工する必要がある日時の組み合わせ
  • passenger_countは両データセットで中央値が1で最大値が9
  • pickup / dropoff_longitute / latituteは、メーターがアクティブ/非アクティブになった座標を表す。
  • store_and_fwd_flagは、データがベンダー(「N」)に直ちに送信されたか、またはサーバー(「Y」)に接続されていなかったためにタクシーのメモリーに格納されたかどうかを示すフラグ。 おそらく、受信の悪い特定の地域と相関関係があるのかもしれない?
  • trip_durationは訓練データのターゲットの特徴量で秒単位で計測される

1.4 欠損値の確認

欠損値を知ることは、データの不明瞭さを示すため重要です。 ほんの数例に基づいて推論を行うことは、賢明ではありません。さらに、欠損値が含まれている場合には、多くのモデリング手順が崩れ、対応する行を完全に削除するか、何らかの方法で補完する必要があります。幸いなことに、私たちのデータは完全であり、欠損値はありません。

R
sum(is.na(train))
## [1] 0

sum(is.na(test))
## [1] 0

1.5 訓練データとテストデータを結合

最終的なモデリングに備えて、訓練データとテストデータを1つにまとめます。一般的に、テストデータを調べすぎないことをお勧めします。これは、モデルをオーバーフィットさせるリスクがあるためです。ただし、2つのデータセット間の単純な一貫性チェックが有利な場合はあります。

R
combine <- bind_rows(train %>% mutate(dset = "train"), 
                     test %>% mutate(dset = "test",
                                     dropoff_datetime = NA,
                                     trip_duration = NA))
combine <- combine %>% mutate(dset = factor(dset))

1.6 特徴量の再形成

以下の分析では、日付と時刻を文字型から日付型に変換します。vendor_idもファクタ型として再コード化します。これにより、これらの特徴量を含む関係を視覚化することが容易になります。

R
train <- train %>%
  mutate(pickup_datetime = ymd_hms(pickup_datetime),
         dropoff_datetime = ymd_hms(dropoff_datetime),
         vendor_id = factor(vendor_id),
         passenger_count = factor(passenger_count))

1.7 一貫性の確認

trip_durationspickup_datetimedropoff_datetimeの間隔と一致するかどうかを調べます。おそらく前者は後者から直接計算されたものですが、あなたは決して知らないでしょう。以下の2つの間隔が一致していない場合、"TRUE"が表示されます。

R
train %>%
  mutate(check = abs(int_length(interval(dropoff_datetime,pickup_datetime)) + trip_duration) > 0) %>%
  select(check, pickup_datetime, dropoff_datetime, trip_duration) %>%
  group_by(check) %>%
  count()

## # A tibble: 1 x 2
## # Groups:   check [1]
##   check       n
##   <lgl>   <int>
## 1 F     1458644

問題なさそうですね。

2 特徴量の視覚化

特徴量の分布の視覚化とそれらの関係は、データセットを理解する上で重要であり、しばしば新しい探求に繋がります。微妙な傾向や相関関係にも気づくために、可能な限り多くの異なる視点からデータを調べることを常にお勧めします。このセクションでは、まず個々のデータ特徴量の分布を見ていきます。
NYCの地図から始め、可能な番号のピックアップ座標をオーバーレイして、問題の場所と距離の概要を表示します。この視覚化では、インタラクティブマップ用のさまざまな素晴らしいなツールを含むleafletパッケージを使用します。 この地図では、ピックアップの場所をズームしてパンすることができます:
※この和訳記事ではインタラクティブではありません。

R
set.seed(1234)
foo <- sample_n(train, 8e3)
leaflet(data = foo) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircleMarkers(~ pickup_longitude, ~pickup_latitude, radius = 1,
                   color = "blue", fillOpacity = 0.3)

スクリーンショット 0030-09-07 23.04.47.png
ほとんどの移動は、マンハッタンでしか行われていなかったことが判明しました。もう一つ注目すべきホットスポットはJFK空港で、市内の南東に伸びています。マップは、いくつかの分布がどのように見えるかを教えてくれます。ターゲットの特徴量であるtrip_durationをプロットすることから始めましょう。

R
train %>%
  ggplot(aes(trip_duration)) +
  geom_histogram(fill = "red", bins = 150) +
  scale_x_log10() +
  scale_y_sqrt()

unnamed-chunk-13-1.png
わかったことは下記の通りです。

  • 大部分の移動は、1000秒、すなわち約17分にピークがある対数正規に見える滑らかな分布に従う。
  • 10秒未満の短い不審な移動がいくつかあります。
  • さらに、trip_durationには、奇妙なデルタ型のピークが1e5秒の直前にあり、それより右のほうにもあります。
R
train %>%
  arrange(desc(trip_duration)) %>%
  select(trip_duration, pickup_datetime, dropoff_datetime, everything()) %>%
  head(10)
## # A tibble: 10 x 11
##    trip_duration pickup_datetime     dropoff_datetime    id      vendor_id
##            <int> <dttm>              <dttm>              <chr>   <fct>    
##  1       3526282 2016-02-13 22:46:52 2016-03-25 18:18:14 id0053… 1        
##  2       2227612 2016-01-05 06:14:15 2016-01-31 01:01:07 id1325… 1        
##  3       2049578 2016-02-13 22:38:00 2016-03-08 15:57:38 id0369… 1        
##  4       1939736 2016-01-05 00:19:42 2016-01-27 11:08:38 id1864… 1        
##  5         86392 2016-02-15 23:18:06 2016-02-16 23:17:58 id1942… 2        
##  6         86391 2016-05-31 13:00:39 2016-06-01 13:00:30 id0593… 2        
##  7         86390 2016-05-06 00:00:10 2016-05-07 00:00:00 id0953… 2        
##  8         86387 2016-06-30 16:37:52 2016-07-01 16:37:39 id2837… 2        
##  9         86385 2016-06-23 16:01:45 2016-06-24 16:01:30 id1358… 2        
## 10         86379 2016-05-17 22:22:56 2016-05-18 22:22:35 id2589… 2        
## # ... with 6 more variables: passenger_count <fct>,
## #   pickup_longitude <dbl>, pickup_latitude <dbl>,
## #   dropoff_longitude <dbl>, dropoff_latitude <dbl>,
## #   store_and_fwd_flag <chr>

これらの記録は24時間以上の移動に対応し、最大でほぼ12日間です。ラッシュアワーは悪いことがあることは知っていますが、この値は少し信じられません。
1年を通して、pickup_datetimeとdropoff_datetimeの分布は次のようになります。

R
p1 <- train %>%
  ggplot(aes(pickup_datetime)) +
  geom_histogram(fill = "red", bins = 120) +
  labs(x = "Pickup dates")

p2 <- train %>%
  ggplot(aes(dropoff_datetime)) +
  geom_histogram(fill = "blue", bins = 120) +
  labs(x = "Dropoff dates")

layout <- matrix(c(1,2),2,1,byrow=FALSE)
multiplot(p1, p2, layout=layout)

unnamed-chunk-15-1.png
かなり均質で、2016年1月から7月までの半年間をカバーしています。1月下旬には2月初旬ごろに興味深い落ち込みがあります。

R
train %>%
  filter(pickup_datetime > ymd("2016-01-20") & pickup_datetime < ymd("2016-02-10")) %>%
  ggplot(aes(pickup_datetime)) +
  geom_histogram(fill = "red", bins = 120)

unnamed-chunk-16-1.png
NYCの冬なので、嵐や天候が悪いのかもしれません?このようなイベントは、おそらく便利な外部データセットを考慮して分析する必要があります。上記のプロットでは、すでに移動の回数に応じて毎日、毎週の調整が行われています。さまざまなコンポーネントでマルチプロットパネルを作成して、passenger_countvendor_idの分布とともにこれらのバリエーションを調べてみましょう。

R
p1 <- train %>%
  group_by(passenger_count) %>%
  count() %>%
  ggplot(aes(passenger_count, n, fill = passenger_count)) +
  geom_col() +
  scale_y_sqrt() +
  theme(legend.position = "none")

p2 <- train %>%
  ggplot(aes(vendor_id, fill = vendor_id)) +
  geom_bar() +
  theme(legend.position = "none")

p3 <- train %>%
  ggplot(aes(store_and_fwd_flag)) +
  geom_bar() +
  theme(legend.position = "none") +
  scale_y_log10()

p4 <- train %>%
  mutate(wday = wday(pickup_datetime, label = TRUE)) %>%
  group_by(wday, vendor_id) %>%
  count() %>%
  ggplot(aes(wday, n, colour = vendor_id)) +
  geom_point(size = 4) +
  labs(x = "Day of the week", y = "Total number of pickups") +
  theme(legend.position = "none")

p5 <- train %>%
  mutate(hpick = hour(pickup_datetime)) %>%
  group_by(hpick, vendor_id) %>%
  count() %>%
  ggplot(aes(hpick, n, color = vendor_id)) +
  geom_point(size = 4) +
  labs(x = "Hour of the day", y = "Total number of pickups") +
  theme(legend.position = "none")

layout <- matrix(c(1,2,3,4,5,5),3,2,byrow=TRUE)
multiplot(p1, p2, p3, p4, p5, layout=layout)

unnamed-chunk-17-1.png

R
p1 <- 1; p2 <- 1; p3 <- 1; p4 <- 1; p5 <- 1

わかったことは下記の通りです。

  • 0人、または7人から9人の乗客がいる移動がありますが、まれな外れ値でしょう。
R
train %>%
  group_by(passenger_count) %>%
  count()

## # A tibble: 10 x 2
## # Groups:   passenger_count [10]
##    passenger_count       n
##    <fct>             <int>
##  1 0                    60
##  2 1               1033540
##  3 2                210318
##  4 3                 59896
##  5 4                 28404
##  6 5                 78088
##  7 6                 48333
##  8 7                     3
##  9 8                     1
## 10 9                     1
  • 大部分の移動は、1人の乗客しかいませんでした。2人の乗客が2番目に多いようです。
  • 乗客数が増えるにつれて、5人から6人の乗客に大きなピークがありますが、3人から4人は順調に減少しています。
  • ベンダー2は、ベンダー1よりもこのデータセットでかなり多くの乗客データを有します(対数y軸に注意してください)。これは、毎週の曜日に当てはまります。
  • 月曜日が最も落ち着いており、金曜日がとても忙しいという興味深いパターンを見つけます。これは2つの異なるベンダーで同じですが、vendor_id == 2は著しく高い照射数を示しています。
  • 直観的に予想しているように、早朝の時間帯には急降下があります。また、2つのベンダーの間に大きな違いはありません。午後4時ごろ、別の落ち込みがあり、夕方に向かっ増えていきます。
  • store_and_fwd_flagは、移動データがベンダー(「N」)に直ちに送信されたか、またはサーバー(「Y」)への接続がなかったためにタクシーのメモリーに保持されたかどうかを示します(y軸は対数です。注意してください)。
R
train %>%
  group_by(store_and_fwd_flag) %>%
  count()

## # A tibble: 2 x 2
## # Groups:   store_and_fwd_flag [2]
##   store_and_fwd_flag       n
##   <chr>                <int>
## 1 N                  1450599
## 2 Y                     8045

これらの数値は、直ちに送信されない移動の約半分に相当します。1日の1時間あたりの移動量は、月によって多少異なり、曜日に強く依存します。

R
p1 <- train %>%
  mutate(hpick = hour(pickup_datetime),
         Month = factor(month(pickup_datetime, label = TRUE))) %>%
  group_by(hpick, Month) %>%
  count() %>%
  ggplot(aes(hpick, n, color = Month)) +
  geom_line(size = 1.5) +
  labs(x = "Hour of the day", y = "count")

p2 <- train %>%
  mutate(hpick = hour(pickup_datetime),
         wday = factor(wday(pickup_datetime, label = TRUE))) %>%
  group_by(hpick, wday) %>%
  count() %>%
  ggplot(aes(hpick, n, color = wday)) +
  geom_line(size = 1.5) +
  labs(x = "Hour of the day", y = "count")

layout <- matrix(c(1,2),2,1,byrow=FALSE)

multiplot(p1, p2, layout=layout)

plot_zoom_png.png

R
p1 <- 1; p2 <- 1

わかったこと。

  • 1月と6月は移動回数が少なく、3月と4月は多いです。この傾向は両方のベンダーIDで確認できます。
  • 週末(土曜日と日曜日、金曜日まで)は、早朝の移動回数が高くなりますが、朝5時から10時の移動回数は少なくなります。これはNYCの営業日と週末の夜の生活の対比に起因する可能性があります。さらに、移動回数は日曜の夕方/夜に減ります。

最後に、pickup/dropoff latitudes and longitudesを視覚化します。

R
p1 <- train %>%
  filter(pickup_longitude > -74.05 & pickup_longitude < -73.7) %>%
  ggplot(aes(pickup_longitude)) +
  geom_histogram(fill = "red", bins = 40)

p2 <- train %>%
  filter(dropoff_longitude > -74.05 & dropoff_longitude < -73.7) %>%
  ggplot(aes(dropoff_longitude)) +
  geom_histogram(fill = "blue", bins = 40)

p3 <- train %>%
  filter(pickup_latitude > 40.6 & pickup_latitude < 40.9) %>%
  ggplot(aes(pickup_latitude)) +
  geom_histogram(fill = "red", bins = 40)

p4 <- train %>%
  filter(dropoff_latitude > 40.6 & dropoff_latitude < 40.9) %>%
  ggplot(aes(dropoff_latitude)) +
  geom_histogram(fill = "blue", bins = 40)

layout <- matrix(c(1,2,3,4),2,2,byrow=FALSE)
multiplot(p1, p2, p3, p4, layout=layout)

unnamed-chunk-21-1.png

R
p1 <- 1; p2 <- 1; p3 <- 1; p4 <- 1; p5 <- 1

ここでは、緯度と経度の値の範囲を制限していました。これは、NYCの外にあるためです。結果として得られた分布は、既に地図上で見たマンハッタンのものと一致しています。
これらはpickup_latitudeの極端な値(外れ値)です。

R
train %>%
  arrange(pickup_latitude) %>%
  select(pickup_latitude, pickup_longitude) %>%
  head(5)

## # A tibble: 5 x 2
##   pickup_latitude pickup_longitude
##             <dbl>            <dbl>
## 1            34.4            -65.8
## 2            34.7            -75.4
## 3            35.1            -71.8
## 4            35.3            -72.1
## 5            36.0            -77.4

train %>%
  arrange(desc(pickup_latitude)) %>%
  select(pickup_latitude, pickup_longitude) %>%
  head(5)

## # A tibble: 5 x 2
##   pickup_latitude pickup_longitude
##             <dbl>            <dbl>
## 1            51.9            -72.8
## 2            44.4            -67.0
## 3            43.9            -71.9
## 4            43.5            -74.2
## 5            43.1            -72.6

偏った分析をしないためにも、これらの値は頭に留めておく必要があります。

3 特徴量の関係

前のセクションでは、主に個々の特徴量の分布について説明しましたが、ここではそれらの特徴量が互いにどのように関連しているか、およびターゲットであるtrip_durationにどのように関連するかを詳しく説明します。

3.1 Pickup date/time vs trip_duration

1日を通じての移動数の変動は、平均旅行時間にどのように影響するのでしょうか?より少ない日と時間はより短い移動につながりますか?ここでは、vendor_idを追加します。さらに、1日のうちに、変動の程度とその不確実性を示すsmoothing layerを追加します。

R
p1 <- train %>%
  mutate(wday = wday(pickup_datetime, label = TRUE)) %>%
  group_by(wday, vendor_id) %>%
  summarise(median_duration = median(trip_duration)/60) %>%
  ggplot(aes(wday, median_duration, color = vendor_id)) +
  geom_point(size = 4) +
  labs(x = "Day of the week", y = "Median trip duration [min]")

p2 <- train %>%
  mutate(hpick = hour(pickup_datetime)) %>%
  group_by(hpick, vendor_id) %>%
  summarise(median_duration = median(trip_duration)/60) %>%
  ggplot(aes(hpick, median_duration, color = vendor_id)) +
  geom_smooth(method = "loess", span = 1/2) +
  geom_point(size = 4) +
  labs(x = "Hour of the day", y = "Median trip duration [min]") +
  theme(legend.position = "none")

layout <- matrix(c(1,2),2,1,byrow=FALSE)
multiplot(p1, p2, layout=layout)

plot_zoom_png.png

わかったことは、

  • 実際には、曜日にはビジネスと同様のパターンがあります。より多くの移動レコードをを持つベンダー2も、ベンダー1よりも一貫して多い回数となっています。vendor_idをモデルに追加して予測する上での重要性をテストする価値があります。
  • 午後早朝にピークを見つけ、午前5時から午後6時と午後8時の間に回数が減る典型的な1日の流れがわかります。平日と時間は、予測するために重要な特徴量であると思わるので、モデルに含める必要があります。

3.2 Passenger count and Vendor vs trip_duration

次の質問は、異なる数の乗客および/または異なる業者が移動の時間と相関しているかどうかです。この問題は、vendor_idsを対比させるファセットラップと共に、passenger_countsのボックスプロットを使用して調べます。

R
train %>%
  ggplot(aes(passenger_count, trip_duration, color = passenger_count)) +
  geom_boxplot() +
  scale_y_log10() +
  theme(legend.position = "none") +
  facet_wrap(~ vendor_id) +
  labs(y = "Trip duration [s]", x = "Number of passengers")

unnamed-chunk-25-1.png
わかったことは、

  • どちらのベンダーも、乗客がいなくても短い時間で移動している。
  • ベンダー1には24時間を超える旅行がすべて含まれていますが、ベンダー2には6人以上の乗客がいる(5つの)レコードがすべて表示され、24時間の制限に近い移動が多くあります。
  • 1人から6人の乗客の平均移動時間は、ベンダー2の場合と非常に似ています。ベンダー1では違いがありますが、小さいものです(y軸は対数です。注意してください)。
R
train %>%
  ggplot(aes(trip_duration, fill = vendor_id)) +
  geom_density(position = "stack") +
  scale_x_log10()

unnamed-chunk-26-1.png
2つのベンダーのtrip_durationの分布の密度を比較すると、中央値は非常に似ていますが、平均は長時間の外れ値を含んでいるベンダー2は偏っていることがわかります。

R
train %>%
  group_by(vendor_id) %>%
  summarise(mean_duration = mean(trip_duration),
            median_duration = median(trip_duration))

## # A tibble: 2 x 3
##   vendor_id mean_duration median_duration
##   <fct>             <dbl>           <dbl>
## 1 1                   845             658
## 2 2                  1059             666

3.3 Store and Forward vs trip_duration

ベンダー1にのみ一時保管データが発生しました。

R
train %>%
  group_by(vendor_id, store_and_fwd_flag) %>%
  count()

## # A tibble: 3 x 3
## # Groups:   vendor_id, store_and_fwd_flag [3]
##   vendor_id store_and_fwd_flag      n
##   <fct>     <chr>               <int>
## 1 1         N                  670297
## 2 1         Y                    8045
## 3 2         N                  780302

unnamed-chunk-29-1.png
私たちは、保存された移動と保存されていない移動の間には大きな違いはないことを見出しました。保存されたものは若干長くなるかもしれませんが、疑わしい長い移動は含まれません。

4 特徴量エンジニアリング

このセクションでは、既存の変数から新しい特徴量を構築し、ターゲット変数のより良い予測子を見つけます。以下の単一のコードブロックに、これらの新しい特徴量をすべて定義し、以下のサブセクションでそれらを研究することにします。これは線形解析には対応していませんが、ノートブックをより読みやすします。
新しい一時的な特徴量(日付、月、日、時間)はpickup_datetimeから得られます。私たちはWikipediaからJFKとLa Guardia空港の座標を得ました。 ブリザード天候は外部気象データに基づいています。

R
jfk_coord <- tibble(lon = -73.778889, lat = 40.639722)
la_guardia_coord <- tibble(lon = -73.872611, lat = 40.77725)

pick_coord <- train %>%
  select(pickup_longitude, pickup_latitude)

drop_coord <- train %>%
  select(dropoff_longitude, dropoff_latitude)

    train$dist <- distCosine(pick_coord, drop_coord)
train$bearing = bearing(pick_coord, drop_coord)

train$jfk_dist_pick <- distCosine(pick_coord, jfk_coord)
train$jfk_dist_drop <- distCosine(drop_coord, jfk_coord)
train$lg_dist_pick <- distCosine(pick_coord, la_guardia_coord)
train$lg_dist_drop <- distCosine(drop_coord, la_guardia_coord)

train <- train %>%
  mutate(speed = dist/trip_duration*3.6,
         date = date(pickup_datetime),
         month = month(pickup_datetime, label = TRUE),
         wday = wday(pickup_datetime, label = TRUE),
         wday = fct_relevel(wday, c("Mon", "Tues", "Wed", "Thurs", "Fri", "Sat", "Sun")),
         hour = hour(pickup_datetime),
         work = (hour %in% seq(8,18)) & (wday %in% c("Mon","Tues","Wed","Thurs","Fri")),
         jfk_trip = (jfk_dist_pick < 2e3) | (jfk_dist_drop < 2e3),
         lg_trip = (lg_dist_pick < 2e3) | (lg_dist_drop < 2e3),
         blizzard = !( (date < ymd("2016-01-22") | (date > ymd("2016-01-29"))) )
         )

4.1 移動の直線距離

ピックアップポイントとドロップオフポイントの座標から、2つのポイント間の直接距離を計算し、それをtrip_durationsと比較することができます。タクシーはカラスではないので、これらの値は可能な最小移動距離に相当します。
これらの距離を計算するために、球面三角法のためのgeosphereパッケージのdistCosine関数を使用しています。この方法は球面上の2点間の最短距離を与えます。この分析の目的のために、地球の形の楕円形の歪みを無視することを選択します。距離vs時間の生の値です。

R
set.seed(4321)

train %>%
  sample_n(5e4) %>%
  ggplot(aes(dist, trip_duration)) +
  geom_point() +
  scale_x_log10() +
  scale_y_log10() +
  labs(x = "Direct distance [m]", y = "Trip duration [s]")

unnamed-chunk-31-1.png
わかったことは、

  • 距離はtrip_durationの増加とともに増加する。
  • 24時間の移動はさらに疑わしくなり、人為的要因によりミスの可能性がさらに高まりました。
  • 非常に短い距離、1mまでの数の移動があり、明らかなtrip_durationsの範囲が広い。

極端な(そして非常に疑わしい)データポイントを削除するためにデータを少しフィルタリングし、データを2-dヒストグラムのためにビン化しましょう。このプロットは、log-log空間で、より大きな距離の場合、trip_durationが線形よりも鈍く増加していることを示しています。

R
train %>%
  filter(trip_duration < 3600 & trip_duration > 120) %>%
  filter(dist > 100 & dist < 100e3) %>%
  ggplot(aes(dist, trip_duration)) +
  geom_bin2d(bins = c(500,500)) +
  scale_x_log10() +
  scale_y_log10() +
  labs(x = "Direct distance [m]", y = "Trip duration [s]")

unnamed-chunk-32-1.png

4.2 移動速度

時間の経過の距離は当然速度と関係があり、タクシーの平均速度を計算することで、偽の値を取り除くことができます。もちろん、移動時間を知る必要があるため、モデルの予測変数としてスピードを使用することはできませんが、トレーニングデータのクリーンアップや予測能力を備えた他の特徴量の検索に役立ちます。これは速度分布です。
unnamed-chunk-33-1.png
極端な値を取り除きましたが、予想したよりも良く見えます。おそらくNYCでは約15km/hの平均速度と考えられそうです。50km/hを超える車は、魔法かなにか(または高速道路の旅)が必要です。また、これは直接距離を指し、実速度は常に高くなっていたことに留意してください。1日と1時間あたりの平均所要時間と同様に、これらの時間帯の平均速度も調べます。

R
p1 <- train %>%
  group_by(wday, vendor_id) %>%
  summarise(median_speed = median(speed)) %>%
  ggplot(aes(wday, median_speed, color = vendor_id)) +
  geom_point(size = 4) +
  labs(x = "Day of the week", y = "Median speed [km/h]")

p2 <- train %>%
  group_by(hour, vendor_id) %>%
  summarise(median_speed = median(speed)) %>%
  ggplot(aes(hour, median_speed, color = vendor_id)) +
  geom_smooth(method = "loess", span = 1/2) +
  geom_point(size = 4) +
  labs(x = "Hour of the day", y = "Median speed [km/h]") +
  theme(legend.position = "none")

p3 <- train %>%
  group_by(wday, hour) %>%
  summarise(median_speed = median(speed)) %>%
  ggplot(aes(hour, wday, fill = median_speed)) +
  geom_tile() +
  labs(x = "Hour of the day", y = "Day of the week") +
  scale_fill_distiller(palette = "Spectral")

layout <- matrix(c(1,2,3,3),2,2,byrow=TRUE)
multiplot(p1, p2, p3, layout=layout)

unnamed-chunk-34-1.png

R
p1 <- 1; p2 <- 1; p3 <- 1

わかったことは、
- タクシーは、週末と月曜日は他の曜日よりも速く移動しているようです。
- 早朝の時間は、8時から6時までは速く、より迅速な移動を可能にします。
- 2つのベンダーの間にはほとんど違いがありません。
- ヒートマップは、これらの傾向がどのように組み合わされて、1日と1週間の関係性において「低速度ゾーン」を作り出すかを可視化します。これに基づいて、私たちはMon-Fri,8am-6pmとして定義する新しい特徴量を作成します。

4.3 ベアリングの方向

直接距離が移動ベクトルの大きさである場合、ベアリングはその方向です。 geosphereパッケージを通じて容易に推定され、北西または南東の方向に出発したかどうかがわかります。 ここでは、ベアリング分布と、trip_duration、直接距離、および速度との関係を視覚化します。

R
p1 <- train %>%
  filter(dist < 1e5) %>%
  ggplot(aes(bearing)) +
  geom_histogram(fill = "red", bins = 75) +
  scale_x_continuous(breaks = seq(-180, 180, by = 45)) +
  labs(x = "Bearing")

p2 <- train %>%
  filter(dist < 1e5) %>%
  ggplot(aes(bearing, dist)) +
  geom_bin2d(bins = c(100,100)) +
  labs(x = "Bearing", y = "Direct distance") +
  scale_y_log10() +
  theme(legend.position = "none") +
  coord_polar() +
  scale_x_continuous(breaks = seq(-180, 180, by = 45))

p3 <- train %>%
  filter(trip_duration < 3600*22) %>%
  filter(dist < 1e5) %>%
  ggplot(aes(bearing, trip_duration)) +
  geom_bin2d(bins = c(100,100)) +
  scale_y_log10() +
  labs(x = "Bearing", y = "Trip duration") +
  coord_polar() +
  scale_x_continuous(breaks = seq(-180, 180, by = 45))

p4 <- train %>%
  filter(speed < 75 & dist < 1e5) %>%
  ggplot(aes(bearing, speed)) +
  geom_bin2d(bins = c(100,100)) +
  labs(x = "Bearing", y = "Speed") +
  coord_polar() +
  scale_x_continuous(breaks = seq(-180, 180, by = 45))

layout <- matrix(c(1,2,3,4),2,2,byrow=TRUE)
multiplot(p1, p2, p3, p4, layout=layout)

unnamed-chunk-35-1.png
わかったことは、

  • ベアリングの方向は、30度と-150度に顕著なピークがある。直観的には、それらはマンハッタンの方向性に関係するのかもしれません。これらの大きなピークの間には、3つまたは4つの小さくて鋭いピークが見える。
  • distancetrip_durationの2Dヒストグラムは、とくに得られるものはありません。2つのピークの時間と距離が短いことが示されていますが、これらのピークの観測数が多いことが原因です。より高いtrip_durationsに向かう淡いクラスタは、-60度と125度の方位付近に見えるかもしれません。
  • distancetrip_durationで見られる明るい「リング」は、対数スケール上の分布のピークを反映しています。
  • スピードは、ベアリング方向に沿って興味深いクラスター構造を示しています。より速い速度の領域は、前述のベアリング方向のピークにおいても、それらに対しておおよそ垂直に(最も顕著には約125度)識別することができます。このデータにマンハッタン通りのグリッドが表示されている可能性があります。ここで、半径スケールは対数ではないことに注意してください。

4.4 空港までの距離

マップと移動のパスでは、NYCの2つの空港(JFKとLa Guardia)のいずれかで多数の移動が開始または終了したことがわかりました。空港は通常市内中心部にはないので、空港からのピックアップ/ドロップオフ距離がより長いtrip_durationsの有用な予測値になると仮定することは合理的です。上の例では、2つの空港の座標を定義し、対応する距離を計算しました。

R
p1 <- train %>%
  ggplot(aes(jfk_dist_pick)) +
  geom_histogram(bins = 30, fill = "red") +
  scale_x_log10() +
  scale_y_sqrt() +
  geom_vline(xintercept = 2e3) +
  labs(x = "JFK pickup distance")

p2 <- train %>%
  ggplot(aes(jfk_dist_drop)) +
  geom_histogram(bins = 30, fill = "blue") +
  scale_x_log10() +
  scale_y_sqrt() +
  geom_vline(xintercept = 2e3) +
  labs(x = "JFK dropoff distance")

p3 <- train %>%
  ggplot(aes(lg_dist_pick)) +
  geom_histogram(bins = 30, fill = "red") +
  scale_x_log10() +
  scale_y_sqrt() +
  geom_vline(xintercept = 2e3) +
  labs(x = "La Guardia pickup distance")

p4 <- train %>%
  ggplot(aes(lg_dist_drop)) +
  geom_histogram(bins = 30, fill = "blue") +
  scale_x_log10() +
  scale_y_sqrt() +
  geom_vline(xintercept = 2e3) +
  labs(x = "La Guardia dropoff distance")

layout <- matrix(c(1,2,3,4),2,2,byrow=FALSE)
multiplot(p1, p2, p3, p4, layout=layout)

unnamed-chunk-36-1.png

R
p1 <- 1; p2 <- 1; p3 <- 1; p4 <- 1; p5 <- 1

これらの数値に基づいて、JFK/La Guardiaの移動は、空港から2km以内にあるピックアップまたはドロップオフを持つと定義することができます。これらの移動の移動時間はどのくらいでしょうか?

R
p1 <- train %>%
  filter(trip_duration < 23*3600) %>%
  ggplot(aes(jfk_trip, trip_duration, color = jfk_trip)) +
  geom_boxplot() +
  scale_y_log10() +
  theme(legend.position = "none") +
  labs(x = "JFK trip")

p2 <- train %>%
  filter(trip_duration < 23*3600) %>%
  ggplot(aes(lg_trip, trip_duration, color = lg_trip)) +
  geom_boxplot() +
  scale_y_log10() +
  theme(legend.position = "none") +
  labs(x = "La Guardia trip")

layout <- matrix(c(1,2),1,2,byrow=FALSE)
multiplot(p1, p2, layout=layout)

unnamed-chunk-37-1.png
仮説は正しく、空港、特に遠方のJFKへの移動は平均移動時間がかなり長いことがわかりました。 これらの2つの特徴量は間違いなく私たちのモデルの一部になるはずです。

5 データクリーニング

モデリングを始める前に、訓練データを整理する必要があります。これまで、問題のある特徴量の値お全体像を得るために、データクリーニングの実施を遅らせていました。ここでの目的は、極端な移動時間や平均速度が非常に遅いなど、不可解な値を持つ特徴量を削除することです。
テストデータには数多くの偽の移動時間があるかもしれませんが、どのような場合でもそれらを予測できません。訓練データから、これらの値を削除することにより、モデルをより強固にし、目に見えないデータを汎用化する可能性が高くなります。これは、常に機械学習の主な目標でもあります。

5.1 極端な移動時間

1日以上の時間がかかる移動距離を視覚化しましょう。誰かがNYCからLAにタクシーを使っていない限り、その値が正確であるとは思われません。ここではmapsパッケージを使用して、ほとんどの移動が始まるまたは終わるマンハッタンの概要を描きます。次に、ピックアップ座標を赤で表示し、ドロップオフ座標を青で表示します。
直接距離をさらに追加するために、geosphereパッケージから別の便利なツールを使用します。gcIntermediateを使用すると、2組の座標の間のパスを補間することができます。私はJonathan Bouchetのカーネルで、このツールを使っているのをみました。もしあなたが彼のカーネルを見たことがないなら、私は読むことを強くお勧めします。

5.1.1 24時間以上の移動

数日間かかる少数の移動から視覚化を始めます。

R
day_plus_trips <- train %>%
  filter(trip_duration > 24*3600)

day_plus_trips %>% select(pickup_datetime, dropoff_datetime, speed)

## # A tibble: 4 x 3
##   pickup_datetime     dropoff_datetime      speed
##   <dttm>              <dttm>                <dbl>
## 1 2016-01-05 00:19:42 2016-01-27 11:08:38 0.0374 
## 2 2016-02-13 22:38:00 2016-03-08 15:57:38 0.0105 
## 3 2016-01-05 06:14:15 2016-01-31 01:01:07 0.00265
## 4 2016-02-13 22:46:52 2016-03-25 18:18:14 0.0203

ny_map <- as.tibble(map_data("state", region = "new york:manhattan"))

tpick <- day_plus_trips %>%
  select(lon = pickup_longitude, lat = pickup_latitude)

tdrop <- day_plus_trips %>%
  select(lon = dropoff_longitude, lat = dropoff_latitude)

p1 <- ggplot() +
  geom_polygon(data=ny_map, aes(x=long, y=lat), fill = "grey60") +
  geom_point(data=tpick,aes(x=lon,y=lat),size=1,color='red',alpha=1) +
  geom_point(data=tdrop,aes(x=lon,y=lat),size=1,color='blue',alpha=1)

for (i in seq(1,nrow(tpick))){
  inter <- as.tibble(gcIntermediate(tpick[i,],  tdrop[i,], n=30, addStartEnd=TRUE))
  p1 <- p1 +  geom_line(data=inter,aes(x=lon,y=lat),color='blue',alpha=.75)
}

p1 + ggtitle("Longer than a day trips in relation to Manhattan")

unnamed-chunk-39-1.png
ここでは問題以外のものは出てきません。飛行機にすぐに搭乗する場合、JFKへの移動は永遠のように見えるかもしれませんが、リアルタイムでこんなに長くなることはまずありません。平均タクシースピードはあまり見えません。

決定:継続的な探索とモデリングのために、これらの値を訓練データセットから削除する必要があります。

5.1.2 24時間に近い移動

誰かがタクシーに乗って、ほぼ一日乗り続けることは想像もつきません。ごくまれに、これが起こるかもしれませんが、もちろん、移動距離は十分に長いものである場合です。
ここでは、1日の旅行を22時間から24時間とすることを定義します。これは、生のtrip_durationの分布の小さなピークをカバーします。これらは、1日の移動の中でトップ5の直接距離(m)です。

R
day_trips <- train %>%
  filter(trip_duration < 24*3600 & trip_duration > 22*3600)

day_trips %>% 
  arrange(desc(dist)) %>%
  select(dist, pickup_datetime, dropoff_datetime, speed) %>%
  head(5)

## # A tibble: 5 x 4
##    dist pickup_datetime     dropoff_datetime    speed
##   <dbl> <dttm>              <dttm>              <dbl>
## 1 60666 2016-06-04 13:54:29 2016-06-05 13:40:30 2.55 
## 2 42401 2016-01-28 21:43:02 2016-01-29 21:33:30 1.78 
## 3 28116 2016-03-14 22:46:05 2016-03-15 22:24:27 1.19 
## 4 22808 2016-06-18 17:41:47 2016-06-19 16:30:41 1.000
## 5 22759 2016-06-01 19:52:42 2016-06-02 19:35:09 0.960

一番上のものは約60km(約37マイル)です。 平均スピードは、いいヒントを示唆しません。これらの移動は地図上でどのように見えるのでしょうか?

R
ny_map <- as.tibble(map_data("state", region = "new york:manhattan"))

set.seed(2017)
day_trips <- day_trips %>%
  sample_n(200)
tpick <- day_trips %>%
  select(lon = pickup_longitude, lat = pickup_latitude)

tdrop <- day_trips %>%
  select(lon = dropoff_longitude, lat = dropoff_latitude)

p1 <- ggplot() +
  geom_polygon(data=ny_map, aes(x=long, y=lat), fill = "grey60") +
  geom_point(data=tpick,aes(x=lon,y=lat),size=1,color='red',alpha=1) +
  geom_point(data=tdrop,aes(x=lon,y=lat),size=1,color='blue',alpha=1)

for (i in seq(1,nrow(tpick))){
  inter <- as.tibble(gcIntermediate(tpick[i,],  tdrop[i,], n=30, addStartEnd=TRUE))
  p1 <- p1 +  geom_line(data=inter,aes(x=lon,y=lat),color='blue',alpha=.25)
}

p1 + ggtitle("Day-long trips in relation to Manhattan")

unnamed-chunk-41-1.png
ここでは地図を読みやすくし、スクリプトを高速に処理するために、約1800のうち200個だけをプロットしています。ピックアップポイントは赤で、ドロップオフポイントは青です。

わかったこと:

  • いくつかの距離が目立ちますが、例外と思われます。2つの主要な移動グループは、マンハッタン内のものと、マンハッタンと空港間のものです。
  • 極端なtrip_durationsが実際のものであることを示唆するものはほとんどありません。
  • 直感的なもう一つの洞察があります。いずれの空港(最も顕著なJFK)への出入りは非常に短いとは言えません。したがって、ピックアップまたは空港へのドロップオフのいずれかの距離は、より長いtrip_durationの有益な予測値となり得ます。

決定:探索的データ分析の結果、モデリングから22時間以上のtrip_durationsを削除します。

5.1.3 数分の移動

trip_durationの分布には、数分しか移動してい記録があります。そのような短時間の移動は可能ですが、現実的であることを確認するために、その時間と速度を確認しましょう。

R
min_trips <- train %>%
  filter(trip_duration < 5*60)

min_trips %>% 
  arrange(dist) %>%
  select(dist, pickup_datetime, dropoff_datetime, speed) %>%
  head(5)

## # A tibble: 5 x 4
##    dist pickup_datetime     dropoff_datetime    speed
##   <dbl> <dttm>              <dttm>              <dbl>
## 1     0 2016-02-29 18:39:12 2016-02-29 18:42:59     0
## 2     0 2016-01-27 22:29:31 2016-01-27 22:29:58     0
## 3     0 2016-01-22 16:13:01 2016-01-22 16:13:20     0
## 4     0 2016-01-18 15:24:43 2016-01-18 15:28:57     0
## 5     0 2016-05-04 22:28:43 2016-05-04 22:32:51     0
5.1.3.1 距離が0

比較的多くの0距離移動があることがわかります。

R
zero_dist <- train %>%
  filter(near(dist,0))

nrow(zero_dist)
## [1] 5897

名目上のトップ移動時間は何でしょうか。

R
What are their nominal top durations?

zero_dist %>%
  arrange(desc(trip_duration)) %>%
  select(trip_duration, pickup_datetime, dropoff_datetime, vendor_id) %>%
  head(5)

## # A tibble: 5 x 4
##   trip_duration pickup_datetime     dropoff_datetime    vendor_id
##           <int> <dttm>              <dttm>              <fct>    
## 1         86352 2016-06-05 01:09:39 2016-06-06 01:08:51 2        
## 2         85333 2016-01-01 15:27:28 2016-01-02 15:09:41 2        
## 3         78288 2016-05-18 13:40:45 2016-05-19 11:25:33 2        
## 4          5929 2016-06-08 16:47:44 2016-06-08 18:26:33 2        
## 5          4683 2016-05-25 17:36:49 2016-05-25 18:54:52 2

実際には、約1日移動していないと思うタクシーがいくつかあります。この場合のデータを信じないことにします。極端なケースを取り除くと、次のようになります。

R
zero_dist %>%
  filter(trip_duration < 6000) %>%
  ggplot(aes(trip_duration, fill = vendor_id)) +
  geom_histogram(bins = 50) +
  scale_x_log10()

unnamed-chunk-45-1.png
わかったことは、

  • 誰かがタクシーに乗ったが、タクシーが動く前に考えが変わったと仮定すると、約1分のtrip_durationsは可能かもしれません。それが「移動」として数えられるかどうかは、別の問題です。しかし、約15分(900秒)の距離をカバーすることのない移動時間はほとんど不可能です。渋滞か、誰が渋滞に止まってタクシーに入るのでなければ…
  • 1分未満のグループのほとんどの移動は、ベンダー1からのものでしたが、10分のグループは主にベンダー2のタクシーで構成されていました。

決定:分析のために1分以上かかる0距離移動を削除します。モデリングからそれらを削除することは、テストデータに同様の移動がある場合には有害である可能性がある。

5.1.3.2 距離が0で短い移動

これっらは、0距離移動を除去した後の平均速度が速いの短い移動時間です。

R
min_trips <- train %>%
  filter(trip_duration < 5*60 & dist > 0)

min_trips %>% 
  arrange(desc(speed)) %>%
  select(trip_duration, dist, pickup_datetime, speed) %>%
  head(10)

## # A tibble: 10 x 4
##    trip_duration   dist pickup_datetime     speed
##            <int>  <dbl> <dttm>              <dbl>
##  1             7  18055 2016-02-13 20:28:30  9285
##  2           282 320484 2016-04-02 20:33:19  4091
##  3             2    784 2016-06-12 06:35:13  1412
##  4            51  19970 2016-06-19 20:18:35  1410
##  5           279 104877 2016-03-20 21:07:56  1353
##  6             2    704 2016-06-23 13:36:48  1267
##  7            20   6919 2016-05-28 15:14:19  1245
##  8             3    914 2016-05-30 17:12:12  1097
##  9             4   1031 2016-01-17 03:11:56   928
## 10             3    762 2016-06-13 16:34:30   914

明らかに、最高速度の値は不可能なものです。サンプル内に一般的な距離分布のマッププロットが含まれています。

R
ny_map <- as.tibble(map_data("state", region = "new york:manhattan"))

set.seed(1234)
foo <- min_trips %>%
  sample_n(600)

tpick <- foo %>%
  select(lon = pickup_longitude, lat = pickup_latitude)
tdrop <- foo %>%
  select(lon = dropoff_longitude, lat = dropoff_latitude)

p1 <- ggplot() +
  geom_polygon(data=ny_map, aes(x=long, y=lat), fill = "grey60") +
  geom_point(data=tpick,aes(x=lon,y=lat),size=1,color='red',alpha=1) +
  geom_point(data=tdrop,aes(x=lon,y=lat),size=1,color='blue',alpha=1)

for (i in seq(1,nrow(tpick))){
  inter <- as.tibble(gcIntermediate(tpick[i,],  tdrop[i,], n=30, addStartEnd=TRUE))
  p1 <- p1 +  geom_line(data=inter,aes(x=lon,y=lat),color='blue',alpha=.25)
}

p1 + ggtitle("Minute-long trips in relation to Manhattan")

unnamed-chunk-47-1.png

5.2 休憩?最も怪しい移動

すべてのデータセットには、ばかばかしいエントリがいくつかあります。ここには、JFK空港から300キロ以上離れた場所のピックアップまたはドロップオフの場所があります。

R
long_dist <- train %>%
  filter( (jfk_dist_pick > 3e5) | (jfk_dist_drop > 3e5) )

long_dist_coord <- long_dist %>%
  select(lon = pickup_longitude, lat = pickup_latitude)

long_dist %>%
  select(id, jfk_dist_pick, jfk_dist_drop, dist, trip_duration, speed) %>%
  arrange(desc(jfk_dist_pick))

## # A tibble: 31 x 6
##    id        jfk_dist_pick jfk_dist_drop      dist trip_duration     speed
##    <chr>             <dbl>         <dbl>     <dbl>         <int>     <dbl>
##  1 id2854272       4128727       4128719      14.8           499    0.107 
##  2 id3777240       4128721       4128726      21.8          1105    0.0711
##  3 id2306955       1253574         21484 1242299             792 5647     
##  4 id1974018       1115646       1115646       0             369    0     
##  5 id0267429        988748        988748       0             961    0     
##  6 id0838705        695767        481141  215468            1131  686     
##  7 id0205460        684133        684133       0             329    0     
##  8 id0978162        674251        941618  315117             875 1296     
##  9 id3525158        665877        665877       0             385    0     
## 10 id1510552        642663        472020  892212             611 5257    

とにかく削除しますが、1分以上の長時間の0距離移動は、どこで起こったものでしょうか?

R
leaflet(long_dist_coord) %>%
  addTiles() %>%
  setView(-92.00, 41.0, zoom = 4) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addMarkers(popup = ~as.character(long_dist$dist), label = ~as.character(long_dist$id))

スクリーンショット 0030-09-08 3.25.50.png
2つの孤独なNYCのタクシーはサンフランシスコの近くや、海に行くタクシーもあります。これらの場所は、予測モデルの汎化性を向上させるために除去すべき異常値です。

5.3 最終的なデータクリーニング

ここでは、前述のクリーニングするデータにフィルタを適用します。このコードブロックは、分析が進むにつれて拡大する可能性があります。

R
train <- train %>%
  filter(trip_duration < 22*3600,
         dist > 0 | (near(dist, 0) & trip_duration < 60),
         jfk_dist_pick < 3e5 & jfk_dist_drop < 3e5,
         trip_duration > 10,
         speed < 100)

6 外部データ

ここでは、分析に必要な追加のデータソースを補足することが奨励されています。これらはKaggleデータセットとして公開され、このディスカッションスレッドで収集されています。Kaggleカーネルに複数のデータソースを含める方法についてもっと知りたければ、この記事をチェックしてください。

6.1 天気予報

まず、Mathijs Waegemakersが提供するNYC気象データのリストをここに組み込むことから始めます。このデータの内容は下記の通りです。

国立気象サービスから収集された天気データ。それは中央公園の気象観測所の2016年から半年をカバーしています。毎日の最低温度、最高温度、平均温度、降水量、新しい降雪量、および現在の積雪深度を含んでいます。温度は華氏で測定され、深さはインチで測定されます。Tは降水量があることを意味します。

ここで特に興味深いのは、雨と雪の秋の統計です。これは1月下旬の移動数の減少と比較することに関心があります。同じデータセットの独立した分析を示すこのカーネルも確認してください。

6.1.1 データのインポートからジョインまで

R
weather <- as.tibble(fread("../input/weather-data-in-new-york-city-2016/weather_data_nyc_centralpark_2016(1).csv"))
glimpse(weather)

## Observations: 366
## Variables: 7
## $ date                  <chr> "1-1-2016", "2-1-2016", "3-1-2016", "4-1...
## $ `maximum temperature` <int> 42, 40, 45, 36, 29, 41, 46, 46, 47, 59, ...
## $ `minimum temperature` <int> 34, 32, 35, 14, 11, 25, 31, 31, 40, 40, ...
## $ `average temperature` <dbl> 38.0, 36.0, 40.0, 25.0, 20.0, 33.0, 38.5...
## $ precipitation         <chr> "0.00", "0.00", "0.00", "0.00", "0.00", ...
## $ `snow fall`           <chr> "0.0", "0.0", "0.0", "0.0", "0.0", "0.0"...
## $ `snow depth`          <chr> "0", "0", "0", "0", "0", "0", "0", "0", ...

日付をrubridateオブジェクトに変換し、雨と雪のデータ(="T")を、小さな数値に変換します。 また、最高温度と最低温度をより短い形式で保存します。

R
weather <- weather %>%
  mutate(date = dmy(date),
         rain = as.numeric(ifelse(precipitation == "T", "0.01", precipitation)),
         s_fall = as.numeric(ifelse(`snow fall` == "T", "0.01", `snow fall`)),
         s_depth = as.numeric(ifelse(`snow depth` == "T", "0.01", `snow depth`)),
         all_precip = s_fall + rain,
         has_snow = (s_fall > 0) | (s_depth > 0),
         has_rain = rain > 0,
         max_temp = `maximum temperature`,
         min_temp = `minimum temperature`)

データを訓練データに結合します。

R
foo <- weather %>%
  select(date, rain, s_fall, all_precip, has_snow, has_rain, s_depth, max_temp, min_temp)

train <- left_join(train, foo, by = "date")

6.1.2 視覚化とtrip_durationへの影響

降雪の統計と1日あたりの移動数を比較してみましょう。

R
p1 <- train %>%
  group_by(date) %>%
  count() %>%
  ggplot(aes(date,n/1e3)) +
  geom_line(size = 1.5, color = "red") +
  labs(x = "", y = "Kilo trips per day")

p2 <- train %>%
  group_by(date) %>%
  summarise(trips = n(),
            snow_fall = mean(s_fall),
            rain_fall = mean(rain),
            all_precip = mean(all_precip)) %>%
  ggplot(aes(date, snow_fall)) +
  geom_line(color = "blue", size = 1.5) +
  labs(x = "", y = "Snowfall") +
  scale_y_sqrt() +
  scale_x_date(limits = ymd(c("2015-12-28", "2016-06-30")))

p3 <- train %>%
  group_by(date) %>%
  summarise(trips = n(),
            snow_depth = mean(s_depth)) %>%
  ggplot(aes(date, snow_depth)) +
  geom_line(color = "purple", size = 1.5) +
  labs(x = "", y = "Snow depth") +
  scale_y_sqrt() +
  scale_x_date(limits = ymd(c("2015-12-29", "2016-06-30")))

p4 <- train %>%
  group_by(date) %>%
  summarise(median_speed = median(speed)) %>%
  ggplot(aes(date, median_speed)) +
  geom_line(color = "orange", size = 1.5) +
  labs(x = "Date", y = "Median speed")

layout <- matrix(c(1,2,3,4),4,1,byrow=FALSE)
multiplot(p1, p2, p3, p4, layout=layout)

unnamed-chunk-55-1.png

R
p1 <- 1; p2 <- 1; p3 <- 1; p4 <- 1

私は知っていました!移動量の減少は、NYC(1月23日)の冬の最大の(そして最初の)雪の降雪に相当します。実際、NYCは吹雪に襲われ、記録的な降雪を経験しました。交通パターンへの影響は甚大であり、速度の中央値は著しく低下しました。(2つ目のプロットのy軸は平方根です。注意してください。)
この大きなスパイクの後、降雪の次の期間は、タクシー移動数またはその速度にはほとんど影響しません。3月中旬に最後の雪が降った際には、しばらく雪が降っていなかったことが原因で、びっくりしたからかもしれません。
今、雪の降った日の移動は、長い時間の移動につながりますか?この質問に答えるために、私たちは、サンプルの全日中の平均移動時間と全降水量(雨+雪)の散布図を見ます。

R
train %>%
  group_by(date, has_snow) %>%
  summarise(duration = mean(trip_duration),
            all_precip = mean(all_precip)) %>%
  ggplot(aes(all_precip, duration, color = has_snow)) +
  geom_jitter(width = 0.04, size = 2) +
  scale_x_sqrt() +
  scale_y_log10() +
  labs(x = "Amount of total precipitation", y = "Average trip duration")

unnamed-chunk-56-1.png
雪の日を除くと、移動がより短くなるように見えます。しかし、これは単純に、乗客がより短い距離を走行する可能性が高いことを意味する可能性があるので、平均速度はどのような関係にあるのでじょうか。

R
p1 <- train %>%
  filter(speed < 50) %>%
  ggplot(aes(has_snow, speed, color = has_snow)) +
  geom_boxplot() +
  theme(legend.position = "none") +
  labs(x = "Snowfall")

p2 <- train %>%
  filter(speed < 50) %>%
  ggplot(aes(has_rain, speed, color = has_rain)) +
  geom_boxplot() +
  theme(legend.position = "none") +
  labs(x = "Rainfall")

layout <- matrix(c(1,2),1,2,byrow=FALSE)
multiplot(p1, p2, layout=layout)

unnamed-chunk-57-1.png

R
p1 <- 1; p2 <- 1

大きな違いはないことがわかります。1月23日(ブリザード)は、最も長いtrip_durationの中央値を持つ雪の日のリストのトップにはいません(ただし、トップ10に入っています)。

R
train %>%
  filter(has_snow == TRUE) %>%
  group_by(date) %>%
  summarise(med_duration = median(trip_duration),
            med_speed = median(speed)) %>%
  arrange(desc(med_duration)) %>%
  head(10)

## # A tibble: 10 x 3
##    date       med_duration med_speed
##    <date>            <dbl>     <dbl>
##  1 2016-01-26          824      10.6
##  2 2016-01-25          805      10.8
##  3 2016-01-27          752      10.9
##  4 2016-01-28          732      11.5
##  5 2016-01-29          705      11.4
##  6 2016-02-11          670      11.7
##  7 2016-03-04          663      12.0
##  8 2016-02-23          662      11.4
##  9 2016-01-23          660      12.4
## 10 2016-01-19          656      12.1

この表では、最も速度が遅い雪の日のトップ5を見ているとき、興味深いパターンを見つけることができます。興味深いことに、ブリザードの2日後(1月25日から27日まで)に全体的に最も低い速度の中央値が観察されています。

R
train %>%
  group_by(date) %>%
  summarise(med_duration = median(trip_duration),
            med_speed = median(speed)) %>%
  arrange(med_speed) %>%
  head(5)

## # A tibble: 5 x 3
##   date       med_duration med_speed
##   <date>            <dbl>     <dbl>
## 1 2016-01-26          824      10.6
## 2 2016-01-25          805      10.8
## 3 2016-01-27          752      10.9
## 4 2016-06-08          768      10.9
## 5 2016-06-01          761      11.2

結論:探索的な観点からは、これはあまり有望ではありません。それでも、雨や雪のデータを暫定的なモデルに入れて、現時点では見ることのできない相互作用の項があるかどうかを調べることもできます。大雪が降った後の(作業中の)週の「ブリザード効果」は、考慮に入れる価値があります。 したがって、特徴量エンジニアリングのコードブロックにブリザードの特徴量を作成します。

6.2 最短距離

別の興味深い外部データは、OSRMを使用したoscarleoによって提供される各旅行の最速ルートの推定データです。データはここで見つけることができ、ピックアップ/ドロップオフ通りと、これら2つのポイント間の合計距離/継続時間と、ターンやハイウェイに入る一連の移動ステップが含まれます。
ディスカッション項目によると、それらは実際には2つのポイント間の最短ではなく最速のルートです。これは交通量を考慮していない可能性が最も高く、1日で最も速いルートが別の曜日に最も速くない可能性があります。 それでも、ここでは私たちの予測目標に価値があるはずの実際の走行距離と所要時間のデータがあります。

6.2.1 データのインポートと概要

データは3つの別々のファイルに分かれ、訓練データ(2ファイル)とテストデータ用のIDに分かれています。ここでは、トレーニングデータに対応するデータを利用します。このデータセットは、元の訓練データの2倍以上のデータです。

R
foo <- as.tibble(fread("../input/new-york-city-taxi-with-osrm/fastest_routes_train_part_1.csv"))
bar <- as.tibble(fread("../input/new-york-city-taxi-with-osrm/fastest_routes_train_part_2.csv"))
foobar <- as.tibble(fread("../input/new-york-city-taxi-with-osrm/fastest_routes_test.csv"))
fastest_route <- bind_rows(foo, bar, foobar)
glimpse(fastest_route)

## Observations: 2,083,777
## Variables: 12
## $ id                   <chr> "id2875421", "id2377394", "id3504673", "i...
## $ starting_street      <chr> "Columbus Circle", "2nd Avenue", "Greenwi...
## $ end_street           <chr> "East 65th Street", "Washington Square We...
## $ total_distance       <dbl> 2009.1, 2513.2, 1779.4, 1614.9, 1393.5, 1...
## $ total_travel_time    <dbl> 164.9, 332.0, 235.8, 140.1, 189.4, 138.8,...
## $ number_of_steps      <int> 5, 6, 4, 5, 5, 5, 2, 13, 6, 9, 4, 4, 10, ...
## $ street_for_each_step <chr> "Columbus Circle|Central Park West|65th S...
## $ distance_per_step    <chr> "0|576.4|885.6|547.1|0", "877.3|836.5|496...
## $ travel_time_per_step <chr> "0|61.1|60.1|43.7|0", "111.7|109|69.9|25....
## $ step_maneuvers       <chr> "depart|rotary|turn|new name|arrive", "de...
## $ step_direction       <chr> "left|straight|right|straight|arrive", "n...
## $ step_location_list   <chr> "-73.982316,40.767869|-73.981997,40.76768...

glimpseの結果を見ると、各ステップごとに、distance_per_steptravel_time_per_stepおよび移動を説明するstep_maneuversstep_directionの詳細な説明と共に名前が付けられます(streets_for_each_step)。 step_maneuversの名前はわかりやすく、データセットの概要にも記載されています。出発と到着は個別のステップです。step_maneuverの"turn"の方向は、step_directionsにリストされています。総距離は、元のデータと同様に、メートル単位で、秒単位で測定されます。

6.2.2 新しいデータの可視化

「最速ルート」のデータをトレーニングデータに結合させる前に、まず数値の分布を視覚化しましょう。

R
p1 <- fastest_route %>%
  ggplot(aes(total_travel_time)) +
  geom_histogram(bins = 150, fill = "red") +
  scale_x_log10() +
  scale_y_sqrt()

p2 <- fastest_route %>%
  ggplot(aes(total_distance)) +
  geom_histogram(bins = 150, fill = "red") +
  scale_x_log10() +
  scale_y_sqrt()

p3 <- fastest_route %>%
  ggplot(aes(number_of_steps)) +
  geom_bar(fill = "red")

p4 <- fastest_route %>%
  ggplot(aes(total_distance, total_travel_time)) +
  geom_bin2d(bins = c(80,80))

layout <- matrix(c(1,2,3,4,4,4),2,3,byrow=TRUE)
multiplot(p1, p2, p3, p4, layout=layout)

unnamed-chunk-62-1.png
わかったことは、

  • total_travel_timeおよびtotal_distanceは、それぞれ約300秒(5分)および2800mの値でピークに達します。これらは中央値です。どちらの特徴量でも高い値はまれですが、2次ピークのサインがあります。 このデータでも非常に短い距離(数メートル)と移動時間(数秒)が短いデータが存在する。
  • 最も頻繁な{number_of_steps`は5です。「depart」と「arrive」は常に存在するため、2ステップ未満の移動はありません。 これらの2つのステップしか含まれていないほとんど60kの旅行があります。これは、直線で移動することを意味します。
R
fastest_route %>% filter(number_of_steps == 2) %>% nrow()
## [1] 83340
  • より多くのステップが存在するが、約20%だけが10を超え、2%未満が20ステップを超える。
  • total_distancetotal_travel_timeは強く相関しています。 おそらく非常に短い距離を除いて。 この線形の周りにはある程度のばらつきがあり、距離と時間の全範囲にわたってかなり均質に見えます。
R
fastest_route %>%
  ggplot(aes(factor(number_of_steps), total_travel_time, color = factor(number_of_steps))) +
  geom_boxplot() +
  labs(x = "number_of_steps") +
  theme(legend.position = "none")

unnamed-chunk-64-1.png
移動時間の中央値はnumber_of_stepsとよく相関していますが、多くの外れ値があるようです。number_of_stepsごとのtotal_distanceの分布は、上記のプロットと非常によく似ています。

6.2.3 結合と結合後の特徴量

結合されたデータを見てみると、元の訓練データセットと組み合わされたtotal_distancetotal_travel_timeに注目します。

R
foo <- fastest_route %>%
  select(id, total_distance, total_travel_time, number_of_steps,
         step_direction, step_maneuvers) %>%
  mutate(fastest_speed = total_distance/total_travel_time*3.6,
         left_turns = str_count(step_direction, "left"),
         right_turns = str_count(step_direction, "right"),
         turns = str_count(step_maneuvers, "turn")
         ) %>%
  select(-step_direction, -step_maneuvers)

train <- left_join(train, foo, by = "id") %>%
  mutate(fast_speed_trip = total_distance/trip_duration*3.6)

さらに、次の新しい特徴量を作成します。

  • fastest_speedは、OSRM値に基づく最速ルートの平均速度であり、fast_speed_tripは、最速ルートと実際の時間を想定した速度です。 両方の特徴量はkm/hの単位です。
  • left_turnsright_turns、およびturnsは、ルートごとのそれらの行動の合計数です。 step_directionsの "left"と "right"は必ずstep_maneuverの"turn"に関連するとは限りませんが、ほとんどの場合、left turnとright turnの最初の推定値として使用します。特に左折すると時間が遅くなることがあります。

それ以外に、使用された道路の詳細なリストやそれぞれのtravel_time_per_stepは、交通の速度が遅い領域を考慮に入れられるかもしれません。

6.2.4 視覚化とtrip_durationへの影響

R
p1 <- train %>%
  ggplot(aes(trip_duration)) +
  geom_density(fill = "red", alpha = 0.5) +
  geom_density(aes(total_travel_time), fill = "blue", alpha = 0.5) +
  scale_x_log10() +
  coord_cartesian(xlim = c(5e1, 8e3))

p2 <- train %>%
  ggplot(aes(dist)) +
  geom_density(fill = "red", alpha = 0.5) +
  geom_density(aes(total_distance), fill = "blue", alpha = 0.5) +
  scale_x_log10() +
  coord_cartesian(xlim = c(2e2, 5e4))

p3 <- train %>%
  ggplot(aes(trip_duration, total_travel_time)) +
  geom_bin2d(bins = c(120,120)) +
  geom_abline(intercept = 0, slope = 1, colour = "red") +
  theme(legend.position = "none")

p4 <- train %>%
  ggplot(aes(dist, total_distance)) +
  geom_bin2d(bins = c(70,70)) +
  geom_abline(intercept = 0, slope = 1, colour = "red") +
  theme(legend.position = "none")

layout <- matrix(c(1,2,3,4),2,2,byrow=TRUE)
multiplot(p1, p2, p3, p4, layout=layout)

unnamed-chunk-66-1.png
ここでは最速の特徴量は青色で、元のtrip_durationと直接距離は赤です。ビニングされた散布図では、赤い線は勾配が1を意味します。

わかったことは、

  • 実際のtrip_durationsは、最速ルートの平均よりもかなり平均的に長くなります。分布の形状も微妙に異なります。最も速い分布では、対数ベースですが、左の歪度よりも右の歪度が大きくなります。しかしながら、分布の幅は広く似ています。
  • trip_durationtotal_travel_time(左下隅)のビニングされた散布図は、trip_durationsに向かって多くの外れ値があることを示しています。興味深いことに、それらのほとんどすべてがtotal_travel_timeの2000秒(約0.55時間)以下に現れます。 これらについては、以下で詳しく説明します。
  • 直接距離は最速の距離よりも長いが、それほどではありません。全体として、直接距離は、ビン分割された散布図に従って、実際の最速の距離に対する良い代理変数です。距離が長くなると、潜在的に問題のある外れ値があります。
  • スピードはどうでしょうか?元の速度推定値をOSRMデータからの最速速度値と比較すると、元の特徴量の値がどれほど現実的であるかがわかります。 この図では、カーネル密度プロットの最速値は青で、ビン分割された散布図の赤線は勾配が1になります。
R
p1 <- train %>%
  ggplot(aes(speed)) +
  geom_density(fill = "red", alpha = 0.5) +
  geom_density(aes(fastest_speed), fill = "blue", alpha = 0.5) +
  coord_cartesian(xlim = c(0, 80))

p2 <- train %>%
  ggplot(aes(speed, fastest_speed)) +
  geom_bin2d(bins = c(50,50)) +
  geom_abline(intercept = 0, slope = 1, colour = "red") +
  theme(legend.position = "none")

layout <- matrix(c(1,2),1,2,byrow=TRUE)
multiplot(p1, p2, layout=layout)

unnamed-chunk-67-1.png
ここではいくつかの効果を見つけました。

  • 交通量を考慮しないため、最速の速度は常に20km/h(12mph)を上回り、実際のデータではそれほど頻繁にはありません。
  • 理論的な速度分布はいくつかのピークを示しており、最も顕著には25〜30km/h(15〜18mph)の間で非常に鋭いピークを示し、40km/h(25mph)付近でより広いピークを示します。NYCの交通管理についてよく知らないと、それらが異なるスピード制限区域と線形の組み合わせに対応していると考えてしまいます。
  • speedfastest_speedの間に相関はありません。これは、trip_durationtotal_travel_timeの散布図は、多くの遅延を反映しており、これは別で説明します。

total_travel_timeに対するターン数の影響は、left_turnsright_turnsの場合と非常によく似ています。 これの背後にある主な理由は、total_travel_timenumber_of_stepsの増加に伴って増加し、ターン数がもちろんnumber_of_stepsと相関します。私たちは1回のターン当たりのターンの割合からより多くを学べます。これらのプロットから、次の結果が得られます。

R
p1 <- train %>%
  ggplot(aes(factor(left_turns), total_travel_time, color = factor(left_turns))) +
  geom_boxplot() +
  labs(x = "Number of left turns") +
  theme(legend.position = "none")

p2 <- train %>%
  ggplot(aes(factor(right_turns), total_travel_time, color = factor(right_turns))) +
  geom_boxplot() +
  labs(x = "Number of right turns") +
  theme(legend.position = "none")

p3 <- train %>%
  ggplot(aes(factor(turns), total_travel_time, color = factor(turns))) +
  geom_boxplot() +
  labs(x = "Total number of turns") +
  theme(legend.position = "none")

layout <- matrix(c(1,2,3,3),2,2,byrow=TRUE)
multiplot(p1, p2, p3, layout=layout)

unnamed-chunk-68-1.png

  • ここでは、fastest_routesのデータから欠落している単一の観測(「id3008062」)によって引き起こされる「NA」レコードがあります。これは私たちの予測に実用的な影響を与えることなく、ごくわずかな問題ですが、この「NA」エントリに気づくのは良いことです。OSRMはこのデータポイントのルートを提供していないことがわかります。
  • 右折は、左折よりもわずか多いです(最大21対19)。
  • トータルのターン数は、total_travel_timeの13ターンで最初の最大値があり、その後、中央値は減少していきます。

6.2.5 trip_duration vs total_travel_time-興味深い遅延
このセクションでは、元のデータのOSRMのtotal_travel_timetaxi trip_durationを比較します。特に、データの大幅な遅延を示唆するケースを理解することに重点を置いています。ここでは、OSRMルートの左折のパーセンテージで色分けされた散布図をプロットします(赤の実線は再び傾き1)。

R
train %>%
  ggplot(aes(trip_duration, total_travel_time, color = left_turns / number_of_steps)) +
  geom_point(alpha = 0.3) +
  geom_abline(intercept = 0, slope = 1, colour = "red") +
  scale_colour_gradient(low = "blue", high = "red") +
  scale_x_log10() +
  scale_y_log10() +
  geom_hline(yintercept = c(50,2000), linetype = 2) +
  labs(color = "Left turns [%]")

unnamed-chunk-70-1.png
わかったことは、

  • trip_duration(タクシー移動の実際のタイミング)が最も速いOSRMルートに基づくtotal_travel_timeよりも低いところにデータポイントの集まりがあります。これらの差は、最も速いルートでは1000秒未満(すなわち約17分間)かかるが、それよりも長い乗り物ではそれほど重要にならない移動にとっては有意義であり得る。
  • total_travel_timeが10以上の「最速」ルートはありません。データクリーニングしたデータの最大値は、約5000sまたは1.5hである。trip_duration(水平の破線で示されている)の重要な外れ値のほとんどは、total_travel_timeが100〜2000sの範囲にあり、 trip_durationが10ksを超えています。
  • より長い時間の移動には、より多くのleft_turnsを含める傾向があるようです。しかしながら、この明らかな効果は、旋回を含まない最短の乗り物と、平均旋回数を含むより長い乗り物(1分を超える)に起因する可能性があります。いずれにしても、外れ値の中には、left_turnsのパーセンテージの明らかな影響はなさそうです。

最も顕著な外れ値は全体のごくわずかしか占めません。

R
train %>%
  filter(total_travel_time > 50 & total_travel_time < 2000) %>%
  mutate(delay = trip_duration > 1e4) %>%
  group_by(delay) %>%
  count()

## # A tibble: 2 x 2
## # Groups:   delay [2]
##   delay       n
##   <lgl>   <int>
## 1 F     1434779
## 2 T         207

バープロットと密度オーバーレイの相対頻度の比較を使用して、これらのデータを調べます。

R
foo <- train %>%
  filter(total_travel_time > 50 & total_travel_time < 2000) %>%
  mutate(delay = trip_duration > 1e4)

p1 <- foo %>%
  ggplot(aes(vendor_id, group = delay)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") + 
  scale_y_continuous(labels=scales::percent) +
  ylab("relative frequencies") +
  facet_grid(~delay) +
  theme(legend.position = "none")

p2 <- foo %>%
  ggplot(aes(jfk_trip, group = delay)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") + 
  scale_y_continuous(labels=scales::percent) +
  ylab("relative frequencies") +
  facet_grid(~delay) +
  theme(legend.position = "none")

p3 <- foo %>%
  ggplot(aes(lg_trip, group = delay)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") + 
  scale_y_continuous(labels=scales::percent) +
  ylab("relative frequencies") +
  facet_grid(~delay) +
  theme(legend.position = "none")

p4 <- foo %>%
  ggplot(aes(work, group = delay)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") + 
  scale_y_continuous(labels=scales::percent) +
  ylab("relative frequencies") +
  facet_grid(~delay) +
  theme(legend.position = "none")

p5 <- foo %>%
  ggplot(aes(total_distance, fill = delay)) +
  geom_density(alpha = 0.5)

layout <- matrix(c(1,2,3,4,5,5),2,3,byrow=TRUE)
multiplot(p1, p2, p3, p4, p5, layout=layout)

unnamed-chunk-72-1.png
わかったことは、

  • ベンダー2は乗り物全体の割合がわずかに高いにもかかわらず、ベンダー2のタクシーではほとんどすべてで長い遅延が発生します。この不一致は、効果の大部分を説明するのに間違いなく重要です。
  • もう一つの興味深い事実は、長く遅延するサンプルでは、​​JFK旅行の割合がかなり高く、ラ・ガーディア旅行の割合は一般の人口よりわずかに高いということです。空港の旅ではtrip_durationには乗客の待ち時間が含まれている可能性があるという仮説を立てることができます。JFKの旅行は、長期間の旅行の約30%を占めます(207のうち60)。
  • 勤務時間中の遅延の割合はわずかに高いが、実際的な意味はないです。
  • さらに、長く遅延する分布は、長距離は、短い距離よりも明らかに好意的です。ここではtotal_travel_time(50〜2000s)の同じ括弧内の比較だけを比較しています。

trip_duration vs total_travel_timeのプロットの中央下部には、異常値のグループがあります。ここではベンダー2のベンダーIDによる旅行の色をエンコードし、右上の長期間のアウトライヤーと右下の数々のアウトライヤーへの影響をさらに強調します。

R
train %>%
  ggplot(aes(trip_duration, total_travel_time, color = vendor_id)) +
  geom_point(alpha = 0.3) +
  geom_abline(intercept = 0, slope = 1, colour = "red") +
  scale_x_log10() +
  scale_y_log10() +
  geom_hline(yintercept = c(10), linetype = 2) +
  geom_vline(xintercept = c(600,5000), linetype = 2)

外れ値の第2のグループは、total_travel_time <10かつtrip_duration = 600~4000です。 通常どおり秒単位で表示されます。 隠された図または図26を参照してください。したがって、それらは非常に短い距離でのトリップで、数秒しかかかりませんでしたが、データでは10分から1.4時間かかりました。

R
train %>%
  filter(total_travel_time < 10 & trip_duration < 10000) %>%
  ggplot(aes(trip_duration)) +
  geom_density(fill = "red") +
  geom_vline(xintercept = c(600,5000), linetype = 2) +
  scale_x_log10()

unnamed-chunk-74-1.png

R
foo <- train %>%
  filter(total_travel_time < 10) %>%
  mutate(delay = trip_duration > 600)

p1 <- foo %>%
  ggplot(aes(vendor_id, group = delay)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") + 
  scale_y_continuous(labels=scales::percent) +
  ylab("relative frequencies") +
  facet_grid(~delay) +
  theme(legend.position = "none")

p2 <- foo %>%
  ggplot(aes(wday, group = delay)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") + 
  scale_y_continuous(labels=scales::percent) +
  ylab("relative frequencies") +
  facet_grid(~delay) +
  theme(legend.position = "none")

p3 <- foo %>%
  ggplot(aes(passenger_count, group = delay)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") + 
  scale_y_continuous(labels=scales::percent) +
  ylab("relative frequencies") +
  facet_grid(~delay) +
  theme(legend.position = "none")

p4 <- foo %>%
  ggplot(aes(total_distance, fill = delay)) +
  geom_density(alpha = 0.5) +
  coord_cartesian(xlim = c(0,100))

layout <- matrix(c(1,1,2,2,2,2,3,3,3,4,4,4),2,6,byrow=TRUE)
multiplot(p1, p2, p3, p4, layout=layout)

unnamed-chunk-75-1.png

わかったことは、

  • ベンダー2は、ベンダー1よりも長いtrip_durationsを持っています。
  • 第2の遅延は、全体的な分布と比較して、週末にもっと頻繁あるように見える。
  • 遅延のために乗客の数も増えていますので、タクシーに出入りする時間(そして荷物)と関係している?
  • 遅延の距離は長いが、100m未満の距離は10分以上かかることはありません。

要約すると、vendor_idで2つのグループに識別でき、空港の乗降車は、長時間にわたる異常な遅延に寄与する可能性が高いことがわかります。wdaypassenger_countのような他の要因は、追加のシグナルに寄与する可能性があります。

7 相関の概要

新しい特徴量を設計し、モデリングを開始する前に、相関行列を使用してパラメタ間の関係を視覚化します。そのためには、すべての特徴量を数値形式に変更する必要があります。視覚化は、corrplot関数を使用します。Corrplotは、プロットのスタイルを操作する上で大きな柔軟性を与えてくれます。
以下に示すのは、2つの特徴の各組み合わせに対する色分けされた相関係数です。簡単に言えば、これは2つの特徴量の関係なので、一方を変更すると予測可能な関係かどうかを示します。この係数がゼロに近づくほど、相関は弱くなります。1と-1の両方が完全相関であり、理想的なケースです(下のプロットでは濃い青と濃い赤)。
ここでは、特徴量がtrip_durationと相関するかどうか、そしてどの程度なのか、興味があります。しかし、潜在的な予測子が相互に関連付けられているかどうかを知りたいので、データセットの共線性を減らし、予測の汎化性を向上させることができます。

R
train %>%
  select(-id, -pickup_datetime, -dropoff_datetime, -jfk_dist_pick,
         -jfk_dist_drop, -lg_dist_pick, -lg_dist_drop, -date) %>%
  mutate(passenger_count = as.integer(passenger_count),
         vendor_id = as.integer(vendor_id),
         store_and_fwd_flag = as.integer(as.factor(store_and_fwd_flag)),
         jfk_trip = as.integer(jfk_trip),
         wday = as.integer(wday),
         month = as.integer(month),
         work = as.integer(work),
         lg_trip = as.integer(lg_trip),
         blizzard = as.integer(blizzard),
         has_snow = as.integer(has_snow),
         has_rain = as.integer(has_rain)) %>%
  select(trip_duration, speed, everything()) %>%
  cor(use="complete.obs", method = "spearman") %>%
  corrplot(type="lower", method="circle", diag=FALSE)

unnamed-chunk-76-1.png
直接関連していない特徴量を削除すると、(有効な)相関行列は次のようになります。 このプロットでは、薄い色は見る価値はあまりありません。

R
foo <- train %>%
  select(-id, -pickup_datetime, -dropoff_datetime, -jfk_dist_pick,
         -jfk_dist_drop, -lg_dist_pick, -lg_dist_drop, -date,
         -store_and_fwd_flag, -hour, -rain, -s_fall, -all_precip,
         -has_rain, -s_depth, -min_temp, -max_temp,
         -wday, -right_turns, -turns, -fast_speed_trip) %>%
  mutate(passenger_count = as.integer(passenger_count),
         vendor_id = as.integer(vendor_id),
         jfk_trip = as.integer(jfk_trip),
         lg_trip = as.integer(lg_trip),
         month = as.integer(month),
         work = as.integer(work),
         blizzard = as.integer(blizzard),
         has_snow = as.integer(has_snow)) %>%
  select(trip_duration, speed, everything())

foo %>
  cor(use="complete.obs", method = "spearman") %>%
  corrplot(type="lower", method="circle", diag=FALSE)

unnamed-chunk-77-1.png
わかったことは、

  • trip_durationとの最も強い相関は、OSRMデータから導出されたtotal_distanceおよびtotal_travel_timeです。すでに上で見たように、distはお互いに、total_travel_timeと相関しています。おそらくnumber_of_stepsを通るターン数もtrip_durationに影響を与えています。
  • trip_durationへの別の影響は、jfk_triplg_tripで見られます。いずれかの空港への移動を示します。同様のことは、空港の移動と平均速度にも当てはまります。
  • ピックアップとドロップオフ座標は相関がありますが、これはやや困惑します。南西から北東に伸びるマンハッタンの形で部分的に説明されるかもしれません。もう1つは、マンハッタンの下部または上部の中での短い移動かもしれません。
  • vendor_idは、ベンダー2が6人以上の乗客を伴う(5つの)すべての移動を有するため、passenger_countと相関。
  • 気象デー中で、最低気温と最高気温はお互いに強い相関関係を示しており、年の月は直感的に理解できます。また、降水量を測定するさまざまな方法は、お互いに自然に相関しています。雪の場合は、月の値があります。 私たちのモデルでは、これらの関係を見ることができるはずです。

8 分類までの道

課題の本質は回帰問題です。一連の特徴量からtrip_durationという連続変数を予測する必要があります。このカーネルでは、引き続き回帰問題として扱いますが、このセクションでは分類問題に変えられる可能性について簡単に説明します。
分類では、目標変数は離散値の(少数の)数しか取ることができません。しばしば、私たちはバイナリの結果に直面します。例えば、人気のあるタイタニックのコンペでは「生存」か「死亡」とかです。 ここでは、どの要因がtrip_durationsの長短に寄与しているかを調べる必要があります。
このノートブックの最初のtrip_durationsの分布を思い出してください。今度は、怪しい偽の移動が取り除かれて、きれいになっています。 ここでは、短、中、長の3つの部分に分けて説明します。

R
train_group <- train %>%
  mutate(tgroup = case_when(trip_duration < 3e2 ~ "fast",
                            trip_duration >= 3e2 & trip_duration <= 1.6e3 ~ "mid",
                            trip_duration > 1.6e3 ~ "slow"))

train_group %>%
  ggplot(aes(trip_duration, fill = tgroup)) +
  geom_histogram(bins = 300) +
  scale_x_log10() +
  scale_y_sqrt()

unnamed-chunk-78-1.png
ここでは、中間に興味がありません。目的は、速いか遅い旅行の特性を比較することです。

R
train_group <- train_group %>%
  filter(tgroup != "mid")

p1 <- train_group %>%
  ggplot(aes(wday, fill = tgroup)) +
  geom_bar(position = "fill") +
  theme(legend.position = "none")

p2 <- train_group %>%
  ggplot(aes(month, fill = tgroup)) +
  geom_bar(position = "fill") +
  theme(legend.position = "none")

p3 <- train_group %>%
  ggplot(aes(hour, fill = tgroup)) +
  geom_bar(position = "fill")

p4 <- train_group %>%
  ggplot(aes(jfk_trip, fill = tgroup)) +
  geom_bar(position = "fill") +
  theme(legend.position = "none")

p5 <- train_group %>%
  ggplot(aes(lg_trip, fill = tgroup)) +
  geom_bar(position = "fill") +
  theme(legend.position = "none")

p6 <- train_group %>%
  ggplot(aes(blizzard, fill = tgroup)) +
  geom_bar(position = "fill") +
  theme(legend.position = "none"

p7 <- train_group %>%
  ggplot(aes(work, fill = tgroup)) +
  geom_bar(position = "fill") +
  theme(legend.position = "none")

layout <- matrix(c(1,1,2,2,3,3,3,3,4,5,6,7),3,4,byrow=TRUE)
multiplot(p1, p2, p3, p4, p5, p6, p7, layout=layout)

unnamed-chunk-79-1.png
わかったことは、
- 空港の特徴量jfk_triplg_tripは、期待どおりに動作しています。また、ブリザードの特徴量も、有望な増加を示しています。
- 週中、週末と平日の間には顕著な違いがあります。このトレンドに応じてwday特徴量をエンコードする価値があります。毎月の傾向は、徐々に「緩やかな」割合を示します。
- 日中は、散布図とヒートマップで前に見た夜間とビジネス時間の差が再びでています。これは特徴量エンジニアリングで部分的に解決されています。
ゆっくりとした移動と速い移動のピックアップ座標を視覚化するために、素晴らしいleafretパッケージを使用して別のマップを作成します。 これらのマップは外部API呼び出しを必要とせず、Kaggleカーネルに完全に含まれています。 私たちはヒートマップを使用して2つのグループの移動を視覚化します。ヒートマップは、リーフレットのエクストラパッケージに追加します。

R
pick_fast <- train_group %>%
  filter(tgroup == "fast") %>%
  select(long = pickup_longitude, lat = pickup_latitude)

pick_slow <- train_group %>%
  filter(tgroup == "slow") %>%
  select(long = pickup_longitude, lat = pickup_latitude)

leaflet(pick_fast) %>% 
  addProviderTiles(providers$CartoDB.DarkMatter, group = "Dark map") %>%
  addProviderTiles(providers$CartoDB.Positron, group = "Light map") %>%
  addScaleBar() %>%
  setView(-73.95, 40.74, zoom = 11) %>%
  addHeatmap(max = 140, gradient = "Reds", radius = 12,
             minOpacity = 0.15, blur = 15, group = "Fast trips") %>%
  addHeatmap(lng = pick_slow$long, lat = pick_slow$lat,
             max = 100, gradient = "Blues", radius = 12, 
             minOpacity = 0.25, blur = 15, group = "Slow trips") %>%
  addLayersControl(
    baseGroups = c("Dark Map", "Light map"),
    overlayGroups = c("Fast trips", "Slow trips"),
    options = layersControlOptions(collapsed = FALSE)
  )

スクリーンショット 0030-09-08 10.47.56.png
わかったこと

  • 予想どおり、空港は低速の移動地として際立っています。興味深いことに、空港周辺の同様な状態です。
  • ハドソン川とジャージーのニューアーク空港のすぐ近くには、数多くの高速の移動のピックアップがあります。
  • 興味深いことに、空港の場所を除いて、高速の移動のピックアップは、遅いピックアップよりも空間的に拡張されているように見えます。 分類問題のために、別の種類の視覚化を説明したいと思います。これらのプロットは、alluvialパッケージを介して利用できるようになっており、フローチャートと棒グラフのミックスの一種であり、ターゲットカテゴリがさまざまな個別の特徴量にどのように関連しているかを示します。 ```R:R allu_train <- train_group %>% group_by(tgroup, work, wday, jfk_trip) %>% count() %>% ungroup

alluvial(allu_train %>% select(-n),
freq=allu_train$n, border=NA,
col=ifelse(allu_train$tgroup == "fast", "red", "blue"),
cex=0.75,
hide = allu_train$n < 150,
ordering = list(
order(allu_train$tgroup=="fast"),
NULL,
NULL,
NULL))
```
unnamed-chunk-81-1.png
このプロットでは、ブロックの垂直サイズとストライプの幅(「alluvia」と呼ばれます)は頻度数に比例します。 同じ名前のパラメータを使って最も低い頻度数のalluviaを隠すことにしました。 カテゴリブロック内では、データセットの特定の側面を強調するために、alluviaの垂直レイヤを並べ替えることができます。 こでは、速いグループと遅いグループがどのようにして異なるカテゴリにどのように関係しているかを見れます。
わかったことは、
- 仕事時間外(仕事==FALSE)の移動の大半は速く、赤いalluviaのより大きな割合から見ることができます。
- 日曜、土曜、日曜日は、平日よりも早い移動の割合が高い。さらに、SatとSunは、対応するボックスにalluviaが接続されていないため、もちろんTRUEには作用しません。これは有用な一貫性チェックです。
- 空港の移動はほとんど遅いです。最初の場所では見えにくいため、(hideパラメータを使用して)いくつかの薄いalluviasを削除しました。
要約すると、この問題は、速い移動と遅い移動の分類の観点から見ると、特定の効果がより明確になり、たとえば空間分布や勤務時間の影響に関する有用な詳細がわかります。さらに、wday(曜日)や月(移動時間に対して無作為に)などの無作為に順序付けされたカテゴリ変数をエンコードする方法に関する有用なヒントを導き出します。これで分類については終わります。

9 シンプルなモデルと予測

最後に最終ステップとして、テストデータのターゲットであるtrip_durationを予測するために、探索的および技術的な洞察いかしてモデリングします。しかし、これは主にEDAカーネルであり、モデリングがメインではありません。このセクションでは、最適化したい人のための改善の余地がたくさんあります。XGBoostを使用しますが、他のアルゴリズムを試すこともできます。
私は、このコンペで時間が進むにつれて、モデルの選択とチューニングの複雑さに焦点が移ることを期待しています。スキルを実証したいマスターとグランドマスターから学ぶことができます。なので、このカーネルはモデルスタッキングには全く触れません。

9.1 準備

9.1.1 準備

テストデータに関連する特徴量を実際に学習していることを確認するために、訓練データとテストデータの時間的および空間的特性を簡単に比較します。 これは別の一貫性チェックです。 探索の前にこれを行うことができましたが、個人的な好みは、分析ができるだけ公平でないように、試験データセットを見る前にトレーニングデータを最初に調べることです。 関連する2つの比較プロットは次のとおりです。

R
foo <- combine %>%
  mutate(date = date(ymd_hms(pickup_datetime))) %>%
  group_by(date, dset) %>%
  count() %>%
  ungroup()

foo %>%
  ggplot(aes(date, n/1e3, color = dset)) +
  geom_line(size = 1.5) +
  labs(x = "", y = "Kilo trips per day")

unnamed-chunk-82-1.png

R
pick_good <- combine %>%
  filter(pickup_longitude > -75 & pickup_longitude < -73) %>%
  filter(pickup_latitude > 40 & pickup_latitude < 42)

pick_good <- sample_n(pick_good, 5e3)

pick_good %>%
  ggplot(aes(pickup_longitude, pickup_latitude, color = dset)) +
  geom_point(size=0.1, alpha = 0.5) +
  coord_cartesian(xlim = c(-74.02,-73.77), ylim = c(40.63,40.84)) +
  facet_wrap(~ dset) +
  #guides(color = guide_legend(override.aes = list(alpha = 1, size = 4))) +
  theme(legend.position = "none")

unnamed-chunk-83-1.png
訓練データとテストデータは、実際には同じ時間と地域をカバーしていることがわかります。

9.1.2 データのコード化

多くの分類器はカテゴリ値を扱うことができないので、ここでは選択された特徴を整数列に変換するようにフォーマットします。エンコーディングでは、ここまで得られた洞察を含む探索的な知識を利用します。これは、大部分の分類器が整数特徴(すなわち、ターゲットへの影響に関して1 <2 <3)に対して自然順序付けを想定することです。別の方法として、ワンホットエンコーディングを使用する方法があります。
訓練データとテストデータとの間で符号化を均一にするために、いくつかの要因レベルに違いがある場合には、組み合わせたデータに対して符号化を実行することが有用な方法です。たとえば、passenger_countが本当であった場合、ここでは訓練とテストデータにはさまざまなレベルが含まれます(テストには7,8回の旅程がないため)。
この再コーディングの段階は、フィーチャエンジニアリングの冒頭で行ったこととほぼ同じです。原理的には、結合データセットですでに作業していた可能性がありましたが、単純なままにしておくために、元の訓練サンプルのみを使用しました。このカーネルを使い、効率を上げるためにカーネルを改善することを大変歓迎します。エンジニアリング/フォーマットコードブロックに魔法のようなものはありません。

R
jfk_coord <- tibble(lon = -73.778889, lat = 40.639722)
la_guardia_coord <- tibble(lon = -73.872611, lat = 40.77725)

# derive distances
pick_coord <- combine %>%
  select(pickup_longitude, pickup_latitude)
drop_coord <- combine %>%
  select(dropoff_longitude, dropoff_latitude)

combine$dist <- distCosine(pick_coord, drop_coord)
combine$bearing = bearing(pick_coord, drop_coord)
combine$jfk_dist_pick <- distCosine(pick_coord, jfk_coord)
combine$jfk_dist_drop <- distCosine(drop_coord, jfk_coord)
combine$lg_dist_pick <- distCosine(pick_coord, la_guardia_coord)
combine$lg_dist_drop <- distCosine(drop_coord, la_guardia_coord)

# add dates
combine <- combine %>%
  mutate(pickup_datetime = ymd_hms(pickup_datetime),
         dropoff_datetime = ymd_hms(dropoff_datetime),
         date = date(pickup_datetime)
  )

# add weather info
foo <- weather %>%
  select(date, rain, s_fall, all_precip, has_snow, has_rain, s_depth, max_temp, min_temp)
combine <- left_join(combine, foo, by = "date")

# add fast routes
foo <- fastest_route %>%
  select(id, total_distance, total_travel_time, number_of_steps,
         step_direction, step_maneuvers) %>%
  mutate(fastest_speed = total_distance/total_travel_time*3.6,
         left_turns = str_count(step_direction, "left"),
         right_turns = str_count(step_direction, "right"),
         turns = str_count(step_maneuvers, "turn")
         ) %>%
  select(-step_direction, -step_maneuvers)

combine <- left_join(combine, foo, by = "id")

# reformat to numerical and recode levels
combine <- combine %>%
  mutate(store_and_fwd_flag = as.integer(factor(store_and_fwd_flag)),
         vendor_id = as.integer(vendor_id),
         month = as.integer(month(pickup_datetime)),
         wday = wday(pickup_datetime, label = TRUE),
         wday = as.integer(fct_relevel(wday, c("Sun", "Sat", "Mon", "Tues", "Wed", "Thurs", "Fri"))),
         hour = hour(pickup_datetime),
         work = as.integer( (hour %in% seq(8,18)) & (wday %in% c("Mon","Tues","Fri","Wed","Thurs")) ),
         jfk_trip = as.integer( (jfk_dist_pick < 2e3) | (jfk_dist_drop < 2e3) ),
         lg_trip = as.integer( (lg_dist_pick < 2e3) | (lg_dist_drop < 2e3) ),
         has_rain = as.integer(has_rain),
         has_snow = as.integer(has_snow),
         blizzard = as.integer( !( (date < ymd("2016-01-22") | (date > ymd("2016-01-29")))) )
         )

分類の結果に従って平日(wday)をコード化することに気付いたかもしれません。都合の良いことに、この月は既に低速移動の割合が増えています。一貫性を確認しておきます。

R
glimpse(combine)
## Observations: 2,083,778
## Variables: 41
## $ id                 <chr> "id2875421", "id2377394", "id3858529", "id3...
## $ vendor_id          <int> 2, 1, 2, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2, 2...
## $ pickup_datetime    <dttm> 2016-03-14 17:24:55, 2016-06-12 00:43:35, ...
## $ dropoff_datetime   <dttm> 2016-03-14 17:32:30, 2016-06-12 00:54:38, ...
## $ passenger_count    <int> 1, 1, 1, 1, 1, 6, 4, 1, 1, 1, 1, 4, 2, 1, 1...
## $ pickup_longitude   <dbl> -73.98215, -73.98042, -73.97903, -74.01004,...
## $ pickup_latitude    <dbl> 40.76794, 40.73856, 40.76394, 40.71997, 40....
## $ dropoff_longitude  <dbl> -73.96463, -73.99948, -74.00533, -74.01227,...
## $ dropoff_latitude   <dbl> 40.76560, 40.73115, 40.71009, 40.70672, 40....
## $ store_and_fwd_flag <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ trip_duration      <int> 455, 663, 2124, 429, 435, 443, 341, 1551, 2...
## $ dset               <fct> train, train, train, train, train, train, t...
## $ dist               <dbl> 1500.1995, 1807.5298, 6392.2513, 1487.1625,...
## $ bearing            <dbl> 99.932546, -117.063997, -159.608029, -172.7...
## $ jfk_dist_pick      <dbl> 22315.02, 20258.98, 21828.54, 21461.51, 236...
## $ jfk_dist_drop      <dbl> 21025.4008, 21220.9775, 20660.3899, 21068.1...
## $ lg_dist_pick       <dbl> 9292.897, 10058.778, 9092.997, 13228.107, 8...
## $ lg_dist_drop       <dbl> 7865.248, 11865.578, 13461.012, 14155.920, ...
## $ date               <date> 2016-03-14, 2016-06-12, 2016-01-19, 2016-0...
## $ rain               <dbl> 0.29, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0...
## $ s_fall             <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0...
## $ all_precip         <dbl> 0.29, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0...
## $ has_snow           <int> 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ has_rain           <int> 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0...
## $ s_depth            <dbl> 0.00, 0.00, 0.01, 0.00, 0.00, 6.00, 0.00, 0...
## $ max_temp           <int> 51, 83, 28, 48, 55, 39, 78, 66, 87, 79, 63,...
## $ min_temp           <int> 40, 62, 16, 30, 38, 28, 63, 54, 73, 63, 50,...
## $ total_distance     <dbl> 2009.1, 2513.2, 11060.8, 1779.4, 1614.9, 13...
## $ total_travel_time  <dbl> 164.9, 332.0, 767.6, 235.8, 140.1, 189.4, 1...
## $ number_of_steps    <int> 5, 6, 16, 4, 5, 5, 5, 17, 2, 13, 6, 9, 4, 4...
## $ fastest_speed      <dbl> 43.86149, 27.25157, 51.87452, 27.16641, 41....
## $ left_turns         <int> 1, 2, 5, 2, 2, 1, 1, 4, 0, 7, 2, 6, 1, 1, 4...
## $ right_turns        <int> 1, 2, 7, 1, 2, 3, 3, 9, 1, 5, 2, 2, 1, 1, 3...
## $ turns              <int> 1, 2, 9, 1, 3, 3, 2, 6, 0, 9, 3, 5, 2, 2, 5...
## $ month              <int> 3, 6, 1, 4, 3, 1, 6, 5, 5, 3, 5, 5, 2, 6, 5...
## $ wday               <int> 3, 1, 6, 4, 2, 2, 5, 2, 5, 7, 6, 1, 5, 4, 5...
## $ hour               <int> 17, 0, 11, 19, 13, 22, 22, 7, 23, 21, 22, 1...
## $ work               <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ jfk_trip           <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ lg_trip            <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ blizzard           <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...

いいですね。 数値以外の特徴量は、idpickup_datetimedropoff_datetimeおよびdateです。これは、どんな場合でも、dsetと一緒に削除されます。これは、訓練データとテストデータを再び分離するために使用します。

9.1.3 特徴量選択、評価指標、バリデーション、交差検証

データのすべての特徴量が役立つわけではありません。ここでは意味のある変数だけを取り上げて、たとえばidや大部分の天気の特徴量を削除します。
原則として、すべての特徴量を含めることができ、有用な特徴量選択をモデリングアルゴリズムに任せることができます。しかし、私たちのケースでは、(a)既存のものからいくつかの特徴量を設計し、(b)強い相関を持つ追加の特徴量を持っています。訓練セット内で著しい共線性を引き起こす可能性があり、個々の特徴量の影響の観点から、モデルの結果を解釈するのがより困難になります。さらに、実用的な観点から見ると、強く相関している特徴はあまり新しい情報を追加せず、統計的理由から削除されるべきです。
特徴量選択は、通常、反復プロセスであり、少数または多数の特徴量のいずれかを使用して初期モデルを実行し、前回の実行の結果に基づいていくつかの特徴量を段階的に追加または削除します。ここではモデル適合ステップを1つだけ示します。あなた自身の分析の出発点としてそれを使用し、それを他のカーネルからの洞察と組み合わせてみてください。探索的分析の結果、最終的な相関行列を基に特徴量選択を行います。

R
# Specific definitions:
#---------------------------------
# predictor features
train_cols <- c("total_travel_time", "total_distance", "hour", "dist",
                "vendor_id", "jfk_trip", "lg_trip", "wday", "month",
                "pickup_longitude", "pickup_latitude", "bearing", "lg_dist_drop")
# target feature
y_col <- c("trip_duration")
# identification feature
id_col <- c("id") 
# auxilliary features
aux_cols <- c("dset")
# cleaning features
clean_cols <- c("jfk_dist_drop", "jfk_dist_pick")
#---------------------------------
# General extraction
#---------------------------------
# extract test id column
test_id <- combine %>%
  filter(dset == "test") %>%
  select_(.dots = id_col)

# all relevant columns for train/test
cols <- c(train_cols, y_col, aux_cols, clean_cols)
combine <- combine %>%
  select_(.dots = cols)

# split train/test
train <- combine %>%
  filter(dset == "train") %>%
  select_(.dots = str_c("-",c(aux_cols)))

test <- combine %>%
  filter(dset == "test") %>%
  select_(.dots = str_c("-",c(aux_cols, clean_cols, y_col)))
#---------------------------------

このタクシーチャレンジでは、評価指標はRMSLE(Root Mean Squared Logarithmic Error)です。 本質的に、これは、ログ空間における予測対データ偏差を最適化することを意味します。 これには、大きな個別の異常値が線形メトリックの場合と同じくらい大きな重みを持たないという利点があります。
モデルフィッティングの評価メトリックを簡単にシミュレートするために、trip_durationを対数で置き換えます。 (未定義のログ(0)を避けるために+1が追加されており、予測ファイルに対してこの1秒を削除することを忘れないようにする必要がありますが、そうしないからとはいえ大きな違いはありません)。

R
train <- train %>%
  mutate(trip_duration = log(trip_duration + 1))

私たちのモデルがどれほどよく汎化しているかを評価するために、クロスバリデーションを実行し、トレーニングデータを訓練データとテストデータに分割します。これにより、アルゴリズムが見ていないサンプルでモデル性能を評価することができます。caretパッケージのツールを使用して、データを80/20に分割しました。私たちはこのような簡単な作業のためにこのパッケージを実際には必要としませんが、私はそれに言及することを含めています。caretは価値のある多目的MLパッケージです。

R
set.seed(4321)

trainIndex <- createDataPartition(train$trip_duration, p = 0.8, list = FALSE, times = 1)

train <- train[trainIndex,]
valid <- train[-trainIndex,]

最後に、訓練データでのみ慎重なクリーニングを実行します。観測データをトレーニングデータから削除することにより、モデルを現実の「理想化された」バージョンにオーバーフィットさせるリスクを常に負います。特に、現在のデータセットでは、多くの値が明白ではないため、この理想化されたモデルは、見た目に似ていないデータには一般化できない可能性があります。このため、バリデーションサンプルは変更していないため、モデルスコアの間に顕著な偏差があれば、あまりフィットしないことがわかります。
なぜきれいにする必要があるのか?私の意見では、モデルを堅牢にするために、外れ値を削除する必要があります。これらの異常値がテストデータにも存在する可能性がありますが、予測する方法があれば、結局は異常値ではない可能性があります。とにかくそれを予測できなければ、とにかくそれを考慮する方法はありません。最後に、解答は、与えられたテストセットを可能な限りうまくフィットさせることを常に目指すとは限りませんが、極端な値を無視できるより実用的なターゲットを持つ場合もある分析(および評価メトリック)の目標にもよります。
ここでは、24時間以上のtrip_durationと、NYCから遠く離れた数のデータポイントだけを削除します。

R
valid <- valid %>%
  select_(.dots = str_c("-",c(clean_cols)))

train <- train %>%
  filter(trip_duration < 24*3600,
         jfk_dist_pick < 3e5 & jfk_dist_drop < 3e5
         ) %>%
  select_(.dots = str_c("-",c(clean_cols)))

9.2 XGBoostのパラメタとフィッティング

モデルフィッティングでは、みんなが好きなXGBoost - eXtreme Gradient Boostingを使用します。これはKaggleの人気度が高いdecision-tree gradient boostingアルゴリズムです。このカーネルでは、かなり単純なフィッティング戦略をとります。パラメタを最適化し、モデルフィッティングにもっと焦点を当てた他のカーネルからインスピレーションを得てください。また、Meganのベースライン「ランダムフォレスト」などの他の分類器を試して、特徴量をどのように扱うかを確認してください。
Gradient Boostingは、弱学習器を前の学習ステップの結果に逐次適用することによって段階的な改善していくものです。Gradient Boostingは、逐次的にアルゴリズムをトレーニングすることにより、ロス関数(評価メトリックによる)を最小化することに焦点を当てています。

R
#convert to XGB matrix
foo <- train %>% select(-trip_duration)
bar <- valid %>% select(-trip_duration)

dtrain <- xgb.DMatrix(as.matrix(foo),label = train$trip_duration)
dvalid <- xgb.DMatrix(as.matrix(bar),label = valid$trip_duration)
dtest <- xgb.DMatrix(as.matrix(test))

ここでは、XGBoostを制御するハイパパラメタを定義します。詳細はこちらを参照してください。 ウォッチリストパラメータは、トレーニングとテストのサンプルメトリクスの両方に注目するようにアルゴリズムに指示します。これらのパラメタはパフォーマンスに最適化されていないため、自由にえらんでください。

R
xgb_params <- list(colsample_bytree = 0.7, #variables per tree 
                   subsample = 0.7, #data subset per tree 
                   booster = "gbtree",
                   max_depth = 5, #tree levels
                   eta = 0.3, #shrinkage
                   eval_metric = "rmse", 
                   objective = "reg:linear",
                   seed = 4321
                   )

watchlist <- list(train=dtrain, valid=dvalid)

ここで分類器をトレーニングします。つまり、トレーニングデータに適合させます。再現性を保証するために、ここにシードを設定します。このモデルをカーネル環境で実行するには、広範なEDA実行時間を考慮して、ラーニングを60回のサンプルラウンドに制限します。

R
set.seed(4321)

gb_dt <- xgb.train(params = xgb_params,
                   data = dtrain,
                   print_every_n = 5,
                   watchlist = watchlist,
                   nrounds = 60)

## [1]  train-rmse:4.227815 valid-rmse:4.226777 
## [6]  train-rmse:0.836443 valid-rmse:0.835214 
## [11] train-rmse:0.445533 valid-rmse:0.443858 
## [16] train-rmse:0.420997 valid-rmse:0.419356 
## [21] train-rmse:0.413746 valid-rmse:0.412055 
## [26] train-rmse:0.410272 valid-rmse:0.408660 
## [31] train-rmse:0.407118 valid-rmse:0.405404 
## [36] train-rmse:0.405107 valid-rmse:0.403487 
## [41] train-rmse:0.403133 valid-rmse:0.401559 
## [46] train-rmse:0.401274 valid-rmse:0.399769 
## [51] train-rmse:0.400081 valid-rmse:0.398557 
## [56] train-rmse:0.398853 valid-rmse:0.397263 
## [60] train-rmse:0.398211 valid-rmse:0.396623

フィッティング後、モデルの性能を評価するために5foldCVを実行しています。また、このステージはより多くのラウンドでKaggleランタイム制限を超えますので、ここでは原則を示すために15回のサンプルラウンドに限定しています。XGBoostのパラメタに応じて、分析には少なくとも100を使用する必要があります。早期停止パラメタは、追加ステップを経てモデルを改善できない場合に停止させるものです。

R
xgb_cv <- xgb.cv(xgb_params,dtrain,early_stopping_rounds = 10, nfold = 5, nrounds=15)

## [1]  train-rmse:4.227837+0.000287    test-rmse:4.227776+0.000829 
## Multiple eval metrics are present. Will use test_rmse for early stopping.
## Will train until test_rmse hasn't improved in 10 rounds.
## 
## [2]  train-rmse:2.978744+0.000522    test-rmse:2.978735+0.001204 
## [3]  train-rmse:2.110670+0.000484    test-rmse:2.110653+0.001094 
## [4]  train-rmse:1.514938+0.003098    test-rmse:1.515065+0.003327 
## [5]  train-rmse:1.107978+0.002013    test-rmse:1.108181+0.002499 
## [6]  train-rmse:0.839585+0.001026    test-rmse:0.839878+0.002332 
## [7]  train-rmse:0.666475+0.001609    test-rmse:0.666924+0.002050 
## [8]  train-rmse:0.560617+0.000942    test-rmse:0.561197+0.002019 
## [9]  train-rmse:0.499846+0.001079    test-rmse:0.500549+0.001802 
## [10] train-rmse:0.466047+0.001115    test-rmse:0.466952+0.001702 
## [11] train-rmse:0.446989+0.000477    test-rmse:0.447998+0.002235 
## [12] train-rmse:0.436274+0.000600    test-rmse:0.437394+0.002416 
## [13] train-rmse:0.429310+0.000337    test-rmse:0.430505+0.002351 
## [14] train-rmse:0.425255+0.000576    test-rmse:0.426569+0.002368 
## [15] train-rmse:0.422719+0.000467    test-rmse:0.424148+0.002391

9.3 変数重要度

トレーニングの後、モデルにとって最も重要な特徴量を確認します。これは、モデルの重要な特徴量を段階的に特定する反復プロセスの出発点となります。ここでは、これらの特徴量を単に視覚化します。

R
imp_matrix <- as.tibble(xgb.importance(feature_names = colnames(train %>% select(-trip_duration)), model = gb_dt))

imp_matrix %>%
  ggplot(aes(reorder(Feature, Gain, FUN = max), Gain, fill = Feature)) +
  geom_col() +
  coord_flip() +
  theme(legend.position = "none") +
  labs(x = "Features", y = "Importance")

unnamed-chunk-94-1.png
意外にも、total_distancetotal_travel_timeのようなOSRMの特徴量は、trip_durationの最も重要な予測変数のようです。それ以外にも、wdaybearingのような特徴量エンジニアリングで設計した特徴量も悪くない。

9.4 予測と提出ファイル

CVスコアに満足したら、モデルを使用してテストデータの予測を行い、提出ファイルを作詞します。そして、Kaggleリーダーボードに提出します。これにより、独立したデータセットに基づくパフォーマンススコアが得られます。

R
test_preds <- predict(gb_dt,dtest)

pred <- test_id %>%
  mutate(trip_duration = exp(test_preds) - 1)

pred %>% write_csv('submit.csv')

この提出ファイルを再度確認して、不正な形式ではないか確認します。その形状と内容を最初に読み込んだサンプル提出ファイルと比較しており、予測値の分布を近似比較のトレーニング値と比較してプロットしています。

R
identical(dim(sample_submit),dim(pred))
## [1] TRUE

glimpse(sample_submit)
## Observations: 625,134
## Variables: 2
## $ id            <chr> "id3004672", "id3505355", "id1217141", "id215012...
## $ trip_duration <int> 959, 959, 959, 959, 959, 959, 959, 959, 959, 959...
glimpse(pred)
## Observations: 625,134
## Variables: 2
## $ id            <chr> "id3004672", "id3505355", "id1217141", "id215012...
## $ trip_duration <dbl> 826.6543, 419.2442, 371.0309, 992.0790, 291.4171...
bind_cols(sample_submit, pred) %>% head(5)
## # A tibble: 5 x 4
##   id        trip_duration id1       trip_duration1
##   <chr>             <int> <chr>              <dbl>
## 1 id3004672           959 id3004672            827
## 2 id3505355           959 id3505355            419
## 3 id1217141           959 id1217141            371
## 4 id2150126           959 id2150126            992
## 5 id1598245           959 id1598245            291

foo <- train %>%
  select(trip_duration) %>%
  mutate(dset = "train",
         trip_duration = exp(trip_duration) - 1)

bar <- pred %>%
  mutate(dset = "predict")

bind_rows(foo, bar) %>%
  ggplot(aes(trip_duration, fill = dset)) +
  geom_density(alpha = 0.5) +
  scale_x_log10()

unnamed-chunk-96-1.png
予測ファイルとtrip_durationの範囲が正常であることがわかりました。現在の状態では、このモデルは、あなたが投稿時に投稿のトップ40%くらいです。明らかに、これはモデリングプロセスを探索し改善するための出発点にすぎません ここにいくつかの提案があります。

  • より多くのラウンドでトレーニングとクロスバリデーションを実行します。 これだけであなたのスコアがすぐに向上します。
  • モデルにさらなる特徴量を追加し、多かれ少なかれクリーニングの影響を探る。
  • 別のXGBoostパラメタを試して、CVスコアにどのように影響するかを確認する。

後はあなたしだいです。広範なEDAが、あなた自身の分析を開始するためのアイデアであり、データの有用な洞察を与えてくれることを願っています。最高の成功をお祈りします。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account