tidy時系列データにおける相関計算 corrr
時系列データに対して相関を出す場面で、毎回どういう変換するんだっけを調べている気がするためメモします。
corrrパッケージでどのようなことができるのかについては、kazutanさんのページが大変参考になります。
kazutan.github.io
データ
店舗ごとの何らかの商品の売り上げのようなデータがあり、各店舗間の売り上げが相関しているのか分析したいということを想定しています。
library(tidyverse) library(lubridate) library(corrr) set.seed(12345) store_sales_tbl <- tibble( date = rep(seq(ymd("2019-01-01"), ymd("2019-12-31"), by = "1 week"), 5), store_name = c(rep("A", 53), rep("B", 53), rep("C", 53), rep("D", 53), rep("E", 53)), sales = rpois(53 * 5, 10) )
データの形式としては、SQLでそのままの抽出してきたような状態が想定されます。
> store_sales_tbl # A tibble: 265 x 3 date store_name sales <date> <chr> <int> 1 2019-01-01 A 11 2 2019-01-08 A 12 3 2019-01-15 A 9 4 2019-01-22 A 8 5 2019-01-29 A 11 6 2019-02-05 A 4 7 2019-02-12 A 11 8 2019-02-19 A 7 9 2019-02-26 A 8 10 2019-03-05 A 11 # … with 255 more rows
相関計算
このようなデータに対して、tidyverse的にstats::corを適用するには意外と面倒です。corrr::correlateが適用できる形に持っていくのが早いです。
ここでは、時系列データの相関について、細かい考慮することはしませんが、変化率を計算する関数などを用意しておきます（ファイナンス関連だと対数変化率などの関数を利用することを想定します）。
cr <- function(x) { (x - lag(x)) / lag(x) * 100 }
あとは、このデータを一度、横持ちのデータに変換してから、変化率の系列に置き換え、corrr::correlateに繋げれば相関行列を算出できます。
store_sales_cor_df <- store_sales_tbl %>% spread(store_name, sales) %>% transmute_if(is.integer, cr) %>% correlate()
算出された相関行列を確認します。
> store_sales_cor_df # A tibble: 5 x 6 rowname A B C D E <chr> <dbl> <dbl> <dbl> <dbl> <dbl> 1 A NA 0.118 -0.0512 -0.102 0.241 2 B 0.118 NA -0.0853 0.0214 0.128 3 C -0.0512 -0.0853 NA -0.304 -0.212 4 D -0.102 0.0214 -0.304 NA 0.158 5 E 0.241 0.128 -0.212 0.158 NA
可視化
ここまでくると、ggplot2で可視化することもできますが、corrr::rplotを使えば、簡単に相関行列を可視化することができます。
store_sales_cor_df %>% rplot(print_cor = TRUE)
corrplotの方が使い慣れているという場合には、corrr::as_matrixを利用すれば可視化まで持っていけます。
library(corrplot) store_sales_cor_df %>% as_matrix() %>% corrplot.mixed(lower = "number", upper = "square")
lubridate::floor_dateで任意の曜日スタートで切り下げ
タイトル通りのメモ。
library(tidyverse) library(lubridate) set.seed(12345) order_tbl <- tibble( order_date = seq(ymd("2019-03-01"), ymd("2019-03-31"), by = "1 day"), order_num = rpois(31, 5) )
日曜日と月曜日スタートする場合を追加します。
> order_tbl %>% + mutate(wd = weekdays(order_date), + sunday_start = floor_date(order_date, unit = "week"), + monday_start = floor_date(order_date, unit = "week", week_start = 1)) # A tibble: 31 x 5 order_date order_num wd sunday_start monday_start <date> <int> <chr> <date> <date> 1 2019-03-01 6 金曜日 2019-02-24 2019-02-25 2 2019-03-02 8 土曜日 2019-02-24 2019-02-25 3 2019-03-03 6 日曜日 2019-03-03 2019-02-25 4 2019-03-04 8 月曜日 2019-03-03 2019-03-04 5 2019-03-05 5 火曜日 2019-03-03 2019-03-04 6 2019-03-06 3 水曜日 2019-03-03 2019-03-04 7 2019-03-07 4 木曜日 2019-03-03 2019-03-04 8 2019-03-08 5 金曜日 2019-03-03 2019-03-04 9 2019-03-09 6 土曜日 2019-03-03 2019-03-04 10 2019-03-10 11 日曜日 2019-03-10 2019-03-04 # … with 21 more rows
つまり、week_start = 1を指定すれば、月曜日スタートで切り下げられます。デフォルトでは、week_start = getOption("lubridate.week.start", 7)が設定されていますので、日曜日スタートになります。
round_dateとcelling_dateを含むfloor_dateの実装場所は以下になります。
github.com
lubridate::rollbackの使い所
日次データが存在して、その日次データとそれを月ごとにまとめ上げた月次データの可視化を行いたい場合があります。
library(tidyverse) library(lubridate) set.seed(12345) order_tbl <- tibble( order_date = seq(ymd("2018-01-01"), ymd("2018-12-31"), by = "1 day"), order_num = rpois(365, 5) )
日次データとしてプロットしたい場合は、以下のようにできます。
order_tbl %>% ggplot(aes(x = order_date, y = order_num)) + geom_point(size = 2, colour = "blue") + geom_line(size = 1, colour = "blue", alpha = 0.5) + scale_x_date(date_breaks = "3 months", date_minor_breaks = "1 months", date_labels = "%Y-%m") + scale_y_continuous(breaks = seq(0, 20, by = 2), minor_breaks = 1:20) + theme(axis.text = element_text(size = 15), axis.title = element_text(size = 20))
次に、月ごとにまとめ上げたデータは以下のように計算できます (年と月のまとめ上げでも可能）。
order_tbl %>% mutate(y = year(order_date), m = month(order_date), ym = str_c(y, "-", sprintf("%02d", m))) %>% group_by(ym) %>% summarise(monthly_order_sum = sum(order_num)) #> # A tibble: 12 x 2 #> ym monthly_order_sum #> <chr> <int> #> 1 2018-01 161 #> 2 2018-02 133 #> 3 2018-03 161 #> 4 2018-04 163 #> 5 2018-05 176 #> 6 2018-06 171 #> 7 2018-07 156 #> 8 2018-08 169 #> 9 2018-09 161 #> 10 2018-10 144 #> 11 2018-11 146 #> 12 2018-12 170
しかし、これだと年月を表すカラムがcharactor型のため、上記の日次データを可視化した際に利用したscale_x_dateの部分を流用して同じスケールで可視化できません。
そこで、lubridate::rollbackを利用すれば、以下のようにできます。
order_tbl %>% mutate(first_date = rollback(order_date, roll_to_first = TRUE)) %>% group_by(first_date) %>% summarise(monthly_order_sum = sum(order_num)) %>% ggplot(aes(x = first_date, y = monthly_order_sum)) + geom_point(size = 2, colour = "blue") + geom_line(size = 1, colour = "blue", alpha = 0.5) + scale_x_date(date_breaks = "3 months", date_minor_breaks = "1 months", date_labels = "%Y-%m") + scale_y_continuous(breaks = seq(0, 200, by = 10), minor_breaks = seq(0, 200, by = 5)) + theme(axis.text = element_text(size = 15), axis.title = element_text(size = 20))
ここでは、その月の最初の日のカラムを追加して、その日で月次にまとめ上げて可視化しています。こうすることで、日次データで利用したscale_x_dateの部分も流用できます。
rollback関数の中身自体は、day関数を使って、日付部分に1を入れるというシンプルな操作を行う便利関数になっています。
github.com
Rでモデルの評価指標を計算する Metrics
機械学習モデルの分類精度を評価したり、時系列モデルのn期先の予測精度を評価したりする際に、地味に評価指標の計算式をネットで調べて自前でRの関数として実装したりしていました。
しかし、（実際にはよく調べていなかっただけですが、）すでに基本的な指標はパッケージとしてまとめられていました。パッケージ名は、至ってそのままでMetricsパッケージになります。
F1、Rrecision、Recallをはじめ、MAPEやRMSEなど基本的な評価指標は用意されています。
実装はシンプルで、関数に入力する実際の値と予測した値の個数が異なる場合でも、警告メッセージは出ますが値は出力される感じなので、その点は少し注意が必要です。
rvestでUser Agentを偽装する
「偽装する」というと感じが悪いですが、いつも使っているChromeやSafariなどのブラウザで表示される同じ状態をスクレイプしたかったということで、rvestでのユーザーエージェントの変更の仕方を調べたメモ。
結論としては、rvest::html_session(url, ...)の...の部分にセットすれば、httr::GET(url = NULL, config = list(), ..., handle = NULL)のconfigにセットするようになっていました。
library(tidyverse) library(rvest) library(httr) URL <- "http://had.co.nz" UA <- "Mozilla/5.0 (Macintosh; Intel Mac OS X 10_14_2) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/71.0.3578.98 Safari/537.36" session <- rvest::html_session(URL, httr::user_agent(UA))
Status CodeとUser Agentを確認します。
> session <session> http://hadley.nz/ Status: 200 Type: text/html Size: 8981 > session$config <request> Options: * useragent: Mozilla/5.0 (Macintosh; Intel Mac OS X 10_14_2) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/71.0.3578.98 Safari/537.36 * autoreferer: 1 > session %>% + html_node(xpath = '//*[@id="page-top"]/footer/div/div/div/p/a[2]') %>% + html_text(trim = TRUE) [1] "@hadleywickham"
当然ながら、User Agentだけではなく、httr::set_cookiesを利用すれば、Cookieもセットすることができます。
COOKIE1 <- "www" COOKIE2 <- "zzz" session <- rvest::html_session(URL, httr::user_agent(UA), httr::set_cookies(w = COOKIE1, z = COOKIE2))
以下のように、Cookie部分が追加されていることが確認できます。
> session$config <request> Options: * useragent: Mozilla/5.0 (Macintosh; Intel Mac OS X 10_14_2) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/71.0.3578.98 Safari/537.36 * cookie: w=www;z=zzz * autoreferer: 1
html_sessionの実装場所は以下になります。
github.com
ggplot2で値の正負によって、色分けする
データの値の正負によって色分けをしたい場合があります。この場合、mutateでTRUE/FALSEの判定をした新しいカラムを追加して可視化するということをよくしていましたが、データが大変多い場合などには、可視化の色分けのためだけにカラムを追加するということはあまりしたくありません。
ggplotの中でなんとかできないか調べていたら、aesの中で比較演算子を使って簡単に色分けができる。ということを今さら知ったので、メモします。
library(tidyverse) sample_tbl <- tibble( type = LETTERS, value = runif(length(LETTERS), -1, 1) ) sample_tbl %>% ggplot(aes(x = type, y = value, fill = value > 0)) + geom_bar(stat = "identity") + labs(x = "", y = "") + theme_bw() + theme(axis.text = element_text(size = 20), legend.position = "none")
value > 0を指定するだけでした...
加えて、こちらはよく調べられていると思われる昇順、降順で並べ替えてのプロットですが、reorderを使ってすぐにできます。
sample_tbl %>% ggplot(aes(x = reorder(type, value), y = value, fill = value > 0)) + geom_bar(stat = "identity") + labs(x = "", y = "") + theme_bw() + theme(axis.text = element_text(size = 20), legend.position = "none")
lubridateのweekとisoweekの違い
日次データを週次データとして売り上げなどをまとめて集計したいということはよくあります。今まで特に何も気にすることなく、lubridateのisoweekを使って週番号を付けて集計していましたが、年を跨いだ場合にどのような週番号が付くのかを確認したのでメモします。
ここでは、年を跨いだ16日分の日次データがあったとします。
> library(tidyverse) > library(lubridate) > sample_tbl <- tibble( + date = seq(as.Date("2018-12-23"), as.Date("2019-01-07"), by = 1), + sales = sample(1:100, 16) + ) > sample_tbl # A tibble: 16 x 2 date sales <date> <int> 1 2018-12-23 19 2 2018-12-24 68 3 2018-12-25 37 4 2018-12-26 36 5 2018-12-27 84 6 2018-12-28 86 7 2018-12-29 59 8 2018-12-30 13 9 2018-12-31 72 10 2019-01-01 40 11 2019-01-02 96 12 2019-01-03 69 13 2019-01-04 23 14 2019-01-05 28 15 2019-01-06 6 16 2019-01-07 4
これに対して、それぞれの週番号と曜日を追加します。
> sample_tbl %>% + mutate(week = week(date)) %>% + mutate(isoweek = isoweek(date)) %>% + mutate(epiweek = epiweek(date)) %>% + mutate(weekdays = weekdays(date)) # A tibble: 16 x 6 date sales week isoweek epiweek weekdays <date> <int> <dbl> <dbl> <dbl> <chr> 1 2018-12-23 73 51 51 52 日曜日 2 2018-12-24 87 52 52 52 月曜日 3 2018-12-25 75 52 52 52 火曜日 4 2018-12-26 86 52 52 52 水曜日 5 2018-12-27 44 52 52 52 木曜日 6 2018-12-28 16 52 52 52 金曜日 7 2018-12-29 31 52 52 52 土曜日 8 2018-12-30 48 52 52 1 日曜日 9 2018-12-31 67 53 1 1 月曜日 10 2019-01-01 91 1 1 1 火曜日 11 2019-01-02 4 1 1 1 水曜日 12 2019-01-03 14 1 1 1 木曜日 13 2019-01-04 65 1 1 1 金曜日 14 2019-01-05 1 1 1 1 土曜日 15 2019-01-06 34 1 1 2 日曜日 16 2019-01-07 40 1 2 2 月曜日
つまり、weekは1月1日からスタートで週番号を付けるので、53週目が存在します。一方、isoweekはISO 8601規格で「年の最初の木曜日を含む週が、その年の第1週」として週番号を付けるので、2018年12月31日も1週目となります（2017年などは年を越しても1月1日は52週目となります）。
大抵のビジネスシーンでは、月曜始まりの7日間の売り上げを週次データとして集計するので、isoweekを使うことになると思います。
おまけとして、epiweekというものもあり、こちらは日曜日を週の始まりとした週番号になります。
