因果推論の基本テクニック、回帰分析は何を意味するのか

因果推論のもっとも基本的なテクニックである回帰分析はよく知られていますが、モデルで仮定している前提や、変数の入れ方によって、結果の解釈が大きく違います。今回は、回帰分析の種類やモデルが意味することを書いてみたいと思います。 何気なく使われる回帰分析も、実はとても奥が深いのです。

因果推論の基本的な考えはこちらをご覧ください。 ちなみに因果推論にも流派がありますが、本ブログではルービンの因果モデルを元にしています。

www.medi-08-data-06.work

因果推論における回帰分析の使い方

まずは、因果推論における一般的な回帰分析を考えてみます。くどいほど使っている広告の有無と売り上げの例を使って、因果推論してみましょう。

データは200人のこんなデータを想定します。

CM	sex	sales
0	0	4642.379
1	1	6123.656
0	0	4530.731
1	0	4473.743
1	0	4781.420
0	0	5165.590
1	0	3992.895

変数は、広告を配信したかどうか(CM: 0 無し,1 有り)、性別(sex : 0 男,1 女)そして、それぞれの人の購入額(sales)です。

まずは広告の有無（X軸）と購買額（Y軸）をグラフにしてみると、こんな感じになります。

広告を配信した群ほど購入額が高い傾向に有りそうです。ここに単回帰分析を用いてモデルを作ってみましょう。式としては以下のようになります。

実際に求めてみると、図のようになりました。

回帰分析が意味すること

さて、この回帰式は何を意味しているのでしょうか。今回の目的は、広告の有無によって、購入額にどれぐらいの差が出るのかという、広告と購入額の因果関係を推論することが目的でした。

前回の記事にも書いたように、現実世界では広告の真の効果というものを推定することができません。そこで、性質が均一であり、唯一の違いは広告の有無のみであると仮定できる(交換可能性:Exchangebilityが成り立つと仮定できる)二つの群の平均値を比較して、その平均値の差を期待される広告の効果ということにしようというのが、ルービン先生の因果推論の考え方です。これを平均因果効果と呼びます。

先ほどの式では、CMは0or1なので、広告無しの群の購入額平均は、広告ありの群の購入額平均はとなり、広告の効果は傾きの係数となることが期待されるということです。少しややこしく書くと

=5336+498CM

=498

となります。縦棒は条件付けの意味なので、はCMが1or0であった場合のsalesの期待値を表していて、その下の式は、広告ありの購入額平均から広告無しの購入額平均を引くとつまり、498になりますよという意味になります。これが回帰分析が意味するところです。

平均因果効果とは難しく聞こえますが、要は平均値ということです。試しにそれぞれの群の平均値を求めてみると

となって、回帰分析の結果とも一致します。二群のt検定を行っても、今回の回帰分析の結果と一致しますので、興味のある方はやってみてください。

交絡因子を調整するということ

ここで、グラフを男女別に書いてみましょう。

よく見てみると広告は女性に多く配信される傾向にあるようです。女性の方が購入しやすい商品であった場合は、性別が交絡因子となり、二群の間で交換可能性が成り立っていません。つまり、二群の平均値は単純には比較できないということです。性別という交絡因子を除くにはどうしたらよいでしょうか？

すぐに思いつくのは、男性と女性の層に分けて回帰分析を行う方法です。男性、女性に分けて回帰分析を行ってみると、以下のようになります。

広告の平均因果効果である傾きは男性で-148、女性では438となりました。

これは何を求めたかというと、性別(sex : 0 男 1 女)という条件のもとでの広告効果という解釈になります。式で表すと

= -148

= 438

となります。ここから言えることは、女性限定という条件のもとでは広告の有無で購買額に違いがみられましたが、男性には効果がないことがわかります。”いつまでも美しく、今なら無料サンプル”といったような広告だったのでしょうか。（ちなみに男性の傾きの95%信頼区間は-338.0871~43.55145であり、有意ではありません。）

ついでに全ての層での平均値も求めておきましょう。 単純にそれぞれのグループで平均値を出すだけです。

• 広告なしかつ男性：5063
• 広告ありかつ男性：4916
• 広告なしかつ女性：6068
• 広告ありかつ女性：6506

重回帰分析は万能ではない

さて、層別化する以外の手法として、有名なものに重回帰分析があります。式としては、

or

となり、これを求めるのが重回帰分析です。実際に求めて、グラフに書いてみると

となり、男性の場合の広告効果は、

= 94

女性の場合の広告効果は

= 94

となりますね。先ほどと比べてどうでしょうか？本来広告効果がないはずの男性においても広告効果があるという結果になり、女性においては本来の広告効果が弱まっています。

ここが、重回帰分析を使う際のもっとも注意すべき前提です。重回帰分析は、要因の効果は全て独立に結果に作用するということを前提とします。

重回帰分析では、広告の効果は男女ともに一定であるという前提があるため、先ほどの例のように、広告の効果が男女で異なる場合は、単純には重回帰分析が使えないわけです。

一応回帰式から求められるそれぞれのグループでの値も出しておきましょう。

• 広告なしかつ男性：CM=0,sex=0 4977
• 広告ありかつ男性：CM=0,sex=0 5072
• 広告なしかつ女性：CM=0,sex=0 6296
• 広告ありかつ女性：CM=0,sex=0 6391

やはり上述した制約のため、層別化して求めた本来の平均値とは異なりますね。

交互作用の威力

さて、重回帰分析が使えないという場合はお手上げでしょうか？

いやいや、そんなことはありません。実は交互作用なるものを使って、対処することができるのです。交互作用とは、変数単体の影響だけでなく、変数同士が掛け合わさった場合にも影響が生じることを言います。今回の場合であれば、広告ありと女性という変数がかけ合わさって初めて、広告効果が出てきましたね。これが交互作用です。

この交互作用を重回帰分析に導入するためには、広告と性別の変数を掛け合わせた変数を式に加えます。

この式のが交互作用の大きさを表しています。これを求めてみると

となります。交互作用であるは585となっていて、女性にのみ585の広告効果が認められることが分かります。これは交互作用を含めないモデルの係数より、女性のみで回帰分析を行った結果に近いですね。

ちなみに交互作用がある場合は、CMの係数である-148やsexの係数である1004の解釈は難しくなります。CMを出せば購買額が148下がるという解釈や、女性の購買額は男性より1004高いなどのは解釈できないので、注意しましょう。

このモデルでも、それぞれのグループでの値を出してみます。

• 広告なしかつ男性：CM=0,sex=0 5063
• 広告ありかつ男性：CM=0,sex=0 4915
• 広告なしかつ女性：CM=0,sex=0 6067
• 広告ありかつ女性：CM=0,sex=0 6505

層別化した場合の平均値と同じ値になりました！ これで、要因の影響が条件間で異なることを考慮した因果推論を行うことができますね。

連続量の扱い方

さて、今までは全ての変数が0ro1の2値の場合を考えてきました。変数が連続量の場合はどうでしょうか？先ほどと少しデータを変えて、広告の有無ではなく広告費、性別ではなく年齢としたデータを考えてみます。

データはこんな感じになります。

広告費が高いほど、購買額は大きくなる傾向にありそうですが、そもそも年齢が高いほど広告費が高いようにも見えます。年齢が交絡バイアスを引き起こしている可能性があります。

ここで、先ほどのように層別化をしようとしてもできませんね。なぜなら、変数が連続量であり、層別化するとなると広告費が4000かつ年齢が25、広告費が4001かつ年齢が25、広告費が4002かつ年齢が26、広告費が4003かつ年齢が26....

など層別化する組み合わせの数が膨大になってしまい、比較ができません。

そこで、

という回帰式から、広告費と購買額の因果関係を推論してみましょう。これを求めると

となって、年齢という交絡を除いても、広告費が1増えると購買額が1増えるという結果になりました。

連続量の解釈

さて、このモデルが意味するところを考えてみます。先ほどの2値データでは、広告ありなしと性別男女の全ての組み合わせデータが存在していましたが、今回の場合は広告費と年齢の全ての組み合わせは存在していません。しかし、モデルからは広告費と年齢がどんな組み合わせであっても、購買額を推定できます。これは、連続量のモデルが、勝手に値を推定しているからです。

例えば、広告費5000かつ年齢25歳の購買額は、広告費4999かつ年齢25歳の購買額と広告費50001かつ年齢25歳の購買額の間ぐらいだろうとなります。

ここが連続量モデルの重要なポイントで、どの年齢層においても広告費と購買額の効果は一定かつ線形的であるという前提があります。このモデルでは、年齢41歳の人に、広告費6000を使うと莫大な売り上げになるという特異的な現象があったとしても、それは無視されてしまいます。これが連続量の変数を使った回帰分析で注意すべき所です。

ちなみに交互作用を使うこともできて、年齢と広告費の相乗効果みたいな現象であれば、表すことができます。

結局回帰分析の意味って？

ごちゃごちゃと書いてきましたが、結局回帰分析を使うということは何を意味しているかというと、条件付きの平均値を求めていることと同等になります。

性別同じであった場合の広告の有無による平均購買額の違い、年齢が同じだった場合の広告費による平均購買額の違いなどなど

条件を統制した上で、要因が結果に及ぼす影響力の平均を求めているのであって、要因が結果に及ぼす全体的な影響力を表してはいないということです。

これを間違えると、結果の解釈に大きな違いが出てきます。特に交互作用のような 影響がない場合は上記の二つは同じ結果になりますが、交互作用がある場合は異なる結果となります。

性別と広告有無モデルで、交互作用を考慮しなかった場合は、広告の効果が過小評価されていましたね。

これはモデルがパラメトリックなモデルか、ノンパラメトリックなモデルかに依存するのですが、詳しくは以下の記事や参考文献（ハーバード）を参考にされると良いかもしれません^^;

回帰分析を使った因果推論の仮定：パラメトリックモデルを使うということ - Unboundedly

まとめ

今回は、因果推論の基本テクニックである回帰分析について書いてみました。何気なく使われる回帰分析も、どんな仮定を置いていて、係数が何を意味するのかはしっかりと理解しておきたいところです。

回帰分析以外にも、まだまだ因果推論の分析手法はたくさんあるので、後々紹介してきます。

※本記事は筆者が個人的に学んだことをまとめた記事なります。数学の記法や詳細な理論、筆者の勘違い等で誤りがあった際はご指摘頂けると幸いです。

www.medi-08-data-06.work

www.medi-08-data-06.work

参考

理論的な背景までしっかりとまとまっていますので、詳しい理論が知りたい方にはオススメです。

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

ハバード大学が無料で公開している因果推論に関するテキストです。RやPythonのコードもあり、これが無料で読めるとは本当にすごいです。本記事の内容もこちらを参考にしています。

https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/

因果推論に関して、おそらく日本で一分かりやすく、しかも丁寧に解説されているブログです。大変参考にさせて頂いてます。

データから因果関係をどう導く？：統計的因果推論の基本、「反事実モデル」をゼロから - Unboundedly

Rのコード

本記事内で使ったデータの作成や回帰分析などのＲコードです。そのうちに実践編としてもまとめますが、ご参考までに

library(tidyverse)
#データセットの作成と可視化
set.seed(123)
CM <- sample(c(0,1),replace = T,200)
sex <- ifelse(CM+rnorm(200,0,1)<0.7,0,1)
sales <- 5000+1000*sex+500*CM*sex+rnorm(200,0,500)
dat <- data_frame(CM=CM,
                  sex=as.factor(sex),
                  sales=sales)

dat %>% 
  ggplot(aes(x=CM,y=sales))+
  geom_smooth(method = "lm",se=F)+
  geom_jitter(width = 0.01)+
  scale_color_manual(values = c("blue","red"))+
  annotate("text",x=0.5,y=7000,label="sales=5336+498CM",size=6)

#男女別の回帰分析
dat %>%
  filter(sex==0) %>% 
  lm(data = .,sales~CM) -> fit0

dat %>%
  filter(sex==1) %>% 
  lm(data = .,sales~CM) ->fit1

dat %>% 
  ggplot(aes(x=CM,y=sales,col=sex))+
  geom_jitter(width = 0.01)+
  scale_color_manual(values = c("blue","red"))+
  geom_smooth(method = "lm",se=F)+
  annotate("text",x=0.5,y=4000,
           label=str_c("sales=",floor(fit0$coefficients[1]),"",floor(fit0$coefficients[2]),"CM"),size=6,col="blue")+
  annotate("text",x=0.5,y=7000,
           label=str_c("sales=",floor(fit1$coefficients[1]),"+",floor(fit1$coefficients[2]),"CM"),size=6,col="red")

#重回帰分析
dat %>% 
  lm(data=.,sales~CM+sex)->fit3

#可視化
dat %>% 
  ggplot(aes(x=CM,y=sales,col=sex))+
  geom_jitter(width = 0.01)+
  scale_color_manual(values = c("blue","red"))+
  geom_abline(slope = fit3$coefficients[2],intercept = fit3$coefficients[1],col="blue")+
  geom_abline(slope = fit3$coefficients[2],intercept = fit3$coefficients[1]+fit3$coefficients[3],col="red")+
  annotate("text",x=0.5,y=7500,        label=str_c("sales=",floor(fit3$coefficients[1]),"+",floor(fit3$coefficients[2]),"CM+",floor(fit3$coefficients[3]),"sex"),size=5,col="black")+
  annotate("text",x=0.5,y=6500,label="sex=1",size=5,col="red")+
   annotate("text",x=0.5,y=5200,label="sex=0",size=5,col="blue")


#それぞれの平均値
new <- data_frame(CM=c(0,1,0,1),sex=as.factor(c(0,0,1,1)))
floor(predict(fit3, newdata = new))

#交互作用のモデル
dat %>% 
  lm(data=.,sales~CM+sex+sex*CM) ->fit4

#平均値
floor(predict(fit4,newdata = new))


#連続量データの作成
set.seed(123)
CM <- rnorm(200,5000,500)
age_mean <- ifelse(CM<5000,30,50)
age <- rnorm(200,age_mean,10)
sales <- 5000+10*age+CM+rnorm(200,0,500)
dat <- data_frame(CM=CM,
                  age = age,
                  sales = sales)
dat %>% 
  ggplot(aes(x=CM,y=sales,col=age))+
  geom_jitter(width = 0.01)+
  scale_color_continuous(low="#bbbbff",high = "#009900")

#回帰分析
dat %>% 
  lm(data=.,sales~CM+age) ->fit5