ブログ

κυσο βλογ

Rによる時系列モデリング入門 第1章 データのグラフ化まで

この文書は

2020年12月11日に岩波書店から出版された、北川源四郎著の『Rによる時系列モデリング入門』について、最初のデータの読み込みから基本的なグラフ化までを確認のために出力してみたものである。やってみて、いくつか気づいた点があるのでメモを残す。

必要なパッケージをロード

本書で必要な関数とデータを全て含むTSSSパッケージをロードする。

library(TSSS)

以下のパッケージはRを用いた現代的なデータ解析には必要不可欠である。

library(tidyverse)

以下のパッケージはtsクラスの変数をデータフレームに変換するために用いる。

library(ggpmisc)

データ取り込み

data("HAKUSAN")
data("Sunspot")
data("Temperature")
data("BLSALLFOOD")
data("WHARD")
data("MYE1F")

以下の3つのデータセットはTSSSパッケージに含まれていないようである。

  • Nikkei225
  • rainfall
  • Haibara

webで検索したら著者が公開しているwebページが見つかった。 恐らく以下のこれらが最新版の元データであろう。

Nikkei225 <- read_csv("http://www.mi.u-tokyo.ac.jp/mds-oudan/lecture_document_2019_math7/%E6%99%82%E7%B3%BB%E5%88%97%E3%83%87%E3%83%BC%E3%82%BF/nikkei225_new.csv")

Rainfall <- read_csv("http://www.mi.u-tokyo.ac.jp/mds-oudan/lecture_document_2019_math7/%E6%99%82%E7%B3%BB%E5%88%97%E3%83%87%E3%83%BC%E3%82%BF/rainfall_new.csv")

Haibara <- read_csv("http://www.mi.u-tokyo.ac.jp/mds-oudan/lecture_document_2019_math7/%E6%99%82%E7%B3%BB%E5%88%97%E3%83%87%E3%83%BC%E3%82%BF/haibara_new.csv")

ただし、これらのURLがいつまで有効かは不明である。

CSVデータの取り込み例は省略する。上のread_csvによるwebからのデータ読み込みで例示は十分だろう。

plot例

本書では素朴なplot関数によるグラフ化が説明されているが、より実践的なggplotを用いてグラフ化を行う。

ggplot(data = HAKUSAN) +
  aes(x = as.numeric(row.names(HAKUSAN)), y = YawRate) +
  xlab("(a) 船舶の方向角速度") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002828p:plain

Sunspot.df <- try_data_frame(Sunspot)
ggplot(data = Sunspot.df) +
  aes(x = time, y = Sunspot) +
  xlab("(b) 太陽黒点数") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002831p:plain

Temperature.df <- data.frame(number = Temperature)
ggplot(data = Temperature.df) +
  aes(x = as.numeric(row.names(Temperature.df)), y = number) +
  xlab("(c) 東京の日最高気温") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002834p:plain

BLSALLFOOD.df <- try_data_frame(BLSALLFOOD)
ggplot(data = BLSALLFOOD.df) +
  aes(x = time, y = BLSALLFOOD) +
  xlab("(d) 食品産業従事者数") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002837p:plain

WHARD.df <- try_data_frame(WHARD)
ggplot(data = WHARD.df) +
  aes(x = time, y = WHARD) +
  xlab("(e) 卸売高") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002841p:plain

MYE1F.df <- data.frame(number = MYE1F)
ggplot(data = MYE1F.df) +
  aes(x = as.numeric(row.names(MYE1F.df)), y = number) +
  xlab("(f) 地震波(東西方向)") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002845p:plain

ggplot(data = Nikkei225) +
  aes(x = as.numeric(row.names(Nikkei225)), y = nikkei225) +
  xlab("(g) 日経 225 平均株価") + ylab("") +
  geom_line()

f:id:XX-Prime:20201215002849p:plain

Haibara.tidy <- Haibara %>% 
  gather(key = measure, value = value, 地下水位, 気圧)
ggplot(Haibara.tidy, aes(x = as.numeric(日時), y=value)) + 
  geom_line() + facet_grid(measure ~ ., scales = "free") + 
  xlab("(h) 地下水位(下段)と気圧(上段)") + ylab("")

f:id:XX-Prime:20201215002852p:plain

HAKUSAN.tidy <- HAKUSAN %>% 
  mutate(index = as.numeric(row.names(.))) %>% 
  gather(key = measure, value = value, Rolling, Pitching, Rudder)
ggplot(HAKUSAN.tidy, aes(x = index, y=value)) + 
  geom_line() + facet_grid(measure ~ ., scales = "free") + 
  xlab("(i) 船舶の縦揺れ(上段), 横揺れ(中段), 舵角(下段)") + ylab("")

f:id:XX-Prime:20201215002855p:plain

気づきなど

tsクラスのデータはそのままではggplotでグラフ化できない

(b)(d)(e)をグラフ化する際に、ggpmiscのtry_data_frame()でデータフレームに変換してからggplotでグラフ化した。

参考

https://www.karada-good.net/analyticsr/r-453

facetで順序を制御するのはどうやら面倒

(h)(i)をグラフ化する際に、もとのデータフレームをtidyな形式に変換した後にfacet_gridで縦に並べて表示したが、順番はアルファベット順になってしまうようで本の順番通りには配置できなかった。ファクター変数でやや無理やり順番を制御することもできるようだが、面倒なのでggplotの表示順をそのまま採用した。

(h)地下水位のデータは長い欠損期間があるようだ

本では行数をx軸のインデックスとして使っているようだが、地下水位データには「日時」という列があり、今回はこれをx軸として図示した所、図のように長い欠損期間が二ヶ所あることが分かった。

本文ではこの長い欠損期間については特に触れられていないようである。

この文書をどうやって作成したか

rmarkdownで作成したが、出力された文書をHatena Blogに載せる方法がわからなかった。以下の解説を参考にmdファイルを生成してコピペした。画像の添付が面倒である。どうにかならんもんか。

参考

mfavoritey.hatenablog.com