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

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

カイ二乗検定のメタアナリシスをやってみた(階層ベイズでも試してみた追記あり&タイトル変更済み)

記事タイトルに反して僕は実は統計的検定が大嫌いなんですが、皆さんいかがお過ごしでしょうか(笑)。ということで、今回はややマニアックなメタアナリシスの話題でもしてみようかと思います。「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


f:id:TJO:20160218213011p:plain


さて、この動線改善には本当にCVRを向上させる効果があったんでしょうか? 各々についてCVRを見てみると、以下のようになっています。


f:id:TJO:20160218223625p:plain


これを見る限りでは、常に「動線改善あり」の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統計量 z_kを用いて*4


 p_{group} = \frac{\displaystyle \sum^n_{k = 1} z_k}{\sqrt{n}}


と表すことができると述べられています。このように異なる"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で以下のようにヒントをもらったので、それに合わせてみました。


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)

f:id:TJO:20160219154425p:plain
f:id:TJO:20160219154439p:plain


とりあえずbの値は0からは正の方向に十分に離れた分布を描いているということで、改善効果ありと結論付けて良さそうです。ただ、個々のサイトの個体差clについては何とも言いようのない感じが。。。heterogeneityの観点から言うとちょっと微妙かも。

*1:最近だと機械学習でやるところが増えてますね、実際に

*2:何となくコンサルっぽいですね笑

*3:つまりrepated measure two-way ANOVAみたいな、というかそもそもone-way ANOVAすら使えない

*4:もちろん正規分布に従うという仮定があるので、それに沿わない場合は使えない

*5:もっともこれもp-value hackingの一種なのかもしれませんが

*6:これは対立仮説が「heterogeneityあり」なので有意でない方が良い検定