他の変数の影響を「調整」する道具としての回帰分析

回帰分析は「予測のための手法」として用いることもできますが、他の変数の影響を「調整」して、見たい効果を見るための道具としても頻繁に使用します。この記事では、単純な集計がはっきり逆向きの答えを返す場面から出発して、調整する道具としての回帰分析の考え方を一つのシミュレーションデータで通して見ていきます。

単純比較が誤るとき

「新機能を使っているユーザーは、使っていないユーザーより、翌月の利用日数が平均3日少ない」。あるサブスク型アプリで、こんな集計結果が出たとします。ユーザーの定着を助けるために作った機能なのに、この数字だけを見ると逆効果だったように見えます。

しかし実際には「新機能は利用日数を2日増やす」というのが正しい効果でした。全員に等しくプラスの効果があるのに、単純な平均の比較でははっきりマイナスに見える。このようなことが起こり得ます。

状況をもう少し具体的にします。ユーザーは5,000人。月に6日ほど使うライト層と、18日ほど使うヘビー層が半々います。ここに定着支援のための新機能をリリースし、アプリ内の案内をライト層に重点的に出しました。もともと使い込んでくれているヘビー層より、離脱が心配なライト層にこそ使ってほしい機能だからです。結果として、ライト層の約6割、ヘビー層の約2割が新機能を使いました。

新機能の真の効果は「翌月の利用日数を2日増やす」と設定します。コードにするとこうです。

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
n = 5000

# 利用頻度セグメント: ライト層とヘビー層が半々
segment = rng.choice(["light", "heavy"], size=n, p=[0.5, 0.5])

# 新機能の利用: ライト層に重点的に案内したため、ライト層ほど使っている
p_use = np.where(segment == "light", 0.6, 0.2)
feature = rng.binomial(1, p_use)

# 翌月の利用日数: ライト層は平均6日、ヘビー層は平均18日。新機能の効果は +2日
base = np.where(segment == "light", 6.0, 18.0)
days = (base + 2.0 * feature + rng.normal(0, 2.5, n)).clip(0, 31)

df = pd.DataFrame({"segment": segment, "feature": feature, "days": days})

新機能を使ったグループと使っていないグループで、翌月の利用日数を単純に比べてみます。

df.groupby("feature")["days"].mean()
# feature
# 0    13.91
# 1    10.91

効果を+2日に設定したはずなのに使ったグループのほうが3日少ない結果が返ってきました。

種明かしは人数の内訳にあります。

人数	新機能なし	新機能あり
ライト層	1,002人	1,543人
ヘビー層	1,958人	497人

案内をライト層に寄せたので、「新機能あり」グループの4分の3はライト層です。逆に「新機能なし」グループの3分の2はヘビー層。つまり「新機能あり vs なし」の比較は、中身を見れば「ほぼライト層 vs ほぼヘビー層」の比較になっています。ライト層はもともと利用日数が少ないので、新機能の+2日はこの層の差に埋もれて、逆向きの結果が出てきます。

セグメントごとに分けて平均を見ると、結果は次のようになります。

翌月の利用日数	新機能なし	新機能あり	差
ライト層	6.08日	8.02日	+1.94日
ヘビー層	17.91日	19.90日	+1.98日

どちらの層でも、新機能を使ったグループのほうが2日ほど多い。このように、グループ分けしたデータのそれぞれで成り立つ関連の向きが、全体を合算すると逆転して見える現象は、シンプソンのパラドックスとして知られています。名前はパラドックスですが、種を明かせば矛盾はどこにもなく、いま見たとおり、比べている2つのグループの中身が偏っていることから起きる算術的な現象です。

いま起きたことを図にすると、こうなります。

利用頻度セグメントが新機能の利用と結果の両方に影響している構造

利用頻度セグメントは、新機能を使うかどうか(案内の偏り)と、翌月の利用日数(もともとの差)の両方に影響しています。このように、原因と結果の両方に影響して比較を歪める変数を交絡変数と呼びます。単純比較が誤るのは、交絡変数の影響が、見たい効果に混ざり込むからです。「相関があっても因果があるとは限らない」と言われるときの代表格、第三の変数が作る擬似相関も、この構造から生まれます。

層別解析

先ほどセグメントごとに分けて計算したように、層の中で比較する方法を層別解析と呼びます。層の中では全員が同じセグメントなので、セグメントの影響は差し引かれ、新機能の効果だけが残ります。「条件が同じ人同士を比べる」という、比較の原則に忠実なやり方です。

混ぜて比べた場合と、層の中で比べた場合を並べると、こうなります。

単純比較では下を向いていた線が、層の中で比べると両方とも上を向く

左は5,000人を新機能の利用だけで2つに分けた比較で、線は下を向いています。右は同じ5,000人をセグメントで分けてから、それぞれの中で比べたものです。ライト層とヘビー層は日数の水準がそもそも12日ほど離れていて、左の比較はこの2つを混ぜたものでした。層の中に入れば比べる相手の水準が揃うので、線は両方とも上を向きます。

層ごとに出た2つの効果は、1つの数字にまとめておくこともできます。全体としての効果を見るなら、人数の多い層の値を重く扱います。ライト層は2,545人(1,002人 + 1,543人)で+1.94日、ヘビー層は2,455人(1,958人 + 497人)で+1.98日でした。人数の割合はライト層が51%、ヘビー層が49%なので、この割合で混ぜます。

1.94日 × 0.51 + 1.98日 × 0.49 ≒ +1.96日

今回は層の人数がほぼ半々なので、2つを単純に平均したのとほとんど変わりません。真の効果である2日に近い値が出ました。交絡の問題は、これで解決したように見えます。

層別解析の限界

ただ、結果に影響していそうな変数は、実際にはセグメントだけではありません。料金プラン、流入チャネル、利用歴。層別解析の理屈をそのまま通すなら、これらすべてが同じ人同士を比べることになります。

df["tenure"] = rng.integers(1, 25, n)  # 利用歴(月)
df["plan"] = rng.choice(["basic", "standard", "premium"], n, p=[0.5, 0.3, 0.2])
df["channel"] = rng.choice(["organic", "ad", "referral"], n, p=[0.4, 0.4, 0.2])

層の数は掛け算で増えます。セグメント2 × プラン3 × チャネル3 = 18層までは、まだ1層に280人ほどいます。ところが利用歴(24通り)を足した瞬間に432層。5,000人を割ると1層あたり12人弱で、その中をさらに新機能あり・なしに分ければ、比べる相手は片側6人ほどです。実際、36の層では片側が1人もおらず、比較そのものが成立していません。

しかも利用歴は、月単位に丸めるという妥協をまだしています。日単位にすれば層はさらに増えますし、累計課金額のような連続値なら「同じ値の人同士」はそもそもほぼ存在しません。

層別の一般化としての重回帰分析

この限界を避けながら、層別と同じ発想を実現するのが重回帰分析です。説明変数を複数持つ回帰分析で、いまの例では式はこうなります。

$\text{利用日数} = \beta_0 + \beta_1 \cdot \text{新機能ダミー} + \beta_2 \cdot \text{ライト層ダミー}$

式に出てくる「新機能ダミー」「ライト層ダミー」は、当てはまれば1、当てはまらなければ0を取る変数です。こうした0か1かの変数をダミー変数と呼びます。回帰分析が受け取れるのは数値だけなので、セグメントのようなカテゴリは、この形に置き換えて式に入れます。なお、ヘビー層のダミーは置きません。ライト層ダミーが0のときがヘビー層なので、1つで足ります。

知りたいのは係数 $\beta_1$ 、つまり新機能ダミーが0から1に変わったときの利用日数の変化です。式の中にライト層ダミーの項が別にあるおかげで、セグメントによる日数の差はそちらの $\beta_2$ が引き受け、 $\beta_1$ には混ざりません。こうして得られる係数の解釈は、「他の変数の値を固定したうえで、対象の変数を1単位動かしたときに、目的変数がどれだけ変動するか」です。いまの例で言えば、 $\beta_1$ は「セグメントが同じユーザー同士で比べたときの、新機能の効果」と読めます。なお、こうして他の変数の影響を差し引くことを、この記事では「調整」と呼んでいますが、心理学では「統制」と呼ぶことが多く、指しているものは同じです。

statsmodels で当てはめてみます。式は 目的変数 ~ 説明変数 の形で書き、セグメントのようなカテゴリの列は C(...) で囲みます。囲んでおくと、さきほどのダミー変数への置き換えを statsmodels がやってくれます。

import statsmodels.formula.api as smf

model = smf.ols("days ~ feature", df).fit()
model.params["feature"]
# -3.00 (単純比較の差と同じ)

model = smf.ols("days ~ feature + C(segment)", df).fit()
model.params["feature"]
# 1.96

説明変数が新機能だけの回帰は、単純比較の差をそのまま返します。セグメントを加えると、係数は+1.96日。層別解析で手計算した重み付き平均の+1.96日と一致します。segment を説明変数に足した回帰は、「層別して、層ごとの差を重み付きで平均する」のとほぼ同じ計算を一度にやってくれているわけです(厳密には回帰の重み付けは人数比と少し異なりますが、効果が層でほぼ同じ今回は、その違いは結果に現れません)。

調整したい変数を増やしたいときには、回帰では式に項を足すだけで済みます。

model = smf.ols(
    "days ~ feature + C(segment) + C(plan) + C(channel) + tenure", df
).fit()
model.params["feature"]
# 1.96

※追加の3変数は結果に影響しないように作ったので、係数は変わりません。

このように回帰でうまくいくのは、層ごとに切り分けて計算しないからです。層別解析は、432層のそれぞれで新機能あり・なしの平均を出すので、864個の平均を数人ずつのデータから求めることになります。一方の回帰が求めるのは、4つの変数すべてを入れても8個の数字だけです。「新機能で何日変わるか」「ライト層だと何日少ないか」「利用歴が1ヶ月増えると何日変わるか」といった数字を、どれも5,000人全員のデータから計算します。1人ひとりの利用日数は、この8個を足し合わせて組み立てるので、細かい組み合わせに人がいなくても計算は止まりません。

ただしこれにも仮定が置かれています。8個で足りるのは、「各変数の影響は足し算で効く」「利用歴は1ヶ月あたり同じだけ効く」という形を決め打ちしているからです。仮定した形が現実と食い違っていれば、その分だけ答えは歪みます。たとえば新機能の効果がライト層とヘビー層で大きく違う場合、上の式はその違いを表現できません。層別が「仮定は少ないがデータを大量に食う」方法だとすれば、回帰は「形の仮定と引き換えに、少ないデータで調整をやりくりする」方法です。

効果が層によって違うとき

新機能の効果がライト層では大きく、ヘビー層ではほとんどない、ということは十分にありえます。定着に困っている層にこそ効く機能なら、むしろそちらが自然かもしれません。ところが先ほどの式では、新機能ダミーの項は1つしかなく、そこに入る係数も1つだけです。「効果は層によらず同じ」と決め打ちした形になっています。

この決め打ちを緩めるには、2つのダミー変数を掛け合わせた項を足します。

$\text{利用日数} = \beta_0 + \beta_1 \cdot \text{新機能ダミー} + \beta_2 \cdot \text{ライト層ダミー} + \beta_3 \cdot (\text{新機能ダミー} \times \text{ライト層ダミー})$

最後の項が交互作用項です。掛け算なので、新機能を使っているライト層のときだけ1になり、それ以外は0になります。

読み方は、ダミーに値を入れて式を書き下すとはっきりします。ヘビー層ならライト層ダミーは0、ライト層なら1です。まずヘビー層から。

$\text{新機能なし} = \beta_0 \qquad \text{新機能あり} = \beta_0 + \beta_1$

ライト層ダミーが0なので、 $\beta_2$ の項も交互作用項も消えます。新機能ありとなしの差は $\beta_1$ です。次にライト層。

$\text{新機能なし} = \beta_0 + \beta_2 \qquad \text{新機能あり} = \beta_0 + \beta_1 + \beta_2 + \beta_3$

こちらは交互作用項も残るので、差は $\beta_1 + \beta_3$ になります。4つを並べると、こうです。

予測される利用日数	新機能なし	新機能あり	新機能の効果
ヘビー層	$\beta_0$	$\beta_0 + \beta_1$	$\beta_1$
ライト層	$\beta_0 + \beta_2$	$\beta_0 + \beta_1 + \beta_2 + \beta_3$	$\beta_1 + \beta_3$

右端の2つをさらに引き算すると $(\beta_1 + \beta_3) - \beta_1 = \beta_3$ で、 $\beta_3$ は「ライト層のほうが効果がどれだけ大きいか」という層間の差そのものです。交互作用項のない前の式は、この $\beta_3$ を0に固定していた、とも読めます。

まず、ここまで使ってきたデータに当てはめてみます。効果を両層とも+2日に設定したデータなので、層間の差はないはずです。

model = smf.ols("days ~ feature * C(segment)", df).fit()

model.params["feature"]                      # 1.98
model.params["feature:C(segment)[T.light]"]  # -0.04

conf_int = model.conf_int()  # 各係数の95%信頼区間
conf_int.loc["feature:C(segment)[T.light]"]  # -0.36 〜 0.27

formula の * は、2つの変数それぞれの項と、その積の項をまとめて置く書き方です。出てくる係数の名前もそれに対応していて、feature:C(segment)[T.light] が交互作用項にあたります。: が掛け合わせ、[T.light] が「ライト層なら1」のダミーなので、名前としては「新機能ダミー × ライト層ダミー」を指しています。

交互作用の係数は-0.04で、95%信頼区間も0をまたいでいます。層によって効果が違うとは言えない、という読み方になります。効果が層によらず同じという仮定が妥当かどうかは、こうして項を足せば確かめられることもあります。

次に、効果が層で違うデータを作って、同じことをやってみます。ライト層では+4日、ヘビー層では効果なしという設定です。

rng2 = np.random.default_rng(123)
segment2 = rng2.choice(["light", "heavy"], size=n, p=[0.5, 0.5])
feature2 = rng2.binomial(1, np.where(segment2 == "light", 0.6, 0.2))
base2 = np.where(segment2 == "light", 6.0, 18.0)
effect2 = np.where(segment2 == "light", 4.0, 0.0)  # 効果が層で違う
days2 = (base2 + effect2 * feature2 + rng2.normal(0, 2.5, n)).clip(0, 31)

df2 = pd.DataFrame({"segment": segment2, "feature": feature2, "days": days2})

層別解析で層ごとの差を見ると、設定した値がそのまま出てきます。

翌月の利用日数	新機能なし	新機能あり	差
ライト層	6.00日	10.00日	+4.00日
ヘビー層	18.00日	18.12日	+0.12日

このデータに、まず交互作用項のない式(前の節と同じ重回帰)を当てはめてみます。

model = smf.ols("days ~ feature + C(segment)", df2).fit()
model.params["feature"]
# 2.55

+2.55日という、ライト層の+4.00日でもヘビー層の+0.12日でもない数字が出てきます。この式は効果を1つだと決めているので、2つの層の効果を1つに均してしまいます。しかもその均し方は人数比での平均とも一致しないため、+2.55日が何の平均なのかを言葉で説明するのも難しくなります。どちらの層のユーザーを見ても、この+2.55日に当てはまる人はいません。

交互作用項を入れると、2つの効果が分かれて出てきます。

model = smf.ols("days ~ feature * C(segment)", df2).fit()

model.params["feature"]                      # 0.12  ヘビー層での効果
model.params["feature:C(segment)[T.light]"]  # 3.88  ライト層はその分大きい
# ライト層での効果 = 0.12 + 3.88 = 4.00

層別解析で出した+0.12日と+4.00日に、そのまま一致します。これは偶然ではありません。セグメントと新機能のすべての組み合わせに項を置いた式は、層ごとに別々の平均を計算するのと同じことをしています。回帰は層別解析の一般化だと書きましたが、その関係が一番はっきり見えるのがこの形です。層別解析は、回帰の特殊な場合として中に含まれています。

それなら考えつく限りの変数に交互作用項を置けばよさそうですが、そうはいきません。掛け合わせる変数を増やすほど項は増え、1つ1つの項の推定に使えるデータは薄くなります。すべての変数のすべての組み合わせに交互作用を置いた式は層別解析そのものですから、そこまで行けば、層別解析と同様の限界に突き当たってしまいます。効果が層によって違いそうだと考える理由のある変数にだけ置く、という使い方になります。

調整してはいけない変数

実際にやってみるときに気になるのは、関係しそうな変数をすべて説明変数に入れてよいのか、という点です。これについては、入れてよい変数といけない変数があります。目安を一つだけ挙げると、調整してよいのは、原因(ここでは新機能の利用)より時間的に前に決まっている変数です。

原因より後に決まる変数を入れると、逆に効果が見えなくなることがあります。例えば、新機能が「通知をオンにするユーザーを増やし、その通知によって利用日数が増える」という経路で効いているとしましょう。このとき、通知の設定まで説明変数に入れてしまうと、通知が同じ状態のユーザー同士を比べることになります。新機能がどれだけ通知をオンにさせていても、その分の影響は通知の項が引き受けてしまい、新機能の係数には残りません。結果として、効果は実際より小さく見えます。仮に、通知を経由しない別の経路でも新機能が効いていたとしても、係数に残るのはその直接の分だけで、知りたかった効果の全体ではなくなります。このような、効果の通り道になっている変数は媒介変数と呼ばれ、交絡変数と違って、調整してはいけない変数です。

2つの構造を並べると、違いがはっきりします。矢印の向きに注目すると、交絡変数は原因に矢印を「向けて」いて、媒介変数は原因から矢印を「受けて」います。

交絡変数は調整し、媒介変数は調整しない

加えて留意が必要なのは、調整できるのは測定できている変数だけ、ということです。今回はシミュレーションなので、交絡がセグメントだけだと知っていましたが、実際には測っていない交絡が残っている可能性を消しきれません。ランダムに割り当てるA/Bテストが実施できるなら、それが一番素直な解決です。回帰による調整が活きるのは、実験ができない・しそびれた場面で、観察データから効果に迫りたいときです。

どの変数なら調整してよいかについて、この記事では「原因より前か後か」という目安だけを紹介しましたが、バックドア基準という厳密な判定基準が知られています(宮川 (2004) が詳しいです)。また、調整のための道具も回帰だけではなく、条件の近いユーザー同士をペアにして比べるマッチングや、たくさんの変数を「新機能を使う確率」という1つの数値に要約してしまう傾向スコアといった発展があります。単純比較では同じ条件の人同士を比べられていない、という今回の問題意識さえ押さえておけば、これらは同じ目的のための別の道具として位置づけられます。

予測の道具、調整の道具

回帰分析の式は、予測に使うときも調整に使うときも同じです。違うのは見る場所で、予測では出力される予測値の精度に関心があり、調整では特定の係数1つの意味に関心があります。関心が違えば変数の選び方も変わります。予測なら精度が上がる変数は何でも歓迎ですが、調整では、いま見たように「原因より前か後か」を考えて選ぶ必要があります。

また、この調整の考え方は、今回扱った重回帰分析に限ったものではなく、回帰分析一般で使えます。たとえば目的変数が利用日数のような数値ではなく、「翌月に解約したかどうか」のような2値なら、モデルはロジスティック回帰に変わります。それでも、交絡変数を説明変数に足して調整するという発想も、「他の変数を固定したうえで」という係数の読み方も、調整してよい変数といけない変数の区別も、そのまま持ち運べます。変わるのは、係数の解釈がオッズ比を通じたものになるといった、モデルごとの読み替えの部分だけです。

「新機能を使ったユーザーのほうが利用日数が3日少ない」という最初の数字が映していたのは、効果ではなく、比べている2つのグループの中身の違いでした。その揃っていない中身を差し引いて、見たい効果だけを取り出すのが、調整する道具としての回帰分析でした。

参考

宮川雅巳 (2004). 『統計的因果推論: 回帰分析の新しい枠組み』朝倉書店.
高橋将宜 (2022). 『統計的因果推論の理論と実装: 潜在的結果変数と欠測データ』共立出版.
星野匡郎, 田中久稔 (2016). 『Rによる実証分析: 回帰分析から因果分析へ』オーム社.