株式会社ホクソエムのブログ

株式会社ホクソエムのブログ

ggplot2 で facet ごとのヒストグラムに平均値の線を引く

ggplot2 で facet ごとのヒストグラムに平均値の線を引きたい。

例えば次のような感じ。

f:id:hoxo_m:20181108192139p:plain

Rコミュニティの Slack である r-wakalang で聞いたところ、即回答がもらえただけでなく、いろいろなやり方を教わったのでメモしておく。 みなさんありがとう。

基本作図

まずはデータを準備する。 ここでは有名な iris データから数値カラム x と facet を作るためのカテゴリカル変数 group を持つデータフレームを作成する。

data(iris) # iris データの読み込み(通常は必要ない)

library(dplyr)

# データの作成
data <- iris %>% select(x = Sepal.Length, group = Species)
# データ内容の確認
head(data)
    x  group
1 5.1 setosa
2 4.9 setosa
3 4.7 setosa
4 4.6 setosa
5 5.0 setosa
6 5.4 setosa

このデータに対して、グループごとにヒストグラムを作成する。

library(ggplot2)
ggplot(data, aes(x)) + 
    geom_histogram() +
    facet_wrap(~group, ncol=1)

f:id:hoxo_m:20181108192206p:plain

この3つのヒストグラムに対して、それぞれのグループの平均値を表す垂直線を追加したいというのが今回の目的。

facet ごとの集計値の扱い

ggplot2 で facet ごとに集計値を扱う場合は、別データとして集計済みのデータを用意する。

例えば今回は平均値の垂直線を引きたいので、次のように平均値を計算したデータフレームを用意する。

data_summary <- data %>%
    group_by(group) %>%
    summarise(mean_x = mean(x))
data_summary
  group      mean_x
1 setosa       5.01
2 versicolor   5.94
3 virginica    6.59

この集計データを用いて平均値をグラフに書き込むことができる(ura_hoxom さん)。

ggplot(data, aes(x)) + 
    geom_histogram() +
    facet_wrap(~group, ncol=1) +
    geom_vline(data = data_summary, aes(xintercept = mean_x), color = "red")

f:id:hoxo_m:20181108192139p:plain

目的は達成された。

他の方法

他の方法も教えてもらった。 ggplot2 では集計データを用意しなくても、集計用の関数を渡すと自動的にデータ集計してくれるらしい。

例えば、次のようにグループごとに平均値を集計する関数を用意する。

mean_by_group = function(d) {
    d %>% group_by(group) %>% summarize(mean_x = mean(x))
}

集計済みデータの代わりに、これを ggplot2 に渡すことができる(heavywatal さん)。

ggplot(data, aes(x)) + 
    geom_histogram() +
    facet_wrap(~group, ncol=1) +
    geom_vline(data = mean_by_group, aes(xintercept = mean_x), color = "red")

f:id:hoxo_m:20181108192139p:plain

また、集計する関数については次のように簡潔に作成できる(yutannihilation さん)*1

mean_by_group = . %>% group_by(group) %>% summarize(mean_x = mean(x))

まとめ

というわけで、みなさんいろいろ教えてくださってありがとうございました。

R-wakalang への参加の仕方については次の資料を参考にしてください。

Enjoy!

*1:これは magrittr の機能らしい。https://github.com/tidyverse/magrittr#building-unary-functions 参照。

ggplot2 で時系列プロットの端点にラベルを表示する

こういう感じで時系列プロットの端点にラベルを表示したい。

f:id:hoxo_m:20180915223055p:plain

この方が時系列とラベルの対応がわかりやすくて良い。

ggrepel パッケージを使うと簡単にできるが、いくつか注意点があるのでここで紹介したい。

まずはデータを用意する。

library(tidyverse)

set.seed(314)
d <- map_dfr(
  c("setosa", "versicolor", "virginica"),
  ~ tibble(
    idx = 1:400,
    value = cumsum(runif(400, -1, 1)),
    type = .
  )
)

ggplot(d, aes(idx, value, colour = type)) + 
  geom_line()

f:id:hoxo_m:20180915223619p:plain

普通にプロットするとこんな感じだが、右端にラベルをつけたい。

ggrepel はいい感じにラベルを配置してくれるパッケージである。

まずはラベルとなるカラムを用意する。 右端にだけラベルをつけたいので ids が最大となる点にだけ type をラベルとして代入し、その他は NA を入れる。

library(ggrepel)

d <- d %>% mutate(label = ifelse(idx == max(idx), type, NA))

ggplot(d, aes(idx, value, colour = type)) +
  geom_line() + 
  geom_label_repel(aes(label = label), na.rm = TRUE)

f:id:hoxo_m:20180915232256p:plain

デフォルトのままだと折れ線にラベルが重なってしまう。 ラベルを右側に出すためには x 軸を右側に広げて nudge_x = Inf を指定すると良い。

ggplot(d, aes(idx, value, colour = type)) +
  geom_line() + 
  geom_label_repel(aes(label = label), na.rm = TRUE, nudge_x = Inf) +
  xlim(NA, 500)

f:id:hoxo_m:20180915232323p:plain

この方法だとラベルはプロットの右端に張り付く。 時系列のすぐ右に出したい場合は、xlim 引数に時系列の端点を指定すると良い。

ggplot(d, aes(idx, value, colour = type)) +
  geom_line() + 
  geom_label_repel(aes(label = label), na.rm = TRUE, xlim = c(400, NA)) + 
  xlim(NA, 500)

f:id:hoxo_m:20180915232359p:plain

謝辞

上記は r-wakalang の Slack 上で @u_ribo さんと @yutannihilation さんに教えていただきました。 ありがとうございます。

参考

RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−

RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−

データ分析のワークフローをdrakeで管理して効率的に作業を進めよう

要約

  • drakeパッケージは、GNU makeのようにあらかじめ定義されたワークフローを自動的に実施する仕組みを、Rユーザに馴染みやすいデータフレーム形式で提供する
  • ワークフローの構築と管理、実行はRの関数として提供され、依存関係を可視化する関数も用意される
    • drakeパッケージを使うことで、データ分析でありがちな「再実行」の負担(再計算、コードの保守)を軽減することが可能となる
    • 各オブジェクトは自動的にキャッシュされ、コードや依存関係に変更のない場合はキャッシュが利用される
    • ワークフローの各処理の状況、依存関係を可視化する関数も用意され、ワークフロー管理が容易になる

f:id:u_ribo:20180905064443j:plain

  • 要約
  • はじめに
    • シーシュポスの岩
    • 既存の解決策
  • drake: Rユーザのためのワークフロー処理パッケージ
    • ワークフロー管理の基礎
    • ワークフローと依存関係の可視化
    • ワークフローの変更
  • 参考URL

はじめに

データ分析の作業は、試行錯誤を繰り返して進めて行くのが一般的です。

Garrett GrolemundとHadley Wickhamの著書 “R for Data Science” (「Rではじめるデータサイエンス」)では、Rを使ったデータ分析作業の工程を以下の図のように整理しています (FontAwesomeのフォントを使って編集しています)。

f:id:u_ribo:20180831061703p:plain

この図の中で、データ整形と可視化、そしてモデリングの工程がひとまとまりになっており、これらの工程を経て分析結果を伝えるための作業に入ります。

データ整形と可視化、モデリングは互いに影響しあう関係にあり、繰り返し実行されることを想定しています。繰り返しの例としてはモデルに用いる変数のアップデートがあります。その際、モデルを実行するコードを修正するだけでなく、それに関連する前後のプロセス、すなわちモデルデータの生成とモデルの結果を利用するプログラムについても見直す必要があるでしょう。

こうした「繰り返し」や「やり直し」の作業は、みなさんご存知の通り、時間がかかる骨の折れる仕事で、退屈なものです。場合によっては、ほとんど1からコードを書き直すなんて事もあるかもしれません。間に実行時間の長い処理が途中に入る時も、処理が終わるまで待つ必要があります。また、変更箇所に見落としがあるとフローは流れなくなってしまいます。変更が多ければ多いほど、人為的な誤りも犯しやすくなる危険があります。

続きを読む

【2018年版】R でハッシュテーブルの速度比較

以前、こういう記事を書いた。

Rでハッシュテーブルを使う方法はいくつかあるが、

  • サイズが 1000 以下ならば名前付きベクトルが速い
  • それ以上なら環境を使った方法が速い
  • hash パッケージは遅いが記述がわかりやすい

ということがわかった。

今回はさらに高速なハッシュテーブルを実装するパッケージを3つ見つけたのでこれらを加えて比較してみる。

使用するデータは前回と同じく次のコードで生成した。

generate_random_labels <- function(num, length = 10) {
  apply(Vectorize(sample, "size")(letters, rep(num, length), replace = TRUE), 1, function(row) paste(row, collapse = ""))
}
generate_random_values <- function(labels, digits = 1) {
  format <- paste0("%s%0", digits, "d")
  sprintf(format, labels, Vectorize(sample, "size")(seq_len(10^digits-1), rep(1,length(labels))))
}

N <- 1000
labels <- generate_random_labels(N)
values <- generate_random_values(labels)
df <- data.frame(labels, values, stringsAsFactors = FALSE)

target <- sample(labels, 5)
target
#> [1] "eguiqhmukh" "yuvlnhfvba" "cjbcajffur" "nttsfcuszh" "ibvgfrjnro"

1. 名前付きベクトル

名前付きベクトルは最もシンプルなハッシュテーブルの実装である。

named_values <- values
names(named_values) <-  labels
named_values[target]
#>    eguiqhmukh    yuvlnhfvba    cjbcajffur    nttsfcuszh    ibvgfrjnro 
#> "eguiqhmukh2" "yuvlnhfvba2" "cjbcajffur1" "nttsfcuszh8" "ibvgfrjnro6" 

2. 環境

R言語徹底解説』第8章 (p.160) には、R における環境 (environment) の特別な使用方法としてハッシュマップとして使えると書いてある。 前回の記事ではやや複雑に書いたが、list2env(setNames()) イディオムを使うとやや簡潔に書ける。

hash_env <- list2env(setNames(as.list(values), labels), hash = TRUE)
unlist(mget(target, envir = hash_env))
#>    eguiqhmukh    yuvlnhfvba    cjbcajffur    nttsfcuszh    ibvgfrjnro 
#> "eguiqhmukh2" "yuvlnhfvba2" "cjbcajffur1" "nttsfcuszh8" "ibvgfrjnro6" 

3. hash パッケージ

hash パッケージを使うと直感的で簡潔な記法でハッシュテーブルを扱える。 結果の順序が保証されないことには注意が必要である。

library(hash)
h <- hash(labels, values)
values(h[target])
#>    cjbcajffur    eguiqhmukh    ibvgfrjnro    nttsfcuszh    yuvlnhfvba 
#> "cjbcajffur1" "eguiqhmukh2" "ibvgfrjnro6" "nttsfcuszh8" "yuvlnhfvba2" 

4. collections パッケージ

collections は高速なコンテナデータ型を提供するパッケージであり、ハッシュテーブルは辞書型 Dict として提供される。 メソッドはベクトル化されていないことに注意が必要である。

library(collections)
d <- Dict$new()
invisible(Vectorize(d$set)(labels, values))
Vectorize(d$get)(target)
#>    eguiqhmukh    yuvlnhfvba    cjbcajffur    nttsfcuszh    ibvgfrjnro 
#> "eguiqhmukh2" "yuvlnhfvba2" "cjbcajffur1" "nttsfcuszh8" "ibvgfrjnro6"

5. datastructures

datastructures パッケージはいくつかの高度なデータ型を提供する。 ハッシュテーブルは hashmap で提供される。

library(datastructures)
hm <- hashmap("character")
hm[labels] <- values
unlist(hm[target])
#> [1] "eguiqhmukh2" "yuvlnhfvba2" "cjbcajffur1" "nttsfcuszh8" "ibvgfrjnro6"

6. hashmap パッケージ

hashmap パッケージはキーとバリューに特定の型しか使えないが環境より速いとのこと。

library(hashmap)
H <- hashmap(labels, values)
H[[target]]
#> [1] "eguiqhmukh2" "yuvlnhfvba2" "cjbcajffur1" "nttsfcuszh8" "ibvgfrjnro6"

速度比較

上で紹介した 6 つの方法について速度を比較してみる。 テーブルサイズ N を 100~100000 まで変化させて速度を調べている。

library(magicfor)

magic_for()

for(N in 10^(2:5)) {
  # Create Data -------------------------------------------------------------
  labels <- generate_random_labels(N)
  values <- generate_random_values(labels)
  df <- data.frame(labels, values, stringsAsFactors = FALSE)
  
  target <- sample(labels, 100)
  
  # Preparation -------------------------------------------------------------
  named_values <- values
  names(named_values) <-  labels
  
  hash_env <- list2env(setNames(as.list(values), labels), hash = TRUE)
  
  library(hash)
  h <- hash(labels, values)
  
  library(collections)
  d <- Dict$new()
  invisible(Vectorize(d$set)(labels, values))
  
  library(datastructures)
  hm <- datastructures::hashmap("character")
  hm[labels] <- values
  
  library(hashmap)
  HH <- hashmap::hashmap(labels, values)
  
  # Benchmark ---------------------------------------------------------------
  library(microbenchmark)
  mb_result <- microbenchmark(
    named_vector = {
      named_values[target]
    }, 
    environment = {
      unlist(mget(target, envir = hash_env))
    },
    hash = {
      hash::values(h[target])
    },
    collections = {
      Vectorize(d$get)(target)
    },
    datastructure = {
      unlist(hm[target])
    },
    hashmap = {
      HH[[target]]
    }
  )
  
  time_df <- data.frame(N=as.character(N), mb_result)
  time_df
}

result <- Reduce(rbind, magic_result()$time_df)

結果を可視化する。

library(ggplot2)
ggplot(result, aes(x=N, y=time, fill=expr)) + 
  geom_boxplot() + scale_y_log10() + facet_wrap(~expr)

f:id:hoxo_m:20180828234546p:plain

library(dplyr)
result_mean <- result %>% 
  group_by(N, expr) %>% 
  summarise(mean=mean(time))
ggplot(result_mean, aes(x=N, y=mean, group=expr, color=expr)) + 
  geom_line() + scale_y_log10()

f:id:hoxo_m:20180828234555p:plain

テーブルサイズが数百程度であれば名前付きベクトルが速い。 それ以上は環境が一番速いようである。 hashmap は環境と遜色ない速度であり、直感的で簡潔な記述ができる。 hash パッケージを使う理由はもはや無いようである。

R言語徹底解説

R言語徹底解説

モデルで扱うデータの前処理をrecipesで行う

ドーモ。ホクソエムの @u_ribo です。本業ではモデリングとは離れたギョームをしています。寂しくなったので、Rのrecipesパッケージについて紹介します。

tidymodels.github.io

モデルに適用するデータの前処理

Rでのモデル式 (model formula) の記述って、利用時に不便を感じることや覚えるのが難しい面が時々ありませんか?例えば、y ~ .」は右辺のドットが、目的変数以外の全ての変数を説明変数として扱うことを示しますが、説明変数に対数変換などの変数変換を行うにはy ~ log(.)という記述はできず、結局、説明変数を「+」でつなげていくことになります。また、交互作用項の指定には「x1 * x2」や「(x1 + x2)^2」、「:」を使う表記が可能ですが、この表記には最初は混乱しませんか?(単に私が不勉強なだけということもあります)

加えて、多くのモデルでは欠損への処理が必要となったり、適用するモデルによって(例えば複数の説明変数を扱うK-NNやSVM)は、変数間の標準化が必要です。

ここで紹介するrecipesパッケージを使うことで、こうしたモデル構築に伴うデータの前処理を楽にすることが期待されます。

  • モデルに適用するデータの前処理
  • 対数変換を行うモデルを例に
    • recipesパッケージによるデータの前処理
  • 交互作用を扱うモデルを例に
  • stepに指定可能な処理
  • まとめ
続きを読む

クロネッカー積でデータを列方向(or行方向)に高速に複製もしくは定数倍する

r-wakalangからの転載です。以下のような質問がありました。

data.frameをカラム・ロウ方向に複製結合したdata.frameを出力させたいのですが、どうも綺麗に書けずです。。アドバイスお願いしますm( )m

この意味は、例えば「行方向に2個・列方向に3個複製」の場合は以下の図です。

f:id:StatModeling:20180824174855p:plain

データフレームが数値のみの場合、これは実はクロネッカー積と呼ばれる行列の特殊な積で表現できます。

  • Rで書く場合

Rではデフォルトで%x%関数(or kronecker関数)が用意されているのでそれを使えば簡単に実装できます。

B <- as.matrix(data.frame(a=1:3, b=4:6))
A <- matrix(1, nrow=2, ncol=3)
res <- A %x% B
  • Pythonで書く場合

Pythonでもnumpykronメソッドが用意されていますので、以下のように書くことができます。

import pandas as pd
import numpy as np

B = pd.DataFrame(np.arange(1,7).reshape(2,3).T, columns=list('ab'))
A = pd.DataFrame(np.ones((2,3)))

kp = pd.DataFrame(np.kron(A,B), columns=pd.MultiIndex.from_product([A,B]))

  * * *  

定数倍も行けます。例えば、以下のように

f:id:StatModeling:20180824163606p:plain

左のデータフレームから右のデータフレームが得たい場合です。

  • Rで書く場合
A <- as.matrix(data.frame(a=1:3, b=4:6))
B <- matrix(1:50, nrow=1)
res <- A %x% B
  • Pythonで書く場合
import pandas as pd
import numpy as np

A = pd.DataFrame(np.arange(1,7).reshape(2,3).T, columns=list('ab'))
B = pd.DataFrame(np.arange(1,51).reshape(1,50))

kp = pd.DataFrame(np.kron(A,B), columns=pd.MultiIndex.from_product([A,B]))

Enjoy!

ggplot2 で時系列の区間に影をつけるのは annotate が便利ぽい

例えば次のような時系列データがあるとします。

library(xts)
ts <- as.xts(Nile)

library(ggplot2)
autoplot(ts)

f:id:hoxo_m:20180821233320p:plain

このプロットの 1900年から1940年までの区間に影をつけたい。

これには annotate() が便利ぽい。

autoplot(ts) + 
  annotate("rect", 
           xmin = as.Date("1900-01-01"), 
           xmax = as.Date("1940-01-01"), 
           ymin = -Inf, ymax = Inf, alpha = 0.2)

f:id:hoxo_m:20180821233547p:plain

Enjoy!

参考