渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

サンプリング時の最適なサンプルサイズをRパッケージ{pwr}で求める

最近、「ビッグデータ」というバズワードに対するアンチテーゼとして叫ばれるようになってきたのが、


統計学ってのは限られたサンプル(抽出)データから、まだ見ぬ全体像を知るためのもの」「だからビッグデータなんて苦労して集める必要はない、サンプリングされたデータだけで十分だ」


という主張。えーと、半分はその通りだと思います。けれども、半分はそうでもないかなぁ、と。


何故なら、レコメンダーとかSPAMフィルタなどのバックエンドシステム開発では、できれば全数データを使って可能な限り精度を上げ続けた方が良いものが多いからです。だからHadoop以下大規模分散処理などの高度な手法を沢山駆使しているわけで、そこでは依然として全数データは非常に重要です。


一方、マーケティングなどでオフライン&アドホックで分析する際には、そこまでやらんでもええやん的な状況は多くあります。手動でデータ分析したいんだけど、全数データ取ってきたらRやPythonがout of memoryで落ちるとか。。。もっとひどい場合は、肝心のデータをHDFSからHiveで抽出しようとしたら、HDFS上でOOM&全力でdump吐いて落ちるとかorz 何でオレこんな重たいデータ相手にしとるんじゃあああみたいな場面は結構ありがちです。


ということで、そういう場合はサンプリングすれば良いじゃないかというのが僕の意見。むしろサンプリングすることで分析の負荷を下げて、代わりに分析の方法やバリエーションを増やしてやれば良いわけです。やり方は色々あって、ユーザーIDを一定のルールに従って無作為抽出するとか、それこそHiveならTABLESAMPLE()でレコード単位で取ってくるとかやっても良いわけですが、最大の問題はサンプルサイズ(標本の個数)をどう決めるか?というところ。


そこで、今回はそのサンプルサイズを決める上で役に立つCRANパッケージ{pwr}を紹介してみます。もう名前の通りで、effect size統計学的に見たその現象の「強さ」)と検出力(個々の統計学的検定が持つ「感度」)に従って最低限必要なサンプルサイズを決めてくれるパッケージです。以下のリポジトリと、解説ページが参考になると思います。


なお、サンプルサイズを決める方法論については『サンプルサイズの決め方』(永田靖・朝倉書店)が最も参考になります。というか、多分日本語でサンプルサイズの決定法を網羅的に「全部」解説している本はこれ一冊しかないかもです。出来る限り、この本でサンプルサイズに関する統計学的原理を(斜め読みでも良いので)網羅的に学習することを強く推奨します。というか、多分勉強しないとこの後の説明を読んでも何を言ってるか分からないかもです。


サンプルサイズの決め方 (統計ライブラリー)

サンプルサイズの決め方 (統計ライブラリー)


ここでは、Webデータ分析とかマーケティングリサーチとかでありがちな、「n1万人の母集団における○○率・○○と××との相関・○○の回帰モデルをサンプリングしてn2人のデータから求めたい」というケースに必要なものをピックアップして紹介します。その他のポイントについては割愛、ということで。


そうそう、effect sizeとstatistical power(検出力)についての説明も今回は割愛します。大事な概念なので、時間のある時に勉強してみて下さい。このパッケージの作者は、Jacob Cohenが提唱するsmall / medium / large effect sizeという概念を参考にしているので、もし興味のある人はその原典を読んでも良いかもです。


Statistical Power Analysis for the Behavioral Sciences

Statistical Power Analysis for the Behavioral Sciences


以下「一般にeffect sizeはこちらの期待に反してものすごーく小さい」と仮定して常にCohen曰くのsmall effect sizeのケースに固定し、statistical powerも常に業界標準値の0.9に固定して、{pwr}パッケージで計算していきます。

{pwr}パッケージのpwr.XYZ.test()関数共通のポイント


知りたい変数まで含めて全て関数の引数として与えて、知りたい変数だけx = NULLのようにNULL指定しておく、というスタイルを取っています。返値ではなく、コマンドラインへの標準出力でしか結果が返ってこないので要注意。


母不良率(○○率)を求めたい場合


「不良率」というのは統計学や品質工学あたりでのタームで、要するに○○率を求めたいというケース。例えば、全体会員数2500万人のサービスで、これまでにソーシャルグラフ内で半分以上のフレンドに対して「いいね!」した人は何%いたか?が知りたい、みたいなパターンを想定してみましょう。この場合は、pwr.2p2n.test()関数を使います。内容は以下の通り。

> pwr.2p2n.test(h=0.2, n1, n2=NULL, sig.level=0.05, power=0.9, alternative="two.sided")
# h: effect size
# n1: 母集団のサイズ
# n2: (ここでは)求めたいサンプルサイズ
# sig.level: 有意水準(0.05で固定)
# power: 検出力(今回は0.9で固定)
# alternative: 検定方法(今回はtwo.sidedで固定)


では実際に、n1を2500万人としてやってみましょう。Cohenはこのケースでのsmall effect sizeは0.2と書いているので、h = 0.2とします。

> pwr.2p2n.test(h=0.2,n1=2.5e8,n2=NULL,sig.level=0.05,power=0.9,alternative="two.sided")

     difference of proportion power calculation for binomial distribution (arcsine transformation) 

              h = 0.2
             n1 = 2.5e+08
             n2 = 262.6858 # これが求めるサンプルサイズ
      sig.level = 0.05
          power = 0.9
    alternative = two.sided

 NOTE: different sample sizes 


最低でも263人のサンプルをランダム抽出してくるべきだ、という結果になりました。これを例えば日本人全員とかに拡大してみると、

> pwr.2p2n.test(h=0.2,n1=1.2e9,n2=NULL,sig.level=0.05,power=0.9,alternative="two.sided")

     difference of proportion power calculation for binomial distribution (arcsine transformation) 

              h = 0.2
             n1 = 1.2e+09
             n2 = 262.6855
      sig.level = 0.05
          power = 0.9
    alternative = two.sided

 NOTE: different sample sizes 


実はほぼ同じ値になります*1。「○○率が知りたければ263人以上!」と覚えればほぼ間違いないかも。


相関係数を求めたい場合


最近は2つの異なるユーザー行動指標間(n日連続ログイン×課金とか)の相関係数を求めたくなるケースも多いと思います。こういう場合はpwr.r.test()関数を使います。

> pwr.r.test(n=NULL, r=0.1, sig.level=0.05, power=0.9, alternative="two.sided")
# n: (ここでは)求めたいサンプルサイズ
# r: 相関係数(ただしこの値はeffect size次第で決めるべきもの)
# sig.level: 有意水準
# power: 検出力
# alternative: 検定方法


相関係数の場合は母集団のサイズは関係ないので、nは1個だけNULLとして引数で与えます。Cohenはsmall effect sizeでの値を見るにはr = 0.1にすべしと書いているので、これを踏襲するとこうなります。

> pwr.r.test(n=NULL,r=0.1,sig.level=0.05,power=0.9,alternative="two.sided")

     approximate correlation power calculation (arctangh transformation) 

              n = 1046.423 # 求めるサンプルサイズ
              r = 0.1
      sig.level = 0.05
          power = 0.9
    alternative = two.sided


最低でも1047人のサンプルを集める必要がある、という結果になりました。


分散分析や重回帰分析を行いたい場合


こいつらは結構なクセモノで、サンプルサイズではなく「群内観測数」「自由度」だけを計算してくれます。なので、ここから先は重回帰分析and/or分散分析の知識が必要です。分からない人はスキップして下さい。

> pwr.anova.test(k, n=NULL, f=0.1, sig.level=0.05, power=0.9)
# **分散分析の場合**
# k: 群の数
# n: (ここでは)求めたい群内観測数
# f: effect size
# sig.level: 有意水準
# power: 検出力

> pwr.f2.test(u, v=NULL, f2=0.02, sig.level=0.05, power=0.9)
# **重回帰分析の場合**
# u: 回帰の自由度
# v: (ここでは)求めたい残差の自由度
# f2: effect size
# sig.level: 有意水準
# power: 検出力


Cohenは分散分析に関してはf=0.1、重回帰分析に関してはf2=0.02にすべしと書いているので、そのまま踏襲して計算するとこうなります。

> pwr.anova.test(k=5, n=NULL, f=0.1, sig.level=0.05, power=0.9)

     Balanced one-way analysis of variance power calculation 

              k = 5
              n = 309.0514 # 群内観測数のサイズ
              f = 0.1
      sig.level = 0.05
          power = 0.9

 NOTE: n is number in each group 

> pwr.f2.test(u=5, v=NULL, f2=0.02, sig.level=0.05, power=0.9)

     Multiple regression power calculation 

              u = 5
              v = 823.0275 # 残差自由度のサイズ
             f2 = 0.02
      sig.level = 0.05
          power = 0.9


直感的には分散分析の方が分かりやすい*2気がするので、何となくですがpwr.anova.test()関数を使った方が良いかもですね。


おわりに


他にもt検定2種、カイ二乗検定、単一集団の母不良率、サンプルサイズが等しい2集団同士の母不良率の比較、正規分布の平均に関するサンプルサイズを求める関数がありますが、webデータ分析で使う場面が思い付かなかったので割愛しました。


ともあれ、例えばHDFSに膨大なデータは突っ込んであるんだけど、いざ全数データに対してHiveクエリ組んで○○率など計算しようとするとOOMになってしまうので無理。。。なんてケースでは、これに基づいて最低限必要なサンプルサイズを計算して、堂々とサンプリングして○○率などを計算すると良いのではないでしょうか。

*1:これは理論を勉強した上でモンテカルロ・シミュレーションやってみればすぐ分かります

*2:群内観測数なので。これが残差自由度だと回帰自由度で割ってから考えないといけない