読者です 読者をやめる 読者になる 読者になる

A/B テストで施策の効果を検証!エンジニアのための R 入門

こんにちは、買物情報事業部でサーバサイドの開発を担当している荒引 (@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.md.tsukuba.ac.jp")' >> ~/.Rprofile

.Renviron は R を実行する上での環境変数を定義するファイルです。

.Rprofile は起動時の処理を記述するファイルです。options 関数の repos 引数に筑波大学のリポジトリを指定しています。 ただ、筑波大学は休日に停電を行うことがあるので、その際は別のリポジトリを指定しないとパッケージをインストールできません。

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_fileaccess_log_file にそれぞれの絶対パスが格納されています。

それでは、eventlogs.tsv と access_logs.tsv を読み込みましょう。 小規模なデータの読み込みには read.table 関数などを使うのが一般的ですが、read.table 系の関数は非常に遅いので data.table パッケージの fread 関数を使います*2。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 に近い感覚で処理を記述できるのでエンジニアにとっては馴染みやすいでしょう。*3

まず、読み込んだデータを tbl_df 関数で dplyr 用の data.frame 表現に変換します。data.table オブジェクトのままでも処理できますが、次の理由から変換することにします。

  1. data.table オブジェクトは POSIXlt オブジェクトを扱えないなどの制約がある
  2. 私の手元の環境では tbl_df で変換した方が若干高速だった
  3. 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 では次のような手順を踏んでいます。

  1. pattern, event でグループ化する
  2. それぞれのグループに関してユニークな 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 パッケージが提供している演算子です*4。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) を行った旨が表示されていますが、連続性の補正の有無で経営判断が変わるほど結果に影響を及ぼすことはまずないと思われるので気にしなくて大丈夫です*5

"X-squared = 26.331, df = 1, p-value = 1.438e-07" はカイ二乗検定のカイ二乗値などを表しています。この中で重要なのが p-value です。この値は、pattern a の CVR と pattern b の CVR が等しい場合に今回のような結果が得られる確率を表しています。CVR が同じである確率と解釈しても A/B テストの検証においては差し支えないでしょう*6。CVR が同じである確率と捉えた場合、同じである確率は 1.44 × 10-7 と極めて低いと言えます。

あとの出力内容は片側検定であることや、95 % 信頼区間の情報やそれぞれの CVR を表しています。

以上の結果より、今回の結果は有意な差と言えそうです。

ちなみに、p-value の目安として 0.05 を切ったら有意な差(5 % 有意水準)とすることが多いですが、経営判断する上では CVR 以外の様々な要因(e.g. メンテナンスコスト)もあるので、あくまで参考程度にするのが良いと思います。

ユーザの定着率(定着数)を出す

新しい LP によって CVR が上がったことを確認できましたが、ユーザの定着率が下がってしまっては意味がありません。 ということで、次は定着率を出してみます。 定着率の定義として、ここではユーザ登録したユーザが n 日後にアクセスした割合を採用します*7

少し複雑な処理になるので、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 では次のような手順を踏んでいます。

  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 の処理は次のようになります。

> # 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))

f:id:a_bicky:20150508114152p:plain

上記の例では pattern によって色を変えて描画していますが、例えばデモグラ情報(年代、性別等の情報)を持ったデータであれば、年代ごとに異なる色を使い、性別ごとに異なる線種を使うようなことも可能です。 ggplot2 の使い方についてはドキュメメントに豊富な例が載っているのでそちらを参照してください。次の資料も ggplot2 を使う上で非常に参考になります。

一粒で3回おいしい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 ライフを!!

※このエントリーは knitr パッケージを使って作成しました

*1:コーディングスタイルに迷ったら R 界で影響力の大きい Hadley Wickham 氏のコーディングスタイルに従うのが良いと思います

*2:手元の環境だと、アクセスログを読み込むのに read.table は data.table::fread の 100 倍近く時間がかかります

*3:まだ枯れていないのでバグを踏むことはありますが・・・

*4:厳密には magrittr パッケージが提供しているものを dplyr パッケージが取り込んでいます

*5:気になるようであれば correct 引数に FALSE を指定することで補正を行わなくすることも可能です

*6:A であれば B が得られる確率が p-value であって、B が得られたから A である確率ではないことを強調しておきます

*7:個人的に、普段の分析では n 週間後にアクセスした割合を採用しています。ソーシャルゲームのように毎日アクセスすることを前提としたサービスでは n 日後にアクセスした割合で良いと思います

/* */ @import "/css/theme/report/report.css"; /* */ /* */ body{ background-image: url('http://cdn-ak.f.st-hatena.com/images/fotolife/c/cookpadtech/20140527/20140527163350.png'); background-repeat: repeat-x; background-color:transparent; background-attachment: scroll; background-position: left top;} /* */ body{ border-top: 3px solid orange; color: #3c3c3c; font-family: 'Helvetica Neue', Helvetica, 'ヒラギノ角ゴ Pro W3', 'Hiragino Kaku Gothic Pro', Meiryo, Osaka, 'MS Pゴシック', sans-serif; line-height: 1.8; font-size: 16px; } a { text-decoration: underline; color: #693e1c; } a:hover { color: #80400e; text-decoration: underline; } .entry-title a{ color: rgb(176, 108, 28); cursor: auto; display: inline; font-family: 'Helvetica Neue', Helvetica, 'ヒラギノ角ゴ Pro W3', 'Hiragino Kaku Gothic Pro', Meiryo, Osaka, 'MS Pゴシック', sans-serif; font-size: 30px; font-weight: bold; height: auto; line-height: 40.5px; text-decoration: underline solid rgb(176, 108, 28); width: auto; line-height: 1.35; } .date a { color: #9b8b6c; font-size: 14px; text-decoration: none; font-weight: normal; } .urllist-title-link { font-size: 14px; } /* Recent Entries */ .recent-entries a{ color: #693e1c; } .recent-entries a:visited { color: #4d2200; text-decoration: none; } .hatena-module-recent-entries li { padding-bottom: 8px; border-bottom-width: 0px; } /*Widget*/ .hatena-module-body li { list-style-type: circle; } .hatena-module-body a{ text-decoration: none; } .hatena-module-body a:hover{ text-decoration: underline; } /* Widget name */ .hatena-module-title, .hatena-module-title a{ color: #b06c1c; margin-top: 20px; margin-bottom: 7px; } /* work frame*/ #container { width: 970px; text-align: center; margin: 0 auto; background: transparent; padding: 0 30px; } #wrapper { float: left; overflow: hidden; width: 660px; } #box2 { width: 240px; float: right; font-size: 14px; word-wrap: break-word; } /*#blog-title-inner{*/ /*margin-top: 3px;*/ /*height: 125px;*/ /*background-position: left 0px;*/ /*}*/ /*.header-image-only #blog-title-inner {*/ /*background-repeat: no-repeat;*/ /*position: relative;*/ /*height: 200px;*/ /*display: none;*/ /*}*/ /*#blog-title {*/ /*margin-top: 3px;*/ /*height: 125px;*/ /*background-image: url('http://cdn-ak.f.st-hatena.com/images/fotolife/c/cookpadtech/20140527/20140527172848.png');*/ /*background-repeat: no-repeat;*/ /*background-position: left 0px;*/ /*}*/