こんにちは、買物情報事業部でサーバサイドの開発を担当している荒引 (@a_bicky) です。
今回のエントリでは R で A/B テストの結果検証を行う方法の一例について紹介します。 エンジニアでも自分の関わった施策の効果検証のために簡単な分析をすることがあるかと思いますが、そんな時にこのエントリが役立てば幸いです。
なお、次のような方は対象外です。
- A/B テストや KPI の設計に興味のある方
- この辺には全く触れません
- プログラミング初心者
- わからない単語が大量に出てくるでしょう
- R で統計学や機械学習の手法をバリバリ使いたい方
- 世の中の “分析” の多くは集計処理がメインです
- Python, Julia など既に分析する上で使い慣れた言語・ツールがある方
- 今回のエントリ程度の内容であればわざわざ乗り換える必要もないでしょう
OS は Mac を前提として説明するので、Windows や Linux では一部動かないものもあるかもしれませんがご了承ください。 また、R のバージョンは現時点で最新バージョンの 3.2.0 であることを前提とします。
何故 R か？
それは私の一番使える言語だからです！というのが一番の理由ですが、他にも次のような理由が挙げられます。
- 無料で使える
- R 関連の書籍なども大量に存在していて情報が豊富
- RStudio や ESS のような素晴らしい IDE が存在する
- パッケージ（Ruby でいう gem）が豊富
- ggplot2 パッケージを使うことで複雑なグラフが手軽に描ける
- data.table, dplyr, stringi パッケージなどのおかげで数百万オーダーのデータでもストレスなく高速に処理できるようになった
ちなみに、R のコーディングスタイルのカオスっぷりに辟易する方もいるでしょうが、そこは耐えるしかないです。*1
アジェンダ
- R の環境構築
- 最低限の設定
- IDE の導入
- R の使い方についてさらっと
- コンソールの使い方
- デバッグ方法
- data.frame について
- A/B テストの結果検証
- コンバージョン率 (CVR) を出す
- CVR の差の検定をする
- ユーザの定着率（定着数）を出す
- 番外編
- R でコホート分析
- 最後に
R の環境構築
Mac の場合は Homebrew 経由でインストールするのが手軽です。Linux の場合もパッケージマネージャ経由で手軽にインストールできます。
% brew tap homebrew/science
% brew install r
これでターミナルから R を起動できるようになったはずです。
% R R version 3.2.0 (2015-04-16) -- "Full of Ingredients" Copyright (C) 2015 The R Foundation for Statistical Computing Platform: x86_64-apple-darwin14.3.0 (64-bit) R は、自由なソフトウェアであり、「完全に無保証」です。 一定の条件に従えば、自由にこれを再配布することができます。 配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。 R は多くの貢献者による共同プロジェクトです。 詳しくは 'contributors()' と入力してください。 また、R や R のパッケージを出版物で引用する際の形式については 'citation()' と入力してください。 'demo()' と入力すればデモをみることができます。 'help()' とすればオンラインヘルプが出ます。 'help.start()' で HTML ブラウザによるヘルプがみられます。 'q()' と入力すれば R を終了します。 >
R を使う際には基本的にこのコンソール上で作業することになります。コマンドラインからスクリプトを実行するのはバッチなどを運用するケースぐらいです。
R を終了するには q 関数を実行するか Ctrl-D を入力します。
> q() Save workspace image? [y/n/c]: n
“Save workspace image?” で y を入力すると現在のセッションの内容がワーキングディレクトリに .RData という名前で保存され、次回起動時に定義したオブジェクトなどが復元されます。
最低限の設定
個人的に次の 2 つの設定だけは欠かせないと思っています。
- 言語の英語化
- エラーメッセージでググった時に参考になる情報が多いので、英語アレルギーでない限りは英語化した方が良いでしょう
- リポジトリの設定
- パッケージをインストールする際にいちいち指定する必要がなくなります
次のコマンドを実行することで上記の設定が可能です。
% # 言語の英語化 % echo 'LANGUAGE=en' >> ~/.Renviron % # パッケージをインストールする際のリポジトリを設定 % echo 'options(repos = "http://cran.ism.ac.jp")' >> ~/.Rprofile
.Renviron は R を実行する上での環境変数を定義するファイルです。
.Rprofile は起動時の処理を記述するファイルです。options 関数の repos 引数に統計数理研究所のリポジトリを指定しています。*2
IDE の導入
前述したように、R を使う際には基本的にコンソール上で作業することになりますが、それなりの処理を行う場合は IDE を利用すると便利です。 IDE を利用することで、スクリプト上のコードをコンソールに送って挙動を確認するようなことも可能になります。 メジャーな IDE は次の 2 つです。
RStudio を使うには、ダウンロードリンクから自分の環境にあったファイルをダウンロードしてインストールするだけです。 Option+Shift+K（Windows は Alt+Shift+K らしい）を押せばショートカットキーが表示されるので、その内容を読むことでどのようなことができるかある程度把握できるでしょう。
ESS は MELPA からインストールできます。例えば init.el に次のような内容を記述して M-x package-list-packages を実行して ess を選択します。melpa-stable だと今は ess-15.3 がインストールされます。
(when (require 'package nil t) (add-to-list 'package-archives '("melpa-stable" . "http://melpa-stable.milkbox.net/packages/") t) (package-initialize))
ESS の初期設定には思うところがあるので私はこのようにカスタマイズしています。
R の使い方についてさらっと
コンソールの使い方
まず、コンソールの使い方についてさらっとご紹介します。再度 R を起動してください。以降の説明は全てコンソール上で行うことを前提とします。
% R -q # -q オプションを指定するとスタートメッセージを表示しない >
コンソール上では文字を入力した状態で TAB を入力すると補完候補が表示されます。
> read
read.csv read.delim read.fortran read.socket readChar readLines
read.csv2 read.delim2 read.ftable read.table readCitationFile readRDS
read.dcf read.DIF read.fwf readBin readline readRenviron
関数名の手前に ? を付けるか、help 関数に関数名等を指定することでドキュメントを参照することができます。
> # print 関数のドキュメントを参照する > ?print > help("print")
この後 data.table::fread のような書き方や & 演算子が出てきますが、これらのドキュメントを見ることもできます。
> # 特殊なシンボルはバッククォートで括る > ?`::` > ?`&`
また、コンソールにオブジェクト名を入力して評価するとそのオブジェクトの内容が表示されます。
> R.version _ platform x86_64-apple-darwin13.4.0 arch x86_64 os darwin13.4.0 system x86_64, darwin13.4.0 status major 3 minor 2.0 year 2015 month 04 day 16 svn rev 68180 language R version.string R version 3.2.0 (2015-04-16) nickname Full of Ingredients
関数オブジェクトを評価することで関数の実装内容を確認することもできます。
> # print 関数の内容を表示 > print function (x, ...) UseMethod("print") <bytecode: 0x7f9782441a38> <environment: namespace:base>
これで、R の基本的な構文を説明しなくてもドキュメントと組み込み関数の実装を頼りにコードを読み解くことができますね！
デバッグ方法
GDB や Ruby の pry-byebug, pry-stack_explorer のように R にも便利なデバッガが存在します。 使い方の対応表は次のとおりです。
| 内容 | R | GDB | pry-byebug/pry-stack_explorer |
|---|---|---|---|
| コールスタックを表示 | where | bt | show-stack |
| ステップオーバー | n | n | next |
| ステップイン | s | s | step |
| ステップアウト | f | fin | step |
| 次のブレークポイントまで実行 | c | c | continue |
| ヘルプを表示 | help | h | help |
| 処理を中断して終了 | Q | q | quit |
| 関数 f にブレークポイントを設定 | debug(f) | break f | break ClassName#f |
| 関数 f のブレークポイントを削除 | undebug(f) | clear f | break –delete <breakpoint number&rt; |
| 現在の位置にブレークポイントを設定 | browser() | binding.pry |
実際に debug 関数を使って関数にブレークポイントを設定してみます。
まず、エラーに直面したら traceback 関数を実行することでバックトレースを確認できるので、その内容から原因となっている関数を特定します。
> f <- function(x) { y <- as.character(x); g(y) } > g <- function(x) sum(x) > f(1:5) Error in sum(x) : invalid 'type' (character) of argument > traceback() 2: g(y) at #1 1: f(1:5)
バックトレースの内容から g 関数の中でエラーが起きていることがわかるので、debug(g) で g 関数にブレークポイントを設定してデバッグすると良いでしょう。
> debug(g) > f(1:5) debugging in: g(y) debug: sum(x) Browse[2]>
R のデバッガでは、コールスタックを辿って親（呼び出し元）のフレームに移動するようなことはできませんが、parent.frame 関数で親のフレームの環境を取得することができます。
evalq 関数を使えば指定した環境でコードを評価することができるので、親のフレームの状態を確認することも可能です。
Browse[2]> # 親フレーム（関数 f）のオブジェクト一覧 Browse[2]> evalq(ls(), parent.frame()) [1] "x" "y" Browse[2]> # 親フレーム（関数 f）の x の値 Browse[2]> evalq(x, parent.frame()) [1] 1 2 3 4 5
これで不具合に直面してもデバッガを駆使して原因を特定できますね！
edit 関数や trace 関数を使うことで、関数の任意の場所にブレークポイントを設定することもできますが割愛します。
data.frame について
R の特徴的なデータ構造に data.frame があります。行列の各列に異なるクラスのデータを保持できるデータ構造で、データベースのテーブルのようなデータ構造をイメージしていただくと良いと思います。
ログの分析では data.frame を使うことがほとんどでしょう。
> # R のデータの例でよく使われる iris データ > head(iris) Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa > # 各列のクラス > sapply(iris, class) Sepal.Length Sepal.Width Petal.Length Petal.Width Species "numeric" "numeric" "numeric" "numeric" "factor" > # 1 行目のデータを取得 > iris[1, ] Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa > # Species の列を取得 > head(iris$Species) [1] setosa setosa setosa setosa setosa setosa Levels: setosa versicolor virginica
A/B テストの結果検証
さくっと R の導入から使い方まで説明したところで、擬似的なログを使って簡単な分析をしてみます。 次のような背景から分析することになったとします。
あるサービスのユーザ数を増やすために、チーム内でランディングページ (以降 LP) の改修を検討しています。既存の LP と新しい LP で 4/1 〜 4/7 の 1 週間 A/B テストを行ったので、新しい LP の効果を検証したいです。
準備
まず、今回の分析で使うパッケージをインストールします。各パッケージの概要については使う際に説明します。
> pkgs <- c("data.table", "dplyr", "tidyr", "ggplot2") > install.packages(pkgs) > # 各パッケージのバージョン > sapply(pkgs, packageVersion, simplify = FALSE) $data.table [1] ‘1.9.4’ $dplyr [1] ‘0.4.1’ $tidyr [1] ‘0.2.0’ $ggplot2 [1] ‘1.0.1’
次に、今回使用する擬似的なログを作成します。ここに擬似データを作成する create_sample_data 関数を定義したスクリプトを配置したので、読み込んだ後に関数を実行してください。
スクリプトの読み込みには source 関数を使います。source 関数は主にローカルのファイルのパスを指定しますが、このように URL を指定することもできます。
> # R 3.2.0 でないと https には対応していないので注意 > source("https://gist.githubusercontent.com/abicky/3a4789c3fd163a71606c/raw/5f1aeb86b8f0eb50caf386aa3fce9bc5354df9b5/create_sample_data.R") > # 擬似データ（40 MB 程度）の作成 > file_paths <- create_sample_data() Creating logs of 'a' pattern .......... done Creating logs of 'b' pattern .......... done Writing files done
create_sample_data 関数を実行することで一時ディレクトリに event_logs.tsv と access_logs.tsv が作成されます。返り値は list で、event_log_file と access_log_file にそれぞれの絶対パスが格納されています。
それでは、event_logs.tsv と access_logs.tsv を読み込みましょう。
小規模なデータの読み込みには read.table 関数などを使うのが一般的ですが、read.table 系の関数は非常に遅いので data.table パッケージの fread 関数を使います*3。data.table パッケージは、data.frame を高速化したような data.table クラスを提供するパッケージですが、今回はデータの読み込みのためだけに使用します。なお、fread 関数で読み込んだデータは data.table クラスのオブジェクトになります。
> # file_paths$event_log_file には event_logs.tsv の絶対パスが格納されている > event_logs <- data.table::fread(file_paths$event_log_file) > event_logs time user_id event pattern 1: 1428373621 201681931 imp a 2: 1428299552 898389685 imp a 3: 1427862703 944675268 imp a 4: 1428109708 660797792 imp a 5: 1428236105 629114044 imp a --- 259289: 1427877582 671321141 cv b 259290: 1428203926 733054604 cv b 259291: 1427948288 590425796 cv b 259292: 1428140625 29272541 cv b 259293: 1428052793 466384837 cv b > # file_paths$access_log_file には access_logs.tsv の絶対パスが格納されている > access_logs <- data.table::fread(file_paths$access_log_file) > access_logs time user_id page 1: 1427862710 944675268 page1 2: 1427862739 944675268 page3 3: 1427862750 944675268 page1 4: 1428236117 629114044 page2 5: 1428236120 629114044 page1 --- 1251843: 1430046543 466384837 page1 1251844: 1430220035 466384837 page2 1251845: 1430220043 466384837 page1 1251846: 1430220055 466384837 page3 1251847: 1430220071 466384837 page2
event_logs の time はそのイベントが発生した unix time、event は “imp” か “cv” でそれぞれ LP の impression と conversion（ここではユーザ登録）を意味しています。pattern は a が新しい LP で b が既存の LP です。
access_logs の time はユーザがアクセスした unix time、page はアクセスしたページを意味しています。今回 page は使いませんがアクセスログっぽさを出すために苦労して付与してみました。
コンバージョン率 (CVR) を出す
CVR の算出など一連のデータ操作には dplyr パッケージを使います。dplyr パッケージは組み込みのデータ操作用の関数に比べて高速な上、SQL に近い感覚で処理を記述できるのでエンジニアにとっては馴染みやすいでしょう。*4
まず、読み込んだデータを tbl_df 関数で dplyr 用の data.frame 表現に変換します。data.table オブジェクトのままでも処理できますが、次の理由から変換することにします。
data.tableオブジェクトはPOSIXltオブジェクトを扱えないなどの制約がある- 私の手元の環境では
tbl_dfで変換した方が若干高速だった data.tableオブジェクトとdata.frameオブジェクト（read.table等で読み込んだ場合のオブジェクト）で dplyr の処理結果が異なることがあるのでどちらかに統一しておきたい
> # パッケージのロード（C++ の using namespace dplyr、Python の from dplyr import * に相当） > library(dplyr) Attaching package: 'dplyr' The following object is masked from 'package:stats': filter The following objects are masked from 'package:base': intersect, setdiff, setequal, union > event_log_df <- tbl_df(event_logs)
さて、CVR を出すには imp と cv の UU を出す必要があります。SQL だと次のようなクエリですね。
SELECT pattern, event, COUNT(DISTINCT(user_id)) uu FROM event_logs GROUP BY pattern, event ;
上記の SQL では次のような手順を踏んでいます。
- pattern, event でグループ化する
- それぞれのグループに関してユニークな user_id をカウントする
これに対応する dplyr の処理は次のようになります。
> event_counts <- event_log_df %>% + # 1. pattern, event でグループ化する + group_by(pattern, event) %>% + # 2. それぞれのグループに関してユニークな user_id をカウントする + summarize(uu = n_distinct(user_id)) > event_counts Source: local data frame [4 x 3] Groups: pattern pattern event uu 1 a cv 30097 2 a imp 100399 3 b cv 28901 4 b imp 99896
%>% は左辺のデータを右辺の関数の第一引数に流す演算子で、dplyr パッケージが提供している演算子です*5。dplyr による処理ではこの演算子を使った記法が好まれます。左側の処理の出力が右側の処理の入力になるので pipe (|) をイメージするとわかりやすいでしょう。
ちなみに、%>% を使わないと次のような書き方になります。
event_counts <- dplyr::summarize( dplyr::group_by(event_log_df, pattern, event), uu = n_distinct(user_id) )
さて、現在次のような縦長のデータになっていて、cv / imp を出すのが少し面倒です。
| pattern | event | uu |
|---|---|---|
| a | cv | 30097 |
| a | imp | 100399 |
| b | cv | 28901 |
| b | imp | 99896 |
これを次のような横長のデータに変換できれば cv の列を imp の列で割るだけで CVR が算出できます。
| pattern | cv | imp |
|---|---|---|
| a | 30097 | 100399 |
| b | 28901 | 99896 |
横長のデータに変換するには tidyr パッケージの spread 関数を使います。tidyr パッケージは縦長のデータを横長にしたり、列を分割したりデータの整形に便利な関数を提供しているパッケージです。
spread 関数の第 1 引数には列に並べたいフィールドを指定します。今回の場合は imp, cv が列に並んでほしいので event を指定することになります。第 2 引数には、第 1 引数に指定したフィールドによって作成される列の各セルに表示する値となるフィールドを指定します。今回の場合は event ごとの uu の値を表示したいので uu を指定します。
> event_counts %>% + tidyr::spread(event, uu) Source: local data frame [2 x 3] pattern cv imp 1 a 30097 100399 2 b 28901 99896
event ごとに列が割り当てられましたね。 あとは cv を imp で割るだけです。
> event_counts_with_cvr <- event_counts %>% + tidyr::spread(event, uu) %>% + # cvr という列を追加 + dplyr::mutate(cvr = cv / imp) > event_counts_with_cvr Source: local data frame [2 x 4] pattern cv imp cvr 1 a 30097 100399 0.2997739 2 b 28901 99896 0.2893109
結果を見る限り、pattern a（新しい LP）の CVR が 30.0 % で pattern b の 28.9 % より高いですね！
CVR の差の検定をする
先ほどの結果より、新しい LP によって CVR が 28.9 % から 30.0 % に向上したことがわかりました。これは有意な差なのでしょうか？
比率の差が有意かどうかを確認するには prop.test 関数を使用します。prop.test 関数ではカイ二乗検定による母比率の差の検定を行います。第 1 引数には成功回数、第 2 引数には試行回数を指定します。今回のケースだと第 1 引数に cv の数、第 2 引数に imp の数を指定することになります。alternative には片側検定をする場合に “greater” か “less” を指定します。指定しない場合は両側検定になります。
> # event_counts_with_cvr の cv 列の内容を取得 > cv <- event_counts_with_cvr$cv > # event_counts_with_cvr の imp 列の内容を取得 > imp <- event_counts_with_cvr$imp > prop.test(cv, imp, alternative = "greater") 2-sample test for equality of proportions with continuity correction data: cv out of imp X-squared = 26.331, df = 1, p-value = 1.438e-07 alternative hypothesis: greater 95 percent confidence interval: 0.007102618 1.000000000 sample estimates: prop 1 prop 2 0.2997739 0.2893109
出力内容についてざっくり説明します。
“2-sample test for equality of proportions with continuity correction” は 2 郡の比率が等しいかどうかを検定したことを意味しています。その際、連続性の補正 (continuity correction) を行った旨が表示されていますが、連続性の補正の有無で経営判断が変わるほど結果に影響を及ぼすことはまずないと思われるので気にしなくて大丈夫です*6。
“X-squared = 26.331, df = 1, p-value = 1.438e-07” はカイ二乗検定のカイ二乗値などを表しています。この中で重要なのが p-value です。この値は、pattern a の CVR と pattern b の CVR が等しい場合に今回のような結果が得られる確率を表しています。CVR が同じである確率と解釈しても A/B テストの検証においては差し支えないでしょう*7。CVR が同じである確率と捉えた場合、同じである確率は 1.44 × 10-7 と極めて低いと言えます。
あとの出力内容は片側検定であることや、95 % 信頼区間の情報やそれぞれの CVR を表しています。
以上の結果より、今回の結果は有意な差と言えそうです。
ちなみに、p-value の目安として 0.05 を切ったら有意な差（5 % 有意水準）とすることが多いですが、経営判断する上では CVR 以外の様々な要因（e.g. メンテナンスコスト）もあるので、あくまで参考程度にするのが良いと思います。
ユーザの定着率（定着数）を出す
新しい LP によって CVR が上がったことを確認できましたが、ユーザの定着率が下がってしまっては意味がありません。 ということで、次は定着率を出してみます。 定着率の定義として、ここではユーザ登録したユーザが n 日後にアクセスした割合を採用します*8。
少し複雑な処理になるので、SQL の書き方にも様々な選択肢がありますが、例えば次のようなクエリで各パターンの経過日数ごとの UU と定着率を算出できます（PostgreSQL）。
SELECT pattern, elapsed_days, uu, uu::float / first_value(uu) OVER (PARTITION BY pattern ORDER BY elapsed_days) retention_ratio FROM ( SELECT pattern, elapsed_days, COUNT(DISTINCT a.user_id) uu FROM ( SELECT user_id, to_timestamp(time)::date - to_timestamp(min(time) OVER (PARTITION BY user_id))::date elapsed_days FROM access_logs ) a INNER JOIN event_logs e ON a.user_id = e.user_id WHERE event = 'cv' AND elapsed_days < 14 GROUP BY pattern, elapsed_days ) t ORDER BY pattern, elapsed_days ;
上記の SQL では次のような手順を踏んでいます。
- access_logs を user_id でグループ化する
- グループごとに最小の unix time を算出する
- 2 で算出した時間と time を日に変換してから差分を算出して elapsed_days とする
- event_logs と user_id で inner join する
- event が ‘cv’ で elapsed_days が 14 より小さいレコードに限定する
- pattern, elapsed_days でグループ化する
- グループごとにユニークな user_id をカウントする
- pattern でグループ化する
- グループごとに elapsed_days が最小のレコードの uu の値で uu を割る
- pattern, elapsed_days でソートする
これに対応する dplyr の処理は次のようになります。
> # dplyr 用の data.frame に変換 > access_log_df <- dplyr::tbl_df(access_logs) > # unix time を date に変換する関数 > time_to_date <- function(time) { + as.Date(as.POSIXct(time, origin = "1970-01-01", tz = "Asia/Tokyo")) + } > uu_transitions <- access_log_df %>% + # 1. access_logs を user_id でグループ化する + group_by(user_id) %>% + # 2. グループごとに最小の unix time を算出する + mutate(min_time = min(time)) %>% + # グループ化を解除（集約処理が終わったら ungroup するのが無難） + ungroup() %>% + # 3. 2 で算出した時間と time を日に変換してから差分を算出して elapsed_days とする + mutate(elapsed_days = time_to_date(time) - time_to_date(min_time)) %>% + # 4. event_logs と user_id で inner join する + inner_join(event_log_df, by = "user_id") %>% + # 5. event が 'cv' で elapsed_days が 14 より小さいレコードに限定する + filter(event == "cv" & elapsed_days < 14) %>% + # 6. pattern, elapsed_days でグループ化する + group_by(pattern, elapsed_days) %>% + # 7. グループごとにユニークな user_id をカウントする + summarize(uu = n_distinct(user_id)) %>% + # 8. pattern でグループ化する + group_by(pattern) %>% + # 9. グループごとに elapsed_days が最小のレコードの uu の値で uu を割る + mutate(retention_ratio = uu / first(uu, order_by = "elapsed_days")) %>% + # グループ化を解除（集約処理が終わったら ungroup するのが無難） + ungroup() %>% + # 10. pattern, elapsed_days でソートする + arrange(pattern, elapsed_days)
CVR を出した時の処理に比べると格段に複雑になりました。わかりにくそうな関数について少し補足説明します。
| 関数 | 概要 |
|---|---|
filter |
条件にマッチするデータのみ抽出する（SQL の WHERE 句に相当） |
arrange |
指定したフィールドでソート（SQL の ORDER BY 句に相当） |
mutate |
指定したフィールドを追加（SQL の SELECT *, new_field に相当） |
group_by 関数を適用した後に mutate 関数を適用すると SQL の PARTION BY 相当の処理が行われます。
filter 関数の引数に & 演算子を使っていますが、これはベクトルの論理積を意味しています。ベクトルの要素ごとに AND 条件で評価し、その結果をベクトルで返す演算子です。なお、ベクトルの論理和には | 演算子を使います。
> # Ruby の (1..5).to_a に相当 > x <- 1:5 > # 2 より大きい要素は TRUE になる > x > 2 [1] FALSE FALSE TRUE TRUE TRUE > # 4 より小さい要素は TRUE になる > x < 4 [1] TRUE TRUE TRUE FALSE FALSE > # ベクトルの論理積。2 より大きく 4 より小さい要素は TRUE になる > x > 2 & x < 4 [1] FALSE FALSE TRUE FALSE FALSE
filter 関数は指定したベクトルの要素が TRUE の場合に対応する行のデータを抽出します。上記の例では event の値が “cv” で elapsed_days が 14 より小さいデータを抽出しています。
ちなみに、上記の処理をもう少し最適化すると次のようになります。不要なデータを早い段階で削除しています。
uu_transitions <- access_log_df %>% group_by(user_id) %>% mutate(min_time = min(time)) %>% ungroup() %>% mutate(elapsed_days = time_to_date(time) - time_to_date(min_time)) %>% filter(elapsed_days < 14) %>% select(user_id, elapsed_days) %>% distinct(user_id, elapsed_days) %>% inner_join( event_log_df %>% filter(event == "cv") %>% select(user_id, pattern), by = "user_id" ) %>% group_by(pattern, elapsed_days) %>% summarize(uu = n_distinct(user_id)) %>% group_by(pattern) %>% mutate(retention_ratio = uu / first(uu, order_by = "elapsed_days")) %>% ungroup() %>% arrange(pattern, elapsed_days)
さて、本題に戻って結果を見てみましょう。
> # print 関数の n 引数に表示する行数を指定することで全データを表示 > print(uu_transitions, n = nrow(uu_transitions)) Source: local data frame [28 x 4] pattern elapsed_days uu retention_ratio 1 a 0 days 30097 1.0000000 2 a 1 days 11270 0.3744559 3 a 2 days 9747 0.3238529 4 a 3 days 9741 0.3236535 5 a 4 days 8830 0.2933847 6 a 5 days 8175 0.2716218 7 a 6 days 7716 0.2563711 8 a 7 days 7223 0.2399907 9 a 8 days 6946 0.2307871 10 a 9 days 6602 0.2193574 11 a 10 days 6437 0.2138751 12 a 11 days 6053 0.2011164 13 a 12 days 5752 0.1911154 14 a 13 days 5555 0.1845699 15 b 0 days 28901 1.0000000 16 b 1 days 12719 0.4400886 17 b 2 days 11379 0.3937234 18 b 3 days 11162 0.3862150 19 b 4 days 10054 0.3478772 20 b 5 days 9377 0.3244524 21 b 6 days 8766 0.3033113 22 b 7 days 8283 0.2865991 23 b 8 days 7944 0.2748694 24 b 9 days 7552 0.2613058 25 b 10 days 7214 0.2496107 26 b 11 days 6929 0.2397495 27 b 12 days 6541 0.2263243 28 b 13 days 6261 0.2166361
初日の uu は pattern a の方が大きいですが、全体的に pattern a の定着率は pattern b に比べて低く、その結果 14 日目 (elapsed_days が 13) の uu は pattern a の方が小さくなっています。
値の推移については数字を眺めるよりも可視化した方がわかりやすいでしょう。可視化には ggplot2 パッケージを使うことにします。
ggplot2 は折れ線グラフ、面グラフ、ヒストグラム等を統一的なインタフェースで描画できるパッケージで、複数のグラフの重ね合わせや軸の分割など複雑なグラフも手軽に描画できることが特徴です。
plot 関数等を使ったグラフの描画に慣れている方にとって ggplot2 は使いにくいかもしれませんが、R を初めて使うエンジニアの方には ggplot2 を使った描画方法を習得することをお勧めします。
> library(ggplot2) Loading required package: methods > # グラフにデータ (uu_transitions) をセットし、横軸を elapsed_days、縦軸を uu にする > g <- ggplot(uu_transitions, aes(x = as.integer(elapsed_days), y = uu)) > # 折れ線グラフ (geom_line) を描画する。その際 pattern によって色を変える > g + geom_line(aes(color = pattern))
上記の例では pattern によって色を変えて描画していますが、例えばデモグラ情報（年代、性別等の情報）を持ったデータであれば、年代ごとに異なる色を使い、性別ごとに異なる線種を使うようなことも可能です。 ggplot2 の使い方についてはドキュメメントに豊富な例が載っているのでそちらを参照してください。次の資料も ggplot2 を使う上で非常に参考になります。
さて、グラフを見ると pattern a は 2 日目以降 pattern b よりも UU が低いことがひと目でわかります。これは LP の内容からユーザが期待する内容とサービスの内容に差異があったと考えられます。
以上を踏まえると、現状のまま新しい LP を導入するのは避けた方が良いという判断になるでしょう。
番外編
R でコホート分析
最近 Google Analytics にコホート分析の機能が追加されましたね。R であのような表を作成してみましょう。縦軸にユーザ登録日、横軸に経過日数、値に定着率を取ることにします。
データは先程も利用した access_log_df を使います。今回はユーザ登録日を付与していることと、レコード数を絞るためにログの期間を 4/7 までに限定していること以外は A/B テストのユーザ定着率を出した時とほとんど変わりません。
> cohort <- access_log_df %>% + group_by(user_id) %>% + mutate(min_time = min(time)) %>% + ungroup() %>% + mutate( + registraion_date = time_to_date(min_time), + elapsed_days = time_to_date(time) - time_to_date(min_time) + ) %>% + filter(time < as.POSIXct("2015-04-08", tz = "Asia/Tokyo")) %>% + group_by(registraion_date, elapsed_days) %>% + summarize(uu = n_distinct(user_id)) %>% + group_by(registraion_date) %>% + mutate(retention_ratio = round(uu / first(uu, order_by = elapsed_days), 3)) > cohort Source: local data frame [36 x 4] Groups: registraion_date registraion_date elapsed_days uu retention_ratio 1 2015-03-31 0 days 532 1.000 2 2015-03-31 1 days 187 0.352 3 2015-03-31 2 days 195 0.367 4 2015-03-31 3 days 188 0.353 5 2015-03-31 4 days 191 0.359 6 2015-03-31 5 days 176 0.331 7 2015-03-31 6 days 145 0.273 8 2015-03-31 7 days 31 0.058 9 2015-04-01 0 days 8351 1.000 10 2015-04-01 1 days 3326 0.398 .. ... ... ... ...
これをおなじみ tidyr::spread 関数で横長のデータに変換します。
> cohort %>% + select(registraion_date, elapsed_days, retention_ratio) %>% + tidyr::spread(elapsed_days, retention_ratio) %>% + arrange(registraion_date) Source: local data frame [8 x 9] registraion_date 0 1 2 3 4 5 6 7 1 2015-03-31 1 0.352 0.367 0.353 0.359 0.331 0.273 0.058 2 2015-04-01 1 0.398 0.358 0.355 0.322 0.297 0.250 NA 3 2015-04-02 1 0.411 0.356 0.350 0.312 0.264 NA NA 4 2015-04-03 1 0.396 0.348 0.361 0.287 NA NA NA 5 2015-04-04 1 0.415 0.359 0.322 NA NA NA NA 6 2015-04-05 1 0.417 0.336 NA NA NA NA NA 7 2015-04-06 1 0.368 NA NA NA NA NA NA 8 2015-04-07 1 NA NA NA NA NA NA NA
簡単ですね！
最後に
以上、単純な例ではありましたが、A/B テストの結果を検証をする際にどのような検証を行うか、R を使う場合にどうすれば良いかがなんとなく理解いただけたのではないかと思います。 もちろん、より単純な検証を求められるケース、より詳細な検証を求められるケース、全く別の検証を求められるケースなど状況によって要求は様々なので、状況に応じて検証方法を検討してください。
今回は単純な集計処理しか行いませんでしたが、これから R で本格的に分析をしようと思っている方は以下の資料にも目を通すことをお勧めします。 Tokyo.R のような R コミュニティにも一度顔を出してみると良いと思います。
- Rデータ自由自在
- R の基本的なデータ処理についてコンパクトにまとめられているので個人的にオススメです
- 同志社大学の金先生のページ
- 統計学や機械学習の手法を R で行う方法について非常に参考になります（かなりの量なのでこういうページがあるということだけでも覚えておくと良いでしょう）
- dplyr の vignettes
- dplyr を使いこなしたい方は一読すべきでしょう
- ハンバーガー統計学にようこそ！
- R とは関係ないですが検定についてわかりやすく解説されています
- アイスクリーム統計学にようこそ！
- 「ハンバーガー統計学にようこそ！」と合わせて読むと良いでしょう
あとは拙作のマニアックなスライドですがいくつか参考になるかもしれません。
- Rデバッグあれこれ
- 今回触れたデバッグ方法よりも詳細な情報を載せています（ただ、情報が古い・・・）
- Rデータフレーム自由自在 〜基本操作からピボットテーブルまで〜
- dplyr を使う場合は見る価値のないスライドです。dplyr がないとどのような処理になるか知りたい方に。
- Rのスコープとフレームと環境と 〜parent.env と parent.frame の違い〜
- デバッグ方法のところで軽く触れましたがフレームなどの話について触れています
- Rのデータ構造とメモリ管理
ではでは快適な R ライフを！！
※このエントリーは knitr パッケージを使って作成しました
*1:コーディングスタイルに迷ったら R 界で影響力の大きい Hadley Wickham 氏のコーディングスタイルに従うのが良いと思います
*2:執筆当時は筑波大学のリポジトリを指定していましたが、2015 年の 6 月に閉鎖したので変更しました（2017-04-22 追記）
*3:手元の環境だと、アクセスログを読み込むのに read.table は data.table::fread の 100 倍近く時間がかかります
*4:まだ枯れていないのでバグを踏むことはありますが・・・
*5:厳密には magrittr パッケージが提供しているものを dplyr パッケージが取り込んでいます
*6:気になるようであれば correct 引数に FALSE を指定することで補正を行わなくすることも可能です
*7:A であれば B が得られる確率が p-value であって、B が得られたから A である確率ではないことを強調しておきます
*8:個人的に、普段の分析では n 週間後にアクセスした割合を採用しています。ソーシャルゲームのように毎日アクセスすることを前提としたサービスでは n 日後にアクセスした割合で良いと思います
