この文書は
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()
Sunspot.df <- try_data_frame(Sunspot) ggplot(data = Sunspot.df) + aes(x = time, y = Sunspot) + xlab("(b) 太陽黒点数") + ylab("") + geom_line()
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()
BLSALLFOOD.df <- try_data_frame(BLSALLFOOD) ggplot(data = BLSALLFOOD.df) + aes(x = time, y = BLSALLFOOD) + xlab("(d) 食品産業従事者数") + ylab("") + geom_line()
WHARD.df <- try_data_frame(WHARD) ggplot(data = WHARD.df) + aes(x = time, y = WHARD) + xlab("(e) 卸売高") + ylab("") + geom_line()
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()
ggplot(data = Nikkei225) + aes(x = as.numeric(row.names(Nikkei225)), y = nikkei225) + xlab("(g) 日経 225 平均株価") + ylab("") + geom_line()
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("")
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("")
気づきなど
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ファイルを生成してコピペした。画像の添付が面倒である。どうにかならんもんか。