概要
前回の記事に続きRの基本的な利用方法を解説します。今回は、効率的なデータ前処理を実現するパッケージ群 tidyverse の使い方を具体的なデータセットを用いて解説します。
必要なパッケージの準備
分析を始める前に必要なパッケージをインストールして読み込みます。
# パッケージが未インストールの場合のみ実行
install.packages("tidyverse") # データ操作や可視化のためのパッケージ群
install.packages("nycflights13") # 分析に使用するデータセット
install.packages("MASS") # 重回帰分析で変数選択を行うために使用
# ライブラリの読み込み
library(tidyverse)
library(nycflights13)
library(MASS)
nycflights13データセット
今回はnycflights13というデータセットを使用します。これは2013年にニューヨーク市の3つの主要空港（JFK, LGA, EWR）から出発した、全フライトのデータが含まれています。データハンドリングや可視化の練習に最適なデータセットです。
このパッケージには、主要なflightsテーブルのほか、関連する複数のテーブルが含まれています。
- flights: 全てのフライト情報（出発・到着時刻、遅延、航空会社など）
- weather: 1時間ごとの気象データ
- planes: 航空機の機体情報
- airports: 空港の位置情報など
- airlines: 航空会社名とコードの対応表
テーブル間の関係
flightsテーブルのcarrier列をキーにしてairlinesテーブルを結合すれば、航空会社コードを正式名称に変換できます。
このデータセットを使えば、「どの航空会社が最も遅延が多いか？」「気象条件はフライト遅延にどう影響するか？」といった分析ができます。
# flightsデータを変数に格納
flights <- nycflights13::flights
flightsテーブルについて
| 列名 | 内容 |
|---|---|
| year, month, day | 出発年、月、日 |
| dep_time, arr_time | 実際の出発時間、到着時間(HHMM or HMM)。ローカルのタイムゾーン |
| sched_dep_time,sched_arr_time | 予定されていた出発時間、到着時間(HHMM or HMM)。ローカルのタイムゾーン |
| dep_delay,arr_delay | 出発、到着の誤差(分単位)。マイナスは予定よりも早めの出発・到着。 |
| carrier | 2文字の航空会社コード。別のテーブルに航空会社マスターあり。 |
| flight | フライト番号 |
| tailnum | 航空機番号。別のテーブルに航空機マスターあり。 |
| origin, dest | 出発地および到着地 |
| air_time | フライト時間(分単位) |
| distance | 到着地までの距離(マイル単位)。 |
| hour, minute | 予定されていた出発時間の時と分。 |
| time_hour | 予定されていたフライト時間をPOSIXct dateで。 |
tidyverseによる基本操作
tidyverseの中核であるdplyrパッケージの関数を使い、データ操作の基本を見ていきます。
%>% (パイプ演算子) を使うことで、処理を左から右へ流れるように記述できます。
1. 特定の列を抽出 (select)
flightsデータからyear列だけを抽出します。
flights %>%
select(year)
## # A tibble: 336,776 x 1
## year
## <int>
## 1 2013
## 2 2013
## 3 2013
## 4 2013
## 5 2013
## 6 2013
## 7 2013
## 8 2013
## 9 2013
## 10 2013
## # … with 336,766 more rows
2. 特定の行を抽出 (filter)
出発地 (origin) が "JFK" のフライトのみを抽出します。
flights %>%
filter(origin == "JFK")
# A tibble: 111,279 × 19
year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time arr_delay carrier flight tailnum origin dest air_time distance hour minute time_hour
<int> <int> <int> <int> <int> <dbl> <int> <int> <dbl> <chr> <int> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dttm>
1 2013 1 1 542 540 2 923 850 33 AA 1141 N619AA JFK MIA 160 1089 5 40 2013-01-01 05:00:00
2 2013 1 1 544 545 -1 1004 1022 -18 B6 725 N804JB JFK BQN 183 1576 5 45 2013-01-01 05:00:00
3 2013 1 1 557 600 -3 838 846 -8 B6 79 N593JB JFK MCO 140 944 6 0 2013-01-01 06:00:00
4 2013 1 1 558 600 -2 849 851 -2 B6 49 N793JB JFK PBI 149 1028 6 0 2013-01-01 06:00:00
5 2013 1 1 558 600 -2 853 856 -3 B6 71 N657JB JFK TPA 158 1005 6 0 2013-01-01 06:00:00
6 2013 1 1 558 600 -2 924 917 7 UA 194 N29129 JFK LAX 345 2475 6 0 2013-01-01 06:00:00
7 2013 1 1 559 559 0 702 706 -4 B6 1806 N708JB JFK BOS 44 187 5 59 2013-01-01 05:00:00
8 2013 1 1 606 610 -4 837 845 -8 DL 1743 N3739P JFK ATL 128 760 6 10 2013-01-01 06:00:00
9 2013 1 1 611 600 11 945 931 14 UA 303 N532UA JFK SFO 366 2586 6 0 2013-01-01 06:00:00
10 2013 1 1 613 610 3 925 921 4 B6 135 N635JB JFK RSW 175 1074 6 10 2013-01-01 06:00:00
# ℹ 111,269 more rows
# ℹ Use `print(n = ...)` to see more rows
3. グループ化して要約 (group_by, summarise)
空港ごと (origin) のフライト便数を集計します。
flights %>%
group_by(origin) %>%
summarise(count=n())
## # A tibble: 3 x 2
## origin count
## <chr> <int>
## 1 EWR 120835
## 2 JFK 111279
## 3 LGA 104662
4. 新しい列を追加 (mutate)
出発地と到着地を結合したorigin_destという新しい列を追加します。
flights %>%
mutate(origin_dest = paste(origin, dest, sep = " -> ")) %>%
select(origin, dest, origin_dest) # 確認のために一部の列だけ表示
## # A tibble: 336,776 x 3
## origin dest origin_dest
## <chr> <chr> <chr>
## 1 EWR IAH EWR -> IAH
## 2 LGA IAH LGA -> IAH
## 3 JFK MIA JFK -> MIA
## # ... with 336,763 more rows
5. データの結合 (left_join)
flightsデータに、carrier列をキーとしてairlinesテーブルを結合し、航空会社名を追加します。
flights %>%
select(carrier, flight) %>% # 分かりやすくするため列を絞る
left_join(airlines, by = "carrier")
## # A tibble: 336,776 x 5
## year origin dest carrier name
## <int> <chr> <chr> <chr> <chr>
## 1 2013 EWR IAH UA United Air Lines Inc.
## 2 2013 LGA IAH UA United Air Lines Inc.
## 3 2013 JFK MIA AA American Airlines Inc.
## 4 2013 JFK BQN B6 JetBlue Airways
## 5 2013 LGA ATL DL Delta Air Lines Inc.
## 6 2013 EWR ORD UA United Air Lines Inc.
## 7 2013 EWR FLL B6 JetBlue Airways
## 8 2013 LGA IAD EV ExpressJet Airlines Inc.
## 9 2013 JFK MCO B6 JetBlue Airways
## 10 2013 LGA ORD AA American Airlines Inc.
## # … with 336,766 more rows
tidyverseとネイティブなRの比較
tidyverseのコードがいかに直感的か、Rの基本的な書き方（ネイティブR）と比較してみます。
列の抽出
ネイティブR
ブラケット [] を使い、[行, 列] の形式で指定します。列名を文字列として引用符で囲む必要があります。
flights[, "year", drop = FALSE]
Tidyverse
select関数を使い、列名を直接指定します。「flightsからyearを選ぶ」と、英文のように読めます。
flights %>%
select(year)
行の抽出
ネイティブR
ブラケット内で条件式を記述します。処理が複雑になると、可読性が低下することがあります
flights[flights$origin == "JFK", ]
Tidyverse
パイプ演算子 %>% を使い、「フィルターをかけて、列を選ぶ」という思考の順序で記述できます。
flights %>%
filter(origin == "JFK")
グループ集計
ネイティブR
aggregate関数を使いますが、引数の指定が少し複雑です。
aggregate(x = list(count = flights$origin),
by = list(origin = flights$origin),
FUN = length)
Tidyverse
「グループ化して (group_by)、要約する (summarise)」という直感的な流れで書けます。
flights %>%
group_by(origin) %>%
summarise(count = n())
列の追加
ネイティブR
$を使って列を追加するか、data.frame関数で新しいデータフレームを作成します。
# 元のデータフレームに列を追加する場合
flights_2 <- flights
flights_2$origin_dest <- paste(flights$origin, "_", flights$dest)
Tidyverse
mutate関数を使い、処理を連結する中で自然に列を追加できます。
flights_2 <- flights %>%
mutate(origin_dest = paste(origin, "_", dest))
tidyverseは、特に複数の処理を連続して行う場合に可読性が高く、コードが簡潔になるというメリットがあります。
応用：線形回帰で遅延要因を分析
flightsの出発遅延 (dep_delay) が、気象 (weather) データとどのような関係にあるかを線形回帰モデルで分析してみます。ここではstepAIC関数を使い、モデルにとって重要な変数を自動的に選択させます。
まず最初に、MASSとdplyrの selectが衝突するのでdplyrを優先させるためのコードを追加します。
# 1. 'conflicted' パッケージをインストールする（初回のみ実行が必要）
# install.packages("conflicted")
# 2. 'conflicted' パッケージを現在のRセッションに読み込む
library(conflicted)
# 3. 関数名の競合が発生した場合の優先順位を設定する
# ここでは、'select' という名前の関数は、'dplyr' パッケージのものを優先して使用するよう指定
conflict_prefer("select", "dplyr")
stepAICで線形モデルを作成し、モデルを評価します。
# flightsとweatherを結合し、モデルに使用するデータを準備
delay_weather <- flights %>%
select(dep_delay, time_hour, origin) %>%
inner_join(weather, by = c("time_hour", "origin")) %>%
select(-time_hour, -origin, -year, -month, -day, -hour) %>% # 不要な列や重複する列を削除
na.omit() # 欠損値を含む行を削除
# 全ての変数を投入したモデルを作成
lm_full <- lm(dep_delay ~ ., data = delay_weather)
# stepAICで変数選択を実行
lm_aic <- stepAIC(lm_full, trace = FALSE) # trace = FALSEで途中経過を非表示
# 結果の要約を表示
summary(lm_aic)
##
## Call:
## lm(formula = dep_delay ~ day + temp + dewp + humid + wind_dir +
## wind_speed + wind_gust + pressure + visib, data = delay_weather)
##
## Residuals:
## Min 1Q Median 3Q Max
## -59.15 -17.25 -11.73 -0.04 935.74
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 215.043479 14.201280 15.143 < 2e-16 ***
## day 0.140517 0.009725 14.449 < 2e-16 ***
## temp 0.435548 0.049726 8.759 < 2e-16 ***
## dewp -0.362058 0.054558 -6.636 3.22e-11 ***
## humid 0.418587 0.030646 13.659 < 2e-16 ***
## wind_dir -0.009423 0.001167 -8.077 6.67e-16 ***
## wind_speed -0.114636 0.034788 -3.295 0.000984 ***
## wind_gust 0.339962 0.030100 11.294 < 2e-16 ***
## pressure -0.223961 0.013137 -17.048 < 2e-16 ***
## visib -1.274536 0.095875 -13.294 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38.57 on 219416 degrees of freedom
## (787561 observations deleted due to missingness)
## Multiple R-squared: 0.02563, Adjusted R-squared: 0.02559
## F-statistic: 641.4 on 9 and 219416 DF, p-value: < 2.2e-16
結果の解釈
各変数の影響: 選択された変数はすべて統計的に有意（Pr(>|t|)の値が非常に小さい）でした。
視程 (visib): 係数が -1.275 と負の値です。これは 視界が良いほど、出発遅延が減少する（定刻通りか、むしろ早まる） 傾向があることを示します。
突風 (wind_gust) や 降水量 (precip): 係数が正の値です。これは突風が強い、または雨が降っているほど、出発が遅延する傾向があることを示します。
モデルの精度: 決定係数 (Multiple R-squared) が 0.026 (約2.6%) と非常に低い値です。これは、 「気象条件は遅延に影響するものの、出発遅延全体のばらつきのわずか2.6%しか説明できていない」 ことを意味します。遅延には、機材繰りや航空管制など、天候以外の要因がはるかに大きく影響していると推測できます。
おわりに
今回は、tidyverseを用いた基本的なデータ操作から、簡単な回帰分析までを紹介しました。ネイティブなRの書き方も重要ですが、tidyverseを使いこなすことで、複雑なデータ処理をより効率的かつ直感的に行うことができます。