Skip to content

Instantly share code, notes, and snippets.

@abicky

abicky/techlife-20150508.Rmd Secret

Last active November 9, 2015 02:00
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save abicky/d3deb617bddda461e2eb to your computer and use it in GitHub Desktop.
Save abicky/d3deb617bddda461e2eb to your computer and use it in GitHub Desktop.
Download ZIP
Raw
techlife-20150508.Rmd
# A/B テストで施策の効果を検証！エンジニアのための R 入門
```{r setup, echo = FALSE}
library(knitr)
opts_chunk$set(comment = NA, tidy = FALSE, prompt = TRUE, collapse = TRUE)
```
こんにちは、買物情報事業部でサーバサイドの開発を担当している荒引 ([@a_bicky](https://twitter.com/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 のコーディングスタイルのカオスっぷりに辟易する方もいるでしょうが、そこは耐えるしかないです。((コーディングスタイルに迷ったら R 界で影響力の大きい [Hadley Wickham 氏のコーディングスタイル](http://adv-r.had.co.nz/Style.html)に従うのが良いと思います))
## アジェンダ
* R の環境構築
- 最低限の設定
- IDE の導入
* R の使い方についてさらっと
- コンソールの使い方
- デバッグ方法
- data.frame について
* A/B テストの結果検証
- コンバージョン率 (CVR) を出す
- CVR の差の検定をする
- ユーザの定着率（定着数）を出す
* 番外編
- R でコホート分析
* 最後に
## R の環境構築
Mac の場合は Homebrew 経由でインストールするのが手軽です。Linux の場合もパッケージマネージャ経由で手軽にインストールできます。
```zsh
% brew tap homebrew/science
% brew install r
```
これでターミナルから 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 を入力します。
```r
> q()
Save workspace image? [y/n/c]: n
```
"Save workspace image?" で y を入力すると現在のセッションの内容がワーキングディレクトリに .RData という名前で保存され、次回起動時に定義したオブジェクトなどが復元されます。
### 最低限の設定
個人的に次の 2 つの設定だけは欠かせないと思っています。
* 言語の英語化
- エラーメッセージでググった時に参考になる情報が多いので、英語アレルギーでない限りは英語化した方が良いでしょう
* リポジトリの設定
- パッケージをインストールする際にいちいち指定する必要がなくなります
次のコマンドを実行することで上記の設定が可能です。
```zsh
% # 言語の英語化
% echo 'LANGUAGE=en' >> ~/.Renviron
% # パッケージをインストールする際のリポジトリを設定
% echo 'options(repos = "http://cran.md.tsukuba.ac.jp")' >> ~/.Rprofile
```
.Renviron は R を実行する上での環境変数を定義するファイルです。
.Rprofile は起動時の処理を記述するファイルです。`options` 関数の `repos` 引数に筑波大学のリポジトリを指定しています。
ただ、筑波大学は休日に停電を行うことがあるので、その際は別のリポジトリを指定しないとパッケージをインストールできません。
### IDE の導入
前述したように、R を使う際には基本的にコンソール上で作業することになりますが、それなりの処理を行う場合は IDE を利用すると便利です。
IDE を利用することで、スクリプト上のコードをコンソールに送って挙動を確認するようなことも可能になります。
メジャーな IDE は次の 2 つです。
* [RStudio](http://www.rstudio.com/products/rstudio/)
* [ESS](http://ess.r-project.org/)（Emacs のメジャーモード）
RStudio を使うには、[ダウンロードリンク](http://www.rstudio.com/products/rstudio/download/)から自分の環境にあったファイルをダウンロードしてインストールするだけです。
Option+Shift+K（Windows は Alt+Shift+K らしい）を押せばショートカットキーが表示されるので、その内容を読むことでどのようなことができるかある程度把握できるでしょう。
ESS は MELPA からインストールできます。例えば init.el に次のような内容を記述して `M-x package-list-packages` を実行して ess を選択します。melpa-stable だと今は ess-15.3 がインストールされます。
```lisp
(when (require 'package nil t)
(add-to-list 'package-archives '("melpa-stable" . "http://melpa-stable.milkbox.net/packages/") t)
(package-initialize))
```
ESS の初期設定には思うところがあるので私は[このように](https://gist.github.com/abicky/5da654e3a7858ec60aed)カスタマイズしています。
## R の使い方についてさらっと
### コンソールの使い方
まず、コンソールの使い方についてさらっとご紹介します。再度 R を起動してください。以降の説明は全てコンソール上で行うことを前提とします。
```zsh
% R -q # -q オプションを指定するとスタートメッセージを表示しない
>
```
コンソール上では文字を入力した状態で TAB を入力すると補完候補が表示されます。
```r
> 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` 関数に関数名等を指定することでドキュメントを参照することができます。
```{r, eval = FALSE}
# print 関数のドキュメントを参照する
?print
help("print")
```
この後 `data.table::fread` のような書き方や `&` 演算子が出てきますが、これらのドキュメントを見ることもできます。
```{r, eval = FALSE}
# 特殊なシンボルはバッククォートで括る
?`::`
?`&`
```
また、コンソールにオブジェクト名を入力して評価するとそのオブジェクトの内容が表示されます。
```{r}
R.version
```
関数オブジェクトを評価することで関数の実装内容を確認することもできます。
```{r}
# print 関数の内容を表示
print
```
これで、R の基本的な構文を説明しなくてもドキュメントと組み込み関数の実装を頼りにコードを読み解くことができますね！
### デバッグ方法
GDB や Ruby の [pry-byebug](https://github.com/deivid-rodriguez/pry-byebug), [pry-stack_explorer](https://github.com/pry/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` 関数を実行することでバックトレースを確認できるので、その内容から原因となっている関数を特定します。
```r
> 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` 関数にブレークポイントを設定してデバッグすると良いでしょう。
```r
> debug(g)
> f(1:5)
debugging in: g(y)
debug: sum(x)
Browse[2]>
```
R のデバッガでは、コールスタックを辿って親（呼び出し元）のフレームに移動するようなことはできませんが、`parent.frame` 関数で親のフレームの環境を取得することができます。
`evalq` 関数を使えば指定した環境でコードを評価することができるので、親のフレームの状態を確認することも可能です。
```r
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}
# R のデータの例でよく使われる iris データ
head(iris)
# 各列のクラス
sapply(iris, class)
# 1 行目のデータを取得
iris[1, ]
# Species の列を取得
head(iris$Species)
```
## A/B テストの結果検証
さくっと R の導入から使い方まで説明したところで、擬似的なログを使って簡単な分析をしてみます。
次のような背景から分析することになったとします。
> あるサービスのユーザ数を増やすために、チーム内でランディングページ (以降 LP) の改修を検討しています。既存の LP と新しい LP で 4/1 〜 4/7 の 1 週間 A/B テストを行ったので、新しい LP の効果を検証したいです。
### 準備
まず、今回の分析で使うパッケージをインストールします。各パッケージの概要については使う際に説明します。
```r
> 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’
```
次に、今回使用する擬似的なログを作成します。[ここ](https://gist.github.com/abicky/3a4789c3fd163a71606c/)に擬似データを作成する `create_sample_data` 関数を定義したスクリプトを配置したので、読み込んだ後に関数を実行してください。
スクリプトの読み込みには `source` 関数を使います。`source` 関数は主にローカルのファイルのパスを指定しますが、このように URL を指定することもできます。
```{r create_sample_data, cache = TRUE}
# 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()
```
`create_sample_data` 関数を実行することで一時ディレクトリに event_logs.tsv と access_logs.tsv が作成されます。返り値は list で、`event_log_file` と `access_log_file` にそれぞれの絶対パスが格納されています。
それでは、eventlogs.tsv と access_logs.tsv を読み込みましょう。
小規模なデータの読み込みには `read.table` 関数などを使うのが一般的ですが、`read.table` 系の関数は非常に遅いので [data.table](https://github.com/Rdatatable/data.table) パッケージの `fread` 関数を使います((手元の環境だと、アクセスログを読み込むのに read.table は data.table::fread の 100 倍近く時間がかかります))。data.table パッケージは、`data.frame` を高速化したような `data.table` クラスを提供するパッケージですが、今回はデータの読み込みのためだけに使用します。なお、`fread` 関数で読み込んだデータは `data.table` クラスのオブジェクトになります。
```{r, cache = TRUE, dependson = "create_sample_data"}
# file_paths$event_log_file には event_logs.tsv の絶対パスが格納されている
event_logs <- data.table::fread(file_paths$event_log_file)
event_logs
# file_paths$access_log_file には access_logs.tsv の絶対パスが格納されている
access_logs <- data.table::fread(file_paths$access_log_file)
access_logs
```
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](https://github.com/hadley/dplyr) パッケージを使います。dplyr パッケージは組み込みのデータ操作用の関数に比べて高速な上、SQL に近い感覚で処理を記述できるのでエンジニアにとっては馴染みやすいでしょう。((まだ枯れていないのでバグを踏むことはありますが・・・))
まず、読み込んだデータを `tbl_df` 関数で dplyr 用の data.frame 表現に変換します。`data.table` オブジェクトのままでも処理できますが、次の理由から変換することにします。
1. `data.table` オブジェクトは `POSIXlt` オブジェクトを扱えないなどの制約がある
2. 私の手元の環境では `tbl_df` で変換した方が若干高速だった
3. `data.table` オブジェクトと `data.frame`オブジェクト（`read.table` 等で読み込んだ場合のオブジェクト）で dplyr の処理結果が異なることがあるのでどちらかに統一しておきたい
```{r}
# パッケージのロード（C++ の using namespace dplyr、Python の from dplyr import * に相当）
library(dplyr)
event_log_df <- tbl_df(event_logs)
```
さて、CVR を出すには imp と cv の UU を出す必要があります。SQL だと次のようなクエリですね。
<!--
CREATE TABLE event_logs (time int, user_id int, event character(3), pattern character(1));
COPY event_logs FROM '/Users/arabiki/work/r/event_logs.tsv' (DELIMITER ' ', FORMAT CSV, HEADER true);
CREATE TABLE access_logs (time int, user_id int, page character(5));
COPY access_logs FROM '/Users/arabiki/work/r/access_logs.tsv' (DELIMITER ' ', FORMAT CSV, HEADER true);
-->
```sql
SELECT
pattern,
event,
COUNT(DISTINCT(user_id)) uu
FROM
event_logs
GROUP BY
pattern,
event
;
```
上記の SQL では次のような手順を踏んでいます。
1. pattern, event でグループ化する
2. それぞれのグループに関してユニークな user_id をカウントする
これに対応する dplyr の処理は次のようになります。
```{r}
event_counts <- event_log_df %>%
# 1. pattern, event でグループ化する
group_by(pattern, event) %>%
# 2. それぞれのグループに関してユニークな user_id をカウントする
summarize(uu = n_distinct(user_id))
event_counts
```
`%>%` は左辺のデータを右辺の関数の第一引数に流す演算子で、dplyr パッケージが提供している演算子です((厳密には magrittr パッケージが提供しているものを dplyr パッケージが取り込んでいます))。dplyr による処理ではこの演算子を使った記法が好まれます。左側の処理の出力が右側の処理の入力になるので pipe (`|`) をイメージするとわかりやすいでしょう。
ちなみに、`%>%` を使わないと次のような書き方になります。
```r
event_counts <- dplyr::summarize(
dplyr::group_by(event_log_df, pattern, event),
uu = n_distinct(user_id)
)
```
さて、現在次のような縦長のデータになっていて、cv / imp を出すのが少し面倒です。
```{r, echo = FALSE, results = "asis"}
kable(event_counts)
```
これを次のような横長のデータに変換できれば cv の列を imp の列で割るだけで CVR が算出できます。
```{r, echo = FALSE, results = "asis"}
kable(event_counts %>% tidyr::spread(event, uu))
```
横長のデータに変換するには [tidyr](https://github.com/hadley/tidyr) パッケージの `spread` 関数を使います。tidyr パッケージは縦長のデータを横長にしたり、列を分割したりデータの整形に便利な関数を提供しているパッケージです。
`spread` 関数の第 1 引数には列に並べたいフィールドを指定します。今回の場合は imp, cv が列に並んでほしいので event を指定することになります。第 2 引数には、第 1 引数に指定したフィールドによって作成される列の各セルに表示する値となるフィールドを指定します。今回の場合は event ごとの uu の値を表示したいので uu を指定します。
```{r}
event_counts %>%
tidyr::spread(event, uu)
```
event ごとに列が割り当てられましたね。
あとは cv を imp で割るだけです。
```{r}
event_counts_with_cvr <- event_counts %>%
tidyr::spread(event, uu) %>%
# cvr という列を追加
dplyr::mutate(cvr = cv / imp)
event_counts_with_cvr
```
```{r, echo = FALSE}
cvr_a <- sprintf("%.1f", (event_counts_with_cvr %>% dplyr::filter(pattern == "a"))$cvr * 100)
cvr_b <- sprintf("%.1f", (event_counts_with_cvr %>% dplyr::filter(pattern == "b"))$cvr * 100)
```
結果を見る限り、pattern a（新しい LP）の CVR が `r cvr_a` % で pattern b の `r cvr_b` % より高いですね！
### CVR の差の検定をする
先ほどの結果より、新しい LP によって CVR が `r cvr_b` % から `r cvr_a` % に向上したことがわかりました。これは有意な差なのでしょうか？
比率の差が有意かどうかを確認するには `prop.test` 関数を使用します。`prop.test` 関数ではカイ二乗検定による母比率の差の検定を行います。第 1 引数には成功回数、第 2 引数には試行回数を指定します。今回のケースだと第 1 引数に cv の数、第 2 引数に imp の数を指定することになります。`alternative` には片側検定をする場合に "greater" か "less" を指定します。指定しない場合は両側検定になります。
```{r}
# 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")
```
```{r, echo = FALSE}
test_result <- prop.test(cv, imp, alternative = "greater")
```
出力内容についてざっくり説明します。
"2-sample test for equality of proportions with continuity correction" は 2 郡の比率が等しいかどうかを検定したことを意味しています。その際、連続性の補正 (continuity correction) を行った旨が表示されていますが、連続性の補正の有無で経営判断が変わるほど結果に影響を及ぼすことはまずないと思われるので気にしなくて大丈夫です((気になるようであれば corret 引数に FALSE を指定することで補正を行わなくすることも可能です))。
"X-squared = 26.331, df = 1, p-value = 1.438e-07" はカイ二乗検定のカイ二乗値などを表しています。この中で重要なのが p-value です。この値は、pattern a の CVR と pattern b の CVR が等しい場合に今回のような結果が得られる確率を表しています。CVR が同じである確率と解釈しても A/B テストの検証においては差し支えないでしょう((A であれば B が得られる確率が p-value であって、B が得られたから A である確率ではないことを強調しておきます))。CVR が同じである確率と捉えた場合、同じである確率は `r signif(test_result$p.value, 3)` と極めて低いと言えます。
あとの出力内容は片側検定であることや、95 % 信頼区間の情報やそれぞれの CVR を表しています。
以上の結果より、今回の結果は有意な差と言えそうです。
ちなみに、p-value の目安として 0.05 を切ったら有意な差（5 % 有意水準）とすることが多いですが、経営判断する上では CVR 以外の様々な要因（e.g. メンテナンスコスト）もあるので、あくまで参考程度にするのが良いと思います。
### ユーザの定着率（定着数）を出す
新しい LP によって CVR が上がったことを確認できましたが、ユーザの定着率が下がってしまっては意味がありません。
ということで、次は定着率を出してみます。
定着率の定義として、ここではユーザ登録したユーザが n 日後にアクセスした割合を採用します((個人的に、普段の分析では n 週間後にアクセスした割合を採用しています。ソーシャルゲームのように毎日アクセスすることを前提としたサービスでは n 日後にアクセスした割合で良いと思います))。
少し複雑な処理になるので、SQL の書き方にも様々な選択肢がありますが、例えば次のようなクエリで各パターンの経過日数ごとの UU と定着率を算出できます（PostgreSQL）。
```sql
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 では次のような手順を踏んでいます。
1. access_logs を user_id でグループ化する
2. グループごとに最小の unix time を算出する
3. 2 で算出した時間と time を日に変換してから差分を算出して elapsed_days とする
4. event_logs と user_id で inner join する
5. event が 'cv' で elapsed_days が 14 より小さいレコードに限定する
6. pattern, elapsed_days でグループ化する
7. グループごとにユニークな user_id をカウントする
8. pattern でグループ化する
9. グループごとに elapsed_days が最小のレコードの uu の値で uu を割る
10. pattern, elapsed_days でソートする
これに対応する dplyr の処理は次のようになります。
```{r}
# 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 条件で評価し、その結果をベクトルで返す演算子です。なお、ベクトルの論理和には `|` 演算子を使います。
```{r}
# Ruby の (1..5).to_a に相当
x <- 1:5
# 2 より大きい要素は TRUE になる
x > 2
# 4 より小さい要素は TRUE になる
x < 4
# ベクトルの論理積。2 より大きく 4 より小さい要素は TRUE になる
x > 2 & x < 4
```
`filter` 関数は指定したベクトルの要素が TRUE の場合に対応する行のデータを抽出します。上記の例では event の値が "cv" で elapsed_days が 14 より小さいデータを抽出しています。
ちなみに、上記の処理をもう少し最適化すると次のようになります。不要なデータを早い段階で削除しています。
```r
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)
```
さて、本題に戻って結果を見てみましょう。
```{r}
# print 関数の n 引数に表示する行数を指定することで全データを表示
print(uu_transitions, n = nrow(uu_transitions))
```
初日の uu は pattern a の方が大きいですが、全体的に pattern b の定着率は pattern a に比べて低く、その結果 14 日目 (elapsed_days が 13) の uu は pattern b の方が小さくなっています。
値の推移については数字を眺めるよりも可視化した方がわかりやすいでしょう。可視化には [ggplot2](http://ggplot2.org/) パッケージを使うことにします。
ggplot2 は折れ線グラフ、面グラフ、ヒストグラム等を統一的なインタフェースで描画できるパッケージで、複数のグラフの重ね合わせや軸の分割など複雑なグラフも手軽に描画できることが特徴です。
`plot` 関数等を使ったグラフの描画に慣れている方にとって ggplot2 は使いにくいかもしれませんが、R を初めて使うエンジニアの方には ggplot2 を使った描画方法を習得することをお勧めします。
```{r}
library(ggplot2)
# グラフにデータ (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 の使い方については[ドキュメメント](http://docs.ggplot2.org/current/)に豊富な例が載っているのでそちらを参照してください。次の資料も ggplot2 を使う上で非常に参考になります。
[一粒で3回おいしいggplot2](http://www.slideshare.net/syou6162/tsukuba)
さて、グラフを見ると pattern b は 2 日目以降 pattern a よりも UU が低いことがひと目でわかります。これは LP の内容からユーザが期待する内容とサービスの内容に差異があったと考えられます。
以上を踏まえると、現状のまま新しい LP を導入するのは避けた方が良いという判断になるでしょう。
## 番外編
### R でコホート分析
最近 Google Analytics にコホート分析の機能が追加されましたね。R であのような表を作成してみましょう。縦軸にユーザ登録日、横軸に経過日数、値に定着率を取ることにします。
データは先程も利用した `access_log_df` を使います。今回はユーザ登録日を付与していることと、レコード数を絞るためにログの期間を 4/7 までに限定していること以外は A/B テストのユーザ定着率を出した時とほとんど変わりません。
```{r}
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
```
これをおなじみ `tidyr::spread` 関数で横長のデータに変換します。
```{r}
cohort %>%
select(registraion_date, elapsed_days, retention_ratio) %>%
tidyr::spread(elapsed_days, retention_ratio) %>%
arrange(registraion_date)
```
簡単ですね！
## 最後に
以上、単純な例ではありましたが、A/B テストの結果を検証をする際にどのような検証を行うか、R を使う場合にどうすれば良いかがなんとなく理解いただけたのではないかと思います。
もちろん、より単純な検証を求められるケース、より詳細な検証を求められるケース、全く別の検証を求められるケースなど状況によって要求は様々なので、状況に応じて検証方法を検討してください。
今回は単純な集計処理しか行いませんでしたが、これから R で本格的に分析をしようと思っている方は以下の資料にも目を通すことをお勧めします。
Tokyo.R のような R コミュニティにも一度顔を出してみると良いと思います。
* [Rデータ自由自在](http://www.amazon.co.jp/gp/product/4621061542)
- R の基本的なデータ処理についてコンパクトにまとめられているので個人的にオススメです
* [同志社大学の金先生のページ](http://mjin.doshisha.ac.jp/R/)
- 統計学や機械学習の手法を R で行う方法について非常に参考になります（かなりの量なのでこういうページがあるということだけでも覚えておくと良いでしょう）
* [dplyr の vignettes](http://cran.r-project.org/web/packages/dplyr/vignettes/)
- dplyr を使いこなしたい方は一読すべきでしょう
* [ハンバーガー統計学にようこそ！](http://kogolab.chillout.jp/elearn/hamburger/)
- R とは関係ないですが検定についてわかりやすく解説されています
* [アイスクリーム統計学にようこそ！](http://kogolab.chillout.jp/elearn/icecream/)
- 「ハンバーガー統計学にようこそ！」と合わせて読むと良いでしょう
あとは拙作のマニアックなスライドですがいくつか参考になるかもしれません。
* [Rデバッグあれこれ](http://www.slideshare.net/abicky/r-9034336)
- 今回触れたデバッグ方法よりも詳細な情報を載せています（ただ、情報が古い・・・）
* [Rデータフレーム自由自在 〜基本操作からピボットテーブルまで〜](http://www.slideshare.net/abicky/r-10128090)
- dplyr を使う場合は見る価値のないスライドです。dplyr がないとどのような処理になるか知りたい方に。
* [Rのスコープとフレームと環境と 〜parent.env と parent.frame の違い〜](http://www.slideshare.net/abicky/r-11657813)
- デバッグ方法のところで軽く触れましたがフレームなどの話について触れています
* [Rのデータ構造とメモリ管理](http://www.slideshare.net/abicky/r-15350295)
ではでは快適な R ライフを！！
※このエントリーは [knitr](http://yihui.name/knitr/) パッケージを使って[作成しました](https://gist.github.com/abicky/d3deb617bddda461e2eb)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment