記事タイトルに反して僕は実は統計的検定が大嫌いなんですが、皆さんいかがお過ごしでしょうか（笑）。ということで、今回はややマニアックなメタアナリシスの話題でもしてみようかと思います。「t-testのメタアナリシス」みたいな、いわゆるRosenthal's methodとして知られているような一般的な話題にもちょっとだけ触れています。
今回想定しているシチュエーション
例えばですが、「あるパーソナライズされた*1動線改善を4つの異なるeコマースサイトで行った」みたいなケースを想定してみます。それぞれのeコマースサイトは独立した事業者であるため、統一したアクセス解析はできないものと仮定します。で、動線改善は「お気に入り→購入」のところに施したとして、ここでA/Bテストを行ったと仮定します。この時、アクセス解析の結果からいわゆる2×2分割表として
サイトA
| 非CV | CV | |
|---|---|---|
| 改善なし | 250 | 125 |
| 改善あり | 300 | 180 |
サイトB
| 非CV | CV | |
|---|---|---|
| 改善なし | 1500 | 650 |
| 改善あり | 1200 | 620 |
サイトC
| 非CV | CV | |
|---|---|---|
| 改善なし | 600 | 300 |
| 改善あり | 550 | 310 |
サイトD
| 非CV | CV | |
|---|---|---|
| 改善なし | 2200 | 1050 |
| 改善あり | 1300 | 700 |
という集計データが得られたものとします。ついでなので、積み上げ棒プロットで図示してみるとこうなります*2。
さて、この動線改善には本当にCVRを向上させる効果があったんでしょうか？ 各々についてCVRを見てみると、以下のようになっています。
これを見る限りでは、常に「動線改善あり」のCVRの方が「なし」を上回っているようです。それゆえ、直感的にはこの動線改善には4つのサイトをまたいで共通したCVR向上効果があるように思われます。しかしながら、既に集計された(tabulated)データゆえ「8通り」「4組」のサンプルしか手元になく、どう見ても検定力が弱い≒サンプルサイズに物を言わせた分析方法は使えません*3。では、どうすれば良いのでしょうか？
まずはR上でデータを整理して、個別にカイ二乗検定してみる
右往左往していても仕方ないので、まずはデータをR上で整理してみましょう。2×2分割表で表せるタイプのデータセットなので、普通にmatrix型にしておくのが良さそうです。
> d1<-matrix(c(250,125,300,180),ncol=2,byrow=T) > d2<-matrix(c(1500,650,1200,620),ncol=2,byrow=T) > d3<-matrix(c(600,300,550,310),ncol=2,byrow=T) > d4<-matrix(c(2200,1050,1300,700),ncol=2,byrow=T) > d1 [,1] [,2] [1,] 250 125 [2,] 300 180 > d2 [,1] [,2] [1,] 1500 650 [2,] 1200 620 > d3 [,1] [,2] [1,] 600 300 [2,] 550 310 > d4 [,1] [,2] [1,] 2200 1050 [2,] 1300 700
まぁここまでは簡単です。次に、一般に2×2分割表で表せるデータに対する効果検証の類は「独立性の検定」即ちカイ二乗検定で行うことができます。以前の記事でも簡単に取り上げたことがありますね。
ということで、取るものもとりあえず4サイトそれぞれについてカイ二乗検定をやってみましょう。
> chisq.test(d1) Pearson's Chi-squared test with Yates' continuity correction data: d1 X-squared = 1.4164, df = 1, p-value = 0.234 > chisq.test(d2) Pearson's Chi-squared test with Yates' continuity correction data: d2 X-squared = 6.4822, df = 1, p-value = 0.0109 > chisq.test(d3) Pearson's Chi-squared test with Yates' continuity correction data: d3 X-squared = 1.3122, df = 1, p-value = 0.252 > chisq.test(d4) Pearson's Chi-squared test with Yates' continuity correction data: d4 X-squared = 3.9182, df = 1, p-value = 0.04777
有意水準p < 0.05とした場合、サイトAは有意でない、Bは有意、Cは有意でない、Dは有意という結果になりました・・・のですが。これでは「4つのサイトそれぞれにとってパーソナライズされた動線改善が有効だったか否か」は分かっても、全体として本当に動線改善の効果があったかどうかは何とも言えません。しかも4サイト中2つでは有意差ありで2つは有意差なしということで、五分五分。これでは結論が出せないことになってしまいます。直感的な印象とは異なり、これはかなりややこしい状況です。
では、どうやってここから「4つのサイトで共通した動線改善の効果の有無」についての蓋然性ある結論を得たら良いのでしょうか？
個々の"study"を統合したいならメタアナリシスで
これは言い換えると「互いに異なるデータセットを持つ個々の閉じたケース（これを"study"と呼ぶ）内で共通した実験的介入を行った」時の実験結果同士をつなぎ合わせるのと同じ取り組みであり、実はそのようなシチュエーションは実験心理学や認知神経科学の実験研究では多く見られることもあってかなり以前からその方法論が整備されています。
例えば、かつて他ならぬ僕自身が某所で出版した論文でも引用したRosenthal (1992)ではn個のstudyが何かしらの統計学的検定によりそれぞれ異なるp値を持つ場合、全体のp値は各々のp値に対応するz統計量を用いて*4
と表すことができると述べられています。このように異なる"study"同士の統計量を統合して全体的な傾向を分析することをメタアナリシスと呼ぶわけですが、例えばこちらのサイトでは様々な統計量を統合する方法が紹介されています。で、ここではどう見ても「異なる個々のカイ事情検定の結果をメタアナリシス的に統合する」手法が必要なわけです。そんなものは本当にあるんでしょうか？
カイ二乗検定のメタアナリシス手法の一つ、Cochran-Mantel-Haenszel testを実践してみる
結論から言うとそんなものは本当にありまして（笑）、その名もCochran-Mantel-Haenszel testと言います。しかもご丁寧にこちらのサイトにその方法がきちんと解説されています。この手法、疫学分野では結構有名らしいのですが、浅学な僕は今まで全く知りませんでしたorz 色々な人に聞いて回った結果たどり着いたのがこちらの手法でして。。。ということで、それをそのままなぞってRでコードを書いてみました。その結果が以下。
> e1<-function(x){ + x[1,1]-(x[1,1]+x[1,2])*(x[1,1]+x[2,1])/sum(x) + } > e2<-function(x){ + ((x[1,1]+x[1,2])*(x[1,1]+x[2,1])*(x[1,2]+x[2,2])*(x[2,1]+x[2,2]))/(sum(x)^3-sum(x)^2) + } > chi2<-((e1(d1)+e1(d2)+e1(d3)+e1(d4))-0.5)^2/(e2(d1)+e2(d2)+e2(d3)+e2(d4)) > 1-pchisq(chi2,1) [1] 0.0002988727 # p値
ということで、有意水準をp < 0.05に置く限りは「4サイト共通(collapsed across 4 web sites)のパーソナライズされた動線改善の効果は有意に見られる」と言えることになりました。これは直感的な印象にも合いますし、それなりに蓋然性ある結果と言えるでしょう*5。
ちなみにこれだけやっておいて何ですが、実はCRANにはちゃんと{metafor}というメタアナリシスに特化したパッケージがあって、それを使っても簡単にCMH testの計算ができます。具体的にはrma.mh関数で算出します。
> install.packages('metafor') > library(metafor) > rma.mh(c(d1[1,1],d2[1,1],d3[1,1],d4[1,1]),c(d1[1,2],d2[1,2],d3[1,2],d4[1,2]),c(d1[2,1],d2[2,1],d3[2,1],d4[2,1]),c(d1[2,2],d2[2,2],d3[2,2],d4[2,2])) Fixed-Effects Model (k = 4) Test for Heterogeneity: Q(df = 3) = 0.4987, p-val = 0.9192 Model Results (log scale): estimate se zval pval ci.lb ci.ub 0.1437 0.0395 3.6361 0.0003 0.0662 0.2212 Model Results (OR scale): estimate ci.lb ci.ub 1.1546 1.0685 1.2476 Cochran-Mantel-Haenszel Test: CMH = 13.0774, df = 1, p-val = 0.0003 Tarone's Test for Heterogeneity: X^2 = 0.4987, df = 3, p-val = 0.9192
同じ手法を用いているので当然ですが、同じp値が得られて「動線改善には有意に効果あり」という結論になりました。なおrma.mh{metafor}関数だと、そもそもこのようなメタアナリシスを行って良いかどうかの指標にもなるheterogeneityのチェックもできるので便利です*6。
また、escalc関数を使えば個々のデータセットについて例えばオッズ比や対数相対リスクと言った指標のeffect sizeを算出した上でさらにオブジェクトとして格納することができます。これは例えばCMH testなどのさらなるメタアナリシスに転用することが可能です。そんなわけで、今後は{metafor}パッケージが有効そうな事例が出てきたら積極的に使ってみようかなと思う次第です。
おまけ
どんな方法論があるのか？と色々な人に聞いて回ったら、友人のBayesianが「だったら共役事前分布としてベータ分布突っ込んで4サイト×2条件間で共通の真のCVR（＝二項分布のパラメータ）を求めるMCMCサンプリングやればいいんじゃね」みたいな力技を薦めてきたんですが、まだやってません（笑）。そのうちやってみます。。。
追記1
分母つきロジスティック回帰でやる方法を思い出したので書いておきます。
> dat.hb<-cbind(rbind(d1[1,],d1[2,],d2[1,],d2[2,],d3[1,],d3[2,],d4[1,],d4[2,]),rep(c(0,1),4),c(rep(1,2),rep(2,2),rep(3,2),rep(4,2))) > dat.hb<-as.data.frame(dat.hb) > names(dat.hb)<-c('ncv','cv','intervention','site') > dat.hb ncv cv intervention site 1 250 125 0 1 2 300 180 1 1 3 1500 650 0 2 4 1200 620 1 2 5 600 300 0 3 6 550 310 1 3 7 2200 1050 0 4 8 1300 700 1 4 > dat.hb$intervention<-as.factor(dat.hb$intervention) > dat.hb$site<-as.factor(dat.hb$site) > dat.hb.glm<-glm(cbind(cv,ncv)~intervention+site,dat.hb,family=binomial) > summary(dat.hb.glm) Call: glm(formula = cbind(cv, ncv) ~ intervention + site, family = binomial, data = dat.hb) Deviance Residuals: 1 2 3 4 5 6 7 8 -0.2024 0.1749 -0.3245 0.3431 0.1683 -0.1686 0.2403 -0.2996 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.67100 0.07492 -8.956 < 2e-16 *** intervention1 0.14368 0.03952 3.636 0.000277 *** site2 -0.15002 0.07923 -1.893 0.058299 . site3 -0.03405 0.08731 -0.390 0.696514 site4 -0.07768 0.07753 -1.002 0.316381 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 20.37475 on 7 degrees of freedom Residual deviance: 0.49884 on 3 degrees of freedom AIC: 69.757 Number of Fisher Scoring iterations: 3
改善あり(intervention1)の偏回帰係数が正でかつ有意に0より大きいということで、改善の効果ありという結果になっています。
追記1の追記
さらに前職の友人でもあるO君から「偏回帰係数のみでやると色々危ないので尤度比検定でやるべきでは」という的確なコメントが後日ありました。はい、尤度比検定とか真面目に普段やってないので完全に失念してましたごめんなさい。。。ということでやってみました。
> summary(dat.hb.glm2) Call: glm(formula = cbind(cv, ncv) ~ site, family = binomial, data = dat.hb) Deviance Residuals: 1 2 3 4 5 6 7 8 -0.9505 0.8328 -1.7558 1.8883 -0.8383 0.8524 -1.2436 1.5748 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.58961 0.07139 -8.259 <2e-16 *** site2 -0.16463 0.07909 -2.082 0.0374 * site3 -0.04445 0.08721 -0.510 0.6103 site4 -0.10354 0.07716 -1.342 0.1796 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 20.375 on 7 degrees of freedom Residual deviance: 13.701 on 4 degrees of freedom AIC: 80.96 Number of Fisher Scoring iterations: 3 > anova(dat.hb.glm,dat.hb.glm2,test='Chisq') Analysis of Deviance Table Model 1: cbind(cv, ncv) ~ intervention + site Model 2: cbind(cv, ncv) ~ site Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 3 0.4988 2 4 13.7013 -1 -13.203 0.0002796 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
やり方はみどりぼん5章pp.101-102を参考にしました。ということで改善あり(intervention1)の効果はあったという結論で良さそうです。
追記2
ついでに色々調べながらStanで階層ベイズロジスティック回帰をやってみました。どういうモデルにするかは僕と同じく岩波DS刊行委員を務める友人の松浦さんからTwitterで以下のようにヒントをもらったので、それに合わせてみました。
@TJO_datasci ブログの記事に関して言えば、カウント数は二項分布に従うとして、その確率pは単純にはp <- inv_logit(a + b*改善 + サイト差)で決めます。そして、サイト差を階層事前分布にするのでOKかと。
— K (@hankagosa) 2016年2月19日
Stanコードは以下の通り。とりあえず"binom_hb_gen.stan"という名前で保存しておきます。N/2とかしてるせいでパーサから怒られると思いますが、一応走ります。
data {
int<lower=0> N;
int<lower=0> ncv[N];
int<lower=0> cv[N];
real<lower=0,upper=1> iv[N];
}
parameters {
real cl[N/2];
real s;
real a;
real b;
}
model {
real p[N];
s~uniform(0,1e4);
cl~normal(0,s);
for (i in 1:N/2){
p[2*i-1]<-inv_logit(a+b*iv[2*i-1]+cl[i]);
p[2*i]<-inv_logit(a+b*iv[2*i]+cl[i]);
}
for (i in 1:N)
cv[i]~binomial(ncv[i]+cv[i],p[i]);
}
個々のサイトの影響clを階層事前分布で与えることでバランスさせた上で、介入の効果の有無を二値説明変数ivにかかる係数bの値から評価しようという形です（これは上の普通のロジスティック回帰と同じ発想）。これをkickするRコードは以下の通り。
> hb.dat $N [1] 8 $ncv [1] 250 300 1500 1200 600 550 2200 1300 $cv [1] 125 180 650 620 300 310 1050 700 $iv [1] 0 1 0 1 0 1 0 1 > hb.fit<-stan(file='binom_hb_gen.stan',data=hb.dat,iter=2000,chains=4) > hb.fit Inference for Stan model: binom_hb_gen. 4 chains, each with iter=2000; warmup=1000; thin=1; post-warmup draws per chain=1000, total post-warmup draws=4000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat cl[1] 0.05 0.01 0.12 -0.11 0.00 0.03 0.07 0.34 149 1.01 cl[2] -0.04 0.01 0.11 -0.24 -0.09 -0.05 -0.01 0.19 147 1.02 cl[3] 0.04 0.01 0.12 -0.13 -0.01 0.02 0.06 0.30 140 1.01 cl[4] 0.01 0.01 0.11 -0.16 -0.03 0.00 0.03 0.27 138 1.01 s 0.12 0.01 0.17 0.01 0.04 0.08 0.13 0.52 125 1.03 a -0.76 0.01 0.11 -1.00 -0.79 -0.75 -0.72 -0.59 135 1.02 b 0.15 0.00 0.04 0.07 0.12 0.15 0.17 0.22 473 1.01 lp__ -7509.93 0.29 3.50 -7517.81 -7511.93 -7509.65 -7507.50 -7504.17 142 1.02 Samples were drawn using NUTS(diag_e) at Fri Feb 19 15:40:42 2016. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1). > hb.fit.coda<-mcmc.list(lapply(1:ncol(hb.fit),function(x) mcmc(as.array(hb.fit)[,x,]))) > plot(hb.fit.coda)
とりあえずbの値は0からは正の方向に十分に離れた分布を描いているということで、改善効果ありと結論付けて良さそうです。ただ、個々のサイトの個体差clについては何とも言いようのない感じが。。。heterogeneityの観点から言うとちょっと微妙かも。
