R で可視化〜時系列データ〜
ggplot2 を使って時系列データの可視化をやってみます。
データは R のベースのデータセットとして用意されている Seatbelts を使います。このデータは1969年から1984年の間で、運転手が重傷または死亡者数やその時のガソリンの値段などが月次で記録された多変量時系列データです。
可視化
データのクラスやどんな変数がいくつあるのかをみてみます。
library(dplyr); library(tidyverse); library(ggplot2) # スタメンライブラリ
ts = datasets::Seatbelts
class(ts); dim(ts); head(ts)# 実行結果
[1] "mts" "ts"
[1] 192 8
DriversKilled drivers front rear kms PetrolPrice VanKilled law
[1,] 107 1687 867 269 9059 0.1029718 12 0
[2,] 97 1508 825 265 7685 0.1023630 6 0
[3,] 102 1507 806 319 9963 0.1020625 12 0
[4,] 87 1385 814 407 10955 0.1008733 8 0
[5,] 119 1632 991 454 11823 0.1010197 10 0
[6,] 106 1511 945 427 12391 0.1005812 13 0ここには表示されていませんが、実はこの ts クラスは日付の情報も含んでる場合が多いです。ggplot2 で可視化する際は日付のカラムが欲しいので、別で作成します。この際に使う lubridate パッケージは、日付に関する操作を大変簡単にしてくれる優れものです。
library(lubridate)
date = seq(ymd('1969/01/01'), by = 'month', length.out = nrow(ts))
data = data.frame(date = date, ts) %>% as_tibble()
data# A tibble: 192 x 9
date DriversKilled drivers front rear kms PetrolPrice VanKilled law
<date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1969-01-01 107 1687 867 269 9059 0.103 12 0
2 1969-02-01 97 1508 825 265 7685 0.102 6 0
3 1969-03-01 102 1507 806 319 9963 0.102 12 0
4 1969-04-01 87 1385 814 407 10955 0.101 8 0
5 1969-05-01 119 1632 991 454 11823 0.101 10 0
6 1969-06-01 106 1511 945 427 12391 0.101 13 0
7 1969-07-01 110 1559 1004 522 13460 0.104 11 0
8 1969-08-01 106 1630 1091 536 14055 0.104 6 0
9 1969-09-01 107 1579 958 405 12106 0.104 10 0
10 1969-10-01 134 1653 850 437 11372 0.103 16 0
# … with 182 more rows1変数の可視化
とりあえず死亡者数をデフォルトでプロットしてみます。
ggplot(data = data, aes(x = date, y = DriversKilled)) +
geom_line()
これを基本として、色々といじっていきます。まず、日付の表示形式を”年/月”にして、表示する日付の数を増やし、文字が被らないように斜めにしてみます。
ggplot(data = data, aes(x = date, y = DriversKilled)) +
geom_line() +
scale_x_date(date_labels = '%Y/%m', date_breaks = '2 year') +
theme(axis.text.x = element_text(angle = 50, hjust = 1))
scale_x_date() の中の date_labels で日付の表示形式を指定し、date_breaks で日付を表示する間隔を指定してます（この場合2年ごと） 。このとき引数に limit = c(as.Date(1969/01/01), as.Date(1975/01/01)) と書けば、表示させる x軸の範囲を変更して局所的に系列をみることができます。theme() の中の element_text() の angle で傾ける角度、hjust で y 軸方向の位置を指定してます。
DriversKilled が大きいところ（>160）だけ図中に日付をプロットし、シートベルト着用が義務付けられた日（1983/01/31）がわかりやすいように縦線を入れてみます。
data_high = data %>% filter(DriversKilled > 160)
ggplot(data = data, aes(x = date, y = DriversKilled)) +
geom_line() +
scale_x_date(date_labels = '%Y/%m/%d', date_breaks = '1 year') +
theme(axis.text.x = element_text(angle = 50, hjust = 1)) +
geom_vline(aes(xintercept = as.Date('1983/02/01')), alpha = 0.7, col = 2, lwd = 1.5) +
geom_text(data = data_high, aes(x = date, y = DriversKilled, label = date))
geom_vline() でどこに縦線を引くかを指定し、geom_text() で図中のどこにラベルを記入するかを指定してます。x 軸がちゃんと対応しているカラムを用意すれば、このように元データからいくつかデータを切り取って、それを図中のラベル表示用に使うことができます。
データに関する印象を述べると、11~12月に死亡事故が多いようです。シートベルト着用義務化の数ヶ月後までは大幅に減少していますが、後の増加具合を見ると「義務化に伴い、それを取り締まるためにパトロールが増えた」といった原因があるのかもしれません。一応義務化以降に過去と比べると水準が下がったように見えますが、義務化後の観測期間が比較的短いのでまだ何とも言い難いです。
月に関する季節成分を丸めてみるために、DriversKilled を年ごとに合計した値をプロットしてみます。点の数が減ったので、geom_point() も足してみやすくしました。
data_year = data %>%
mutate(year = floor_date(`date`, unit = 'year')) %>%
group_by(year) %>%
summarise(DriversKilled_year = sum(DriversKilled))
ggplot(data = data_year, aes(x = year, y = DriversKilled_year)) +
geom_line() +
geom_point() +
scale_x_date(date_labels = '%Y', date_breaks = '1 year')
年でまとめてみると、シートベルト義務化以降の1983年、1984年は劇的に減少してますね。
次に、月毎の特徴を見るため、月別の死者数を棒グラフでみてみます。
data_month = data %>%
mutate(month = month(date, label = TRUE)) %>%
group_by(month) %>%
summarise(DriversKilled_byMonth = sum(DriversKilled))
ggplot(data = data_month, aes(x = month, y = DriversKilled_byMonth)) +
geom_bar(stat = 'identity')
時系列で見た印象通り、12月に近づくにつれて死亡者数が多くなっています。ところが1月になると急に下がっていることがわかります。
理由は分かりませんが、「路面凍結によるスリップ等で事故の数が増える」「日没が早くなるので、視界が悪い時間帯が長くなるため」など、色々と仮説を立てるのには役立ちそうです。
多変量を一気に表示
あまり多くてもごちゃごちゃしてしまうので、今回は front と rear を同時にプロットしてみます。ちなみに、株価のシミュレーションなど、ある系列（確率過程など）がどのような振る舞いをするかなどをみたい場合は多数の系列をいっぺんにプロットしますが、同様の方法でできます。
library(reshape2)
data_melt = data[,c('date', 'front', 'rear')] %>% melt(id = 'date')
ggplot(data = data_melt, aes(x = date, y = value)) +
geom_line(aes(col = variable))
こうみるとわかりやすいですが、前の座席と比べて後部座席の死亡または重症者数はシートベルト着用義務化後もそれほど変わっているようには見えません。これはおそらくシートベルト着用義務が前の座席のみに対するものだったからでしょう。