第2章 介入効果を測るための回帰分析：Huluの推薦ログから効果を取り出す

第1歩：バイアスの原因が分かっているなら、それを揃えて比べればよい

第1章で見た通り、RCTでないデータの単純な平均の差には、効果とセレクションバイアスが混ざっている。

セレクションバイアスは、「介入を受けた人と受けなかった人で、もともとの傾向が違う」ことから生まれる。そして多くの場合、その違いを生んでいる変数には見当がつく。上のHuluの例なら、訴求対象を決めるときに使った「アニメを視聴していたか」がそれである。

原因が分かっているなら、対処は素朴でよい。 その変数の値が同じ人同士で比べればよい 。アニメ視聴者の中だけで訴求あり・なしを比べ、非視聴者の中だけで訴求あり・なしを比べれば、「訴求群にアニメ視聴者が多い」という偏りは比較の中から消える。この方法を 層別比較 と呼ぶ。

小さな例で確かめる。8人のプロフィールについて、 $X$ をアニメ視聴者か、 $Z$ を『名探偵コナン』を訴求したか、 $Y$ を7日間のコナン視聴分数とする（この表は2.2.2でも再び使う）。

全員をまぜて比べると、訴求あり4人の平均は70分、なし4人の平均は30分で、差は40分である。一方、層ごとに比べると次になる。

層ごとに比べると、どちらの層でも差は20分である。単純比較の40分のうち20分は、訴求群にアニメ視聴者が多かったことによる見せかけの差であった。 同じ条件の人同士で比べることで、セレクションバイアスが消えた 。これが、この章全体でやりたいことの原型である。

第2歩：層別比較は、揃えたい条件が増えるとすぐに破綻する

ところが、現実には揃えたい条件は一つではない。訴求対象を決めるルールが「過去90日のアニメ視聴日数が3日以上、かつ契約1か月以上、かつキッズプロフィールでない」だったとしよう。すると、揃えるべき変数は「アニメ視聴日数（0〜90日の整数）」「契約月数」「キッズか」の三つになる。さらに、視聴傾向に関係しそうな「総視聴時間」や「利用デバイス」も揃えたくなる。

条件を厳密に揃えるということは、「視聴日数37日、契約14か月、キッズでない、総視聴時間12.3時間、iOS」のような細かい層を作り、その中で訴求あり・なしを比べることである。しかし、そのような層には人がほとんどいない。訴求ありとなしの両方がそろっている層はさらに少なく、多くの層では差そのものが計算できない。変数が連続値であれば、値がぴったり一致する人はそもそも存在しない。

つまり層別比較は、考え方は正しいが、 条件の組み合わせが増えると、比較に使えるデータが枯渇して実行できなくなる 。

第3歩：層ごとの平均に「形」を仮定すると、人がいない層でも比べられる

行き詰まりを抜ける発想は、層ごとにバラバラに平均を計算するのをやめて、 「条件 $X$ と介入 $Z$ が与えられたときの $Y$ の平均」に、簡単な式の形を仮定する ことである。

$E[Y\mid X,Z]$ は「 $X$ と $Z$ がこの値の人たちの $Y$ の平均」であり、層別比較で計算していた「層ごとの平均」そのものである。この式は、その層ごとの平均が次の二つの規則に従うと仮定している。

• $X$ が1増えると、 $Y$ の平均は $\beta_1$ だけ増える（ $X$ の影響は直線的である）。
• $Z$ が0から1になると、 $Y$ の平均は $\beta_Z$ だけ増える（介入の効果はどの層でも同じである）。

この仮定が、なぜ「人がいない層でも比べられる」ことにつながるのか。層別比較と回帰で、データから決めなければならない数の個数を比べると分かる。

層別比較では、「層の数 × 2」個の平均を、それぞれ別々に求める必要がある。たとえば $X$ が0から3までの4通りなら、 $X=0$ の訴求なし平均、 $X=0$ の訴求あり平均、 $X=1$ の訴求なし平均、というように8個の平均を求め、各層でその差を取る。8個の平均はお互いに無関係に計算されるので、8個のマスすべてに人がいなければならない。一つでも空のマスがあれば、その層の差は計算できない。

一方、回帰では、決めなければならない数は $\beta_0,\beta_1,\beta_Z$ の3個だけである。しかも、この3個は8個のマスに散らばった全員のデータを一緒に使って決める。 $\beta_1$ は「 $X$ が1増えると $Y$ の平均がいくら増えるか」という規則なので、どの層の人からも情報を得られる。 $\beta_Z$ は「訴求で $Y$ の平均がいくら増えるか」なので、訴求あり・なしの両方がいる層がどこか一つでもあれば決まる。3個の数が決まれば、式に $X$ と $Z$ を代入することで、 人がいないマスの平均も計算できる 。これが「隣の層の情報から式で補える」という意味である。

ワークド例：空のマスを式で埋める

$X$ を「過去90日のアニメ視聴日数」（0〜3日）とし、8人のデータを次のように持っているとする。話を単純にするため、誤差のないデータとする。

これを層別比較の表に並べ替えると、空のマスが二つある。視聴日数0日の人は誰も訴求されておらず、3日の人は全員が訴求されている。訴求ルールが「視聴日数が多い人ほど訴求する」というものであれば、このように両端の層が片側だけになるのは自然なことである。

層別比較では、 $X=1$ と $X=2$ の層からしか差が得られない。 $X=0$ と $X=3$ の4人は、比較に使えずに捨てることになる。

回帰では、8人全員を使って3個の数を決める。この表のデータは $Y=20+10X+20Z$ にぴったり乗っているので、最小2乗法は $\hat{\beta}_0=20$ 、 $\hat{\beta}_1=10$ 、 $\hat{\beta}_Z=20$ を返す（2.1.1で最小2乗法の計算を確認したあと、2.1.5の ols() にこの表を入れれば確かめられる）。

• $\hat{\beta}_1=10$ は、訴求なしの列が $X=0,1,2$ で20、30、40と10ずつ増え、訴求ありの列が $X=1,2,3$ で50、60、70と10ずつ増えていることから決まる。 どちらの列も「1日あたり10分」という同じ規則に従う という仮定があるから、二つの列の情報を合わせて一つの $\hat{\beta}_1$ を決められる。
• $\hat{\beta}_Z=20$ は、両群がいる $X=1$ と $X=2$ の層内の差20から決まる。

3個の数が決まれば、空のマスを式で埋められる。 $X=3$ の訴求なしの平均は $20+10\times3=50$ 、 $X=0$ の訴求ありの平均は $20+0+20=40$ と計算できる。

「 $X=3$ の訴求なし平均は50」という数は、 $X=3$ の訴求なしの人を一人も観測していないのに得られている。 $X=0,1,2$ の訴求なしの人たちが示す「1日あたり10分」という傾きを、 $X=3$ まで延長したからである。これが「隣の層の情報から補う」ということの正確な意味である。

同じことを、横軸に $X$ 、縦軸に $Y$ を取った散布図で見ると、回帰が何をしているかがはっきりする。訴求なしの人たちと訴求ありの人たちに、同じ傾きの直線を1本ずつ引き、人がいない $X$ の位置まで直線を延長しているのである。

この例では $X$ の値が4通りしかないので、層別比較でもなんとか差が求まった。しかし共変量が「視聴日数（0〜90）」「総視聴時間」「契約月数」の三つになると、マスの数は数万個に膨れ上がり、ほとんどのマスは空になる。それでも回帰で決めるべき数は $\beta_0$ と各共変量の係数と $\beta_Z$ の5個だけであり、5個を決めるのに必要なのは「全体として十分な人数」であって「すべてのマスに人がいること」ではない。 これが、条件が多くても比較を続けられる理由である 。

ただし、この「補う」操作は仮定に依存していることも忘れてはいけない。上の例で $X=3$ の訴求なし平均を50と置けたのは、「1日あたり10分という増え方が $X=3$ でも続く」と仮定したからである。もし実際には視聴日数3日の人の視聴分数が頭打ちになっていて45分だったなら、補った値はずれ、 $\hat{\beta}_Z$ にもずれが入る。だからこそ2.1.2では「共変量の関数の形は分析者が決めなければならない」という話が出てくる。また、層ごとの差が本当に同じかどうかは、2.4.5で改めて問い直す。

そして、この式のもとで「 $X$ が同じ人の中での、訴求ありとなしの平均の差」を計算すると、 $X$ に関する項は打ち消し合う。

つまり $\beta_Z$ は、 層別比較で求めていた「同じ条件の人同士の差」を、すべての層についてまとめて一つの数で表したもの である。8人の例にこの式を当てると $\beta_Z=20$ が得られ、層別比較の答えと一致する（計算は2.2.2で確認する）。層ごとの差が同じでない場合に $\beta_Z$ が何を表すかは、2.4.5で扱う。

この節の残りで確認すること

ここまでで、回帰分析が効果検証に登場する理由は説明できた。ただし、この議論は「 $E[Y\mid X,Z]$ に式の形を仮定すれば、そのパラメータを推定できる」という点に依存している。では、その式のパラメータは、手元のデータからどうやって求めるのか。そして、求めた値はなぜ「条件付き期待値」の推定と言えるのか。この二つを確認するのが、教科書が最初に置いている単回帰の説明である。したがって、以下の2.1.1は「線の引き方の説明」ではなく、 「回帰で引いた線は条件付き期待値の推定である」ことを確立するための準備 として読むとよい。それが確立できれば、2.1.3で「条件付き期待値の差 $=\beta_Z$ 」という結論に、そのままつながる。

ワークド例：4人のデータから傾きと切片を手計算する

抽象的な式を、手を動かして確かめる。Huluのあるプロフィール4人について、 $X$ を「過去4週間の視聴日数」、 $Y$ を「今週の総視聴分数」とする。

平均は $\bar X=5$ 、 $\bar Y=60$ である。単回帰の傾きは「 $X$ と $Y$ の共分散を $X$ の分散で割ったもの」であり、分母と分子の $1/N$ は打ち消し合うので、ずれの積の合計をずれの2乗の合計で割ればよい。

切片は、線が点 $(\bar X,\bar Y)$ を必ず通ることから求まる。

したがって、推定された線は $\hat Y=5+11X$ である。「過去の視聴日数が1日多いと、今週の視聴分数は平均して11分多い」と読める。予測値と誤差を確かめると、Aは $5+22=27$ で誤差は $+3$ 、Bは $49$ で $+1$ 、Cは $71$ で $-11$ 、Dは $93$ で $+7$ となり、誤差の合計は0になる。これは最小2乗法の解が必ず満たす性質である。

この計算を図にすると次のようになる。4つの点に対して、赤い縦線（残差）の2乗の合計がもっとも小さくなる位置に直線が置かれている。

回帰は「近似の線」であるだけでなく「条件付き期待値」の推定である

ここが2.1.1の難所である。回帰分析は、手元のデータに線を当てるだけの操作ではない。手元のデータで得た $\hat{\beta}$ は、 母集団で回帰を行ったときに得られるパラメータ $\beta$ の推定値 である。そして母集団では、 $Y$ と $X$ の関係を次のように書くことができる。

$E[Y\mid X]$ は「 $X$ がある値のときの $Y$ の平均」であり、 条件付き期待値 と呼ぶ。この式が言っているのは、次の二つのことである。

• どんな $Y$ も、「 $X$ で決まる平均の部分」と「そこからのずれ $u$ 」に分解できる。
• 条件付き期待値 $E[Y\mid X]$ は、 $Y$ を予測する関数の中で、2乗誤差がもっとも小さくなる関数である。

したがって、最小2乗法で線を引くという操作は、条件付き期待値 $E[Y\mid X]$ を線形の形で近似する操作に一致する。上のワークド例で言えば、「 $\hat Y=5+11X$ 」は「過去の視聴日数が $X$ 日の人たちの、今週の平均視聴分数はおよそ $5+11X$ 分である」という条件付き期待値の推定である。この見方が、次節で回帰係数を「効果」として解釈するための土台になる。

共変量の「関数の形」は分析者が決める

重回帰では、 $X$ の2乗のような変数を追加できる。

なぜこのようなことが必要になるのか。たとえば、教育方法の年収への効果を推定するときに、共変量として年齢を入れる場合を考える。年収は経験や年齢によって変わるが、その関係は直線ではない。年齢が上がるほど1年分の経験が年収に与える影響は小さくなるとすれば、その関係は2次関数で表される。このような場合、 $X$ だけでなく $X^2$ も加える必要がある。

Huluでも同じことが起きる。「過去の総視聴時間」と「今週の視聴確率」の関係は直線ではなく、視聴時間が長い人ほど伸びが鈍る（天井に近づく）ことが多い。したがって、視聴時間をそのまま入れるのではなく、対数を取る、2乗項を入れる、あるいは区分（0時間、1〜5時間、5時間以上など）のダミー変数にするといった工夫が必要になる。どの形が適切かは、データを眺めたうえで分析者が事前知識から決める。

一方で、 介入変数 $Z$ は0か1しか取らないため、関数の形を考える必要がなく、常に線形のまま扱えばよい 。0と1の2点を通る関係は必ず直線で表せるからである。この点は2.4.2でもう一度出てくる。

ワークド例：具体的な数値で差を確認する

教科書の数値をそのまま使う。手元のデータで回帰を行い、 $\hat{\beta}_0=0.5$ 、 $\hat{\beta}_1=0.2$ 、 $\hat{\beta}_2=-0.02$ 、 $\hat{\beta}_Z=0.4$ が得られたとする。共変量が $X=5$ の人について計算すると、次になる。

差は $0.4$ である。次に $X=10$ の人では、 $0.5+2.0-2.0+0.4=0.9$ と $0.5+2.0-2.0=0.5$ で、やはり差は $0.4$ になる。共変量の値がいくらであっても、介入ありとなしの差は $\hat{\beta}_Z=0.4$ で一定である。これが「介入グループでは $Y$ の平均が非介入グループより0.4高い」という解釈の意味である。

$X$ を横軸に取って二つの条件付き期待値を描くと、共変量との関係が曲線であっても、介入ありの曲線は介入なしの曲線を縦にちょうど $\hat{\beta}_Z$ だけ持ち上げたものになっていることが分かる。

なお、この「差が $X$ に依存しない」という性質は、モデルに $Z$ と $X$ の交互作用項を入れていないことの結果である。効果が人によって違う場合の扱いは、2.4.5と追加手法の節で扱う。

最小2乗法を自分で書く

計画行列（説明変数を横に並べた行列）を $\mathbf{X}$ 、目的変数のベクトルを $\mathbf{y}$ とすると、最小2乗法の解は次の閉じた式で求まる。

標準誤差は、誤差の分散の推定値 $\hat\sigma^2=\text{RSS}/(N-k)$ （RSSは残差の2乗和、 $k$ はパラメータ数）を使って、次の行列の対角成分の平方根として得られる。

これはRの lm() の既定と同じ計算であり、誤差の分散がどの人でも同じ（等分散）という仮定に基づいている。この仮定が怪しいときの対処は追加手法の節で述べる。

from dataclasses import dataclass
import numpy as np
import polars as pl
from scipy import stats

@dataclass
class OLSResult:
    terms: list[str]
    coef: np.ndarray
    se: np.ndarray
    residuals: np.ndarray
    n: int
    k: int

    @property
    def tidy(self) -> pl.DataFrame:
        # R の broom::tidy() に相当する表を返す
        dof = self.n - self.k
        t_value = np.where(self.se > 0, self.coef / self.se, np.inf)
        p_value = 2 * stats.t.sf(np.abs(t_value), df=dof)
        return pl.DataFrame({
            "term": self.terms, "estimate": self.coef, "std_error": self.se,
            "statistic": t_value, "p_value": p_value,
        })

    def coef_of(self, term: str) -> float:
        return float(self.coef[self.terms.index(term)])


def design_matrix(frame: pl.DataFrame, columns: list[str]) -> tuple[np.ndarray, list[str]]:
    """切片つきの計画行列を作る。文字列列は最初の水準を基準にしたダミー変数へ展開する。"""
    parts = [np.ones((frame.height, 1))]
    names = ["(Intercept)"]
    for col in columns:
        series = frame.get_column(col)
        if series.dtype == pl.String:
            levels = sorted(series.unique().to_list())
            for level in levels[1:]:
                parts.append((series == level).cast(pl.Float64).to_numpy().reshape(-1, 1))
                names.append(f"{col}{level}")
        else:
            parts.append(series.cast(pl.Float64).to_numpy().reshape(-1, 1))
            names.append(col)
    return np.hstack(parts), names


def ols(frame: pl.DataFrame, outcome: str, predictors: list[str]) -> OLSResult:
    used = frame.select([outcome, *predictors]).drop_nulls()   # R の na.omit と同じ
    x, names = design_matrix(used, predictors)
    y = used.get_column(outcome).cast(pl.Float64).to_numpy()
    xtx_inv = np.linalg.pinv(x.T @ x)
    beta = xtx_inv @ x.T @ y
    residuals = y - x @ beta
    n, k = x.shape
    sigma2 = max(float(residuals @ residuals) / (n - k), 0.0)
    se = np.sqrt(np.clip(np.diag(sigma2 * xtx_inv), 0.0, None))
    return OLSResult(names, beta, se, residuals, n, k)

channel のような文字列の列は、Rの lm() が自動で行うのと同じように ダミー変数 に展開する。水準が「Multichannel」「Phone」「Web」の三つなら、最初の「Multichannel」を基準にして channelPhone と channelWeb の二列を作る。三つすべてを列にすると切片と完全に重なってしまい、解が一つに決まらないからである。

Rの操作との対応は次の通りである。

バイアスのあるデータで回帰を実行する

データの準備は第1章と同じである。男性向けメール群とメールなし群に絞った male_df がRCTデータであり、そこから「購買意欲が高そうな人を介入群に多く残す」抽出を行った biased_df がバイアスのあるデータである。

rng = np.random.default_rng(1)
biased_df = (
    male_df.with_columns(
        ((pl.col("history") > 300) | (pl.col("recency") < 6)
         | (pl.col("channel") == "Multichannel")).alias("high_intent"),
        pl.Series("random_number", rng.random(male_df.height)),
    )
    .with_columns(
        pl.when(pl.col("treatment") == 1)
        .then(pl.when(pl.col("high_intent")).then(1.0).otherwise(0.5))
        .otherwise(pl.when(pl.col("high_intent")).then(0.5).otherwise(1.0))
        .alias("keep_probability")
    )
    .filter(pl.col("random_number") < pl.col("keep_probability"))
)

biased_reg = ols(biased_df, "spend", ["treatment", "history"])
print(biased_reg.tidy)

推定するモデルは、目的変数に購買額、介入変数にメールの有無、共変量に過去の購買額を使った次の式である。

実行結果は次の通りである。

教科書のRの結果では treatment の推定値は0.9026、標準誤差は0.1743であった。PythonとRでは乱数生成器が異なるため、同じシード1でも抽出される行が違い、値が少しずれる。しかし読み方は同じである。介入変数 treatment の推定値は約0.85であり、p値は非常に小さいので、「メール配信の効果はない」という帰無仮説は棄却される。この値は「メールを送ることで売上が平均0.85ほど増加する」と解釈できる。

なぜこの式が成り立つのか

教科書では「知られている」とだけ書かれているが、導出は簡単なので確認しておく。モデルBに、 $X_{\text{omit}}$ の回帰式を代入する。

最後の括弧の中は、 $Z_i$ と相関しない誤差である（ $\epsilon_i$ は $Z$ に回帰した残差なので $Z$ と無相関、 $e_i$ はモデルBの誤差なので $Z$ と無相関）。したがってこの式は「 $Y$ を $Z$ だけに回帰したモデル」の形になっており、その $Z$ の係数は $\beta_1+\beta_2\gamma_1$ である。これがモデルAの $\alpha_1$ にほかならない。

この余分な項 $\gamma_1\beta_2$ を 脱落変数バイアス（Omitted Variable Bias、OVB） と呼ぶ。モデルAで効果を分析すると、得られる推定結果は「本来の効果 + OVB」になってしまう。一方、 $X_{\text{omit}}$ を加えたモデルBでは、OVBが取り除かれる。

OVBの直観

$\beta_2$ はモデルBにおける $X_{\text{omit}}$ と $Y$ の関係であり、 $\gamma_1$ は $Z$ と $X_{\text{omit}}$ の関係である。つまりOVBは、「 $X_{\text{omit}}$ と $Y$ の関係」に「 $Z$ と $X_{\text{omit}}$ の関係」を掛けたものである。直観的には、 省略された共変量が $Y$ に与える影響が、 $Z$ との相関を通じて、あたかも $Z$ の効果のように見えてしまっている ということである。

ワークド例：Huluの8人でOVBを手計算する

抽象的な式を、小さなデータで確かめる。2.1の冒頭で層別比較に使った8人の表をもう一度使う。 $X$ は「アニメ視聴者か（1/0）」、 $Z$ は「『名探偵コナン』を上段トレイで訴求したか（1/0）」、 $Y$ は「7日間のコナン視聴分数」である。訴求はアニメ視聴者に集中させたので、アニメ視聴者4人のうち3人、非視聴者4人のうち1人が訴求を受けたとする。真の効果は全員 $+20$ 分で、訴求がなければアニメ視聴者は60分、非視聴者は20分見るとする。誤差はない。

モデルA（単純比較） ：訴求あり4人の平均は $(80+80+80+40)/4=70$ 、なし4人の平均は $(60+20+20+20)/4=30$ である。 $Z$ だけの単回帰の係数は平均差に等しいので、 $\alpha_1=70-30=40$ となる。真の効果20に対して、2倍に見えている。

モデルB（ $X$ を追加） ：データは誤差なしで $Y=20+40X+20Z$ に従っているので、回帰は $\beta_1=20$ 、 $\beta_2=40$ を返す。 $\beta_2=40$ は「アニメ視聴者は非視聴者より、訴求の有無をそろえて40分多く見る」という関係である。

モデルC（ $X$ を $Z$ に回帰） ：訴求あり4人のうちアニメ視聴者は $3/4=0.75$ 、なし4人のうちは $1/4=0.25$ である。したがって $\gamma_1=0.75-0.25=0.5$ となる。「訴求を受けた群には、アニメ視聴者が50ポイント多い」という関係である。

OVBの式を確かめる。

一致した。単純比較で得た40は、真の効果20に、「訴求群にアニメ視聴者が多い（ $\gamma_1=0.5$ ）」と「アニメ視聴者はもともと40分多く見る（ $\beta_2=40$ ）」の積20が上乗せされたものである。第1章では同じ差を「効果 + セレクションバイアス」と言葉で分解したが、OVBの式はそのセレクションバイアスを、二つの相関の積として定量化してくれる。

この例はPythonでも確認できる。

hulu8 = pl.DataFrame({
    "anime_fan":   [1, 1, 1, 1, 0, 0, 0, 0],
    "recommended": [1, 1, 1, 0, 1, 0, 0, 0],
    "minutes":     [80, 80, 80, 60, 40, 20, 20, 20],
})
model_a = ols(hulu8, "minutes", ["recommended"])
model_b = ols(hulu8, "minutes", ["recommended", "anime_fan"])
model_c = ols(hulu8, "anime_fan", ["recommended"])
alpha1 = model_a.coef_of("recommended")            # 40.0
beta1, beta2 = model_b.coef_of("recommended"), model_b.coef_of("anime_fan")  # 20.0, 40.0
gamma1 = model_c.coef_of("recommended")            # 0.5
print(gamma1 * beta2, alpha1 - beta1)             # 20.0 20.0

有意でない共変量を外してはいけない理由

2.1.6で「共変量の検定結果は解釈しない」と述べた。OVBの式を見ると、その理由がより強い形で分かる。OVBの値 $\gamma_1\beta_2$ の中には、「 $\beta_2$ が統計的に有意か」という情報はどこにも入っていない。有意でなくても $\beta_2$ が0でなければ、その変数を外したときにOVBが発生する。したがって、「有意でないから」という理由で共変量を除外すると、バイアスを含んだ推定結果を得ることになる。共変量の選択は、検定結果ではなく、次に述べる「介入の割り当てがどう決まったか」という知識に基づいて行う。

共変量の選択とCIA

回帰分析のモデルは、次のステップで作る。

1. 介入の割り当てがどのように決定されたかを考える。 たとえば「去年の購買量が多い人にメールを送っただろう」「最近購買した人に送っただろう」と想定する。Huluなら「アニメ視聴日数が多いプロフィールに訴求しただろう」「新規会員は対象外だろう」と想定する。
2. 想定される決定方法を表現できる共変量を選ぶ。 昨年の購買量、直近の購買。Huluなら過去90日のジャンル別視聴日数、契約からの経過月数。
3. 選んだ共変量と $Y$ の関係を考慮して、それぞれの関数の形を決める。 線形か、対数か、2乗項が必要か、区分のダミーにするか。

推定された効果の妥当性を主張するには、作ったモデルの共変量がCIAを満たしていると考える必要がある。しかし、これを主張するうえで二つの問題がある。

問題1：バイアスの評価ができない

一つは、 得られた推定値にどれだけバイアスが残っているかを評価する方法がない という問題である。もし評価できるなら、共変量の組み合わせを試してバイアスが最小になるものを選べばよい。しかし現実にはそれができない。

OVBはバイアスの評価になっているように見えるかもしれない。しかし、OVBは「モデル間でのバイアスの変化」を示すものであり、「残っているバイアスの大きさ」を示すものではない。上の例で history を追加してOVBが0.041減ったことは分かるが、history を追加した後のモデルにどれだけバイアスが残っているかは分からない。

このため、分析者は「セレクションバイアスがどのような理由で発生し、どの変数を含めればコントロールでき、どの組み合わせでCIAが満たされるか」を、データではなく知識に基づいて仮定する必要がある。

問題2：必要な共変量がデータにない

もう一つは、 手持ちのデータに含まれる変数だけでは、バイアスを十分に減らせない可能性がある という問題である。バイアスを生む要因を理解していても、それを表すデータが定義できない、あるいは定義できても入手できない場合がある。

たとえばHuluで、あるトレイの訴求対象が「運用担当者の目視による選定」で決まっていた場合、担当者の判断基準は記録されていない。あるいは、「そのプロフィールがコナンをどれくらい楽しみにしているか」は、視聴履歴から一部は推測できても、完全には測れない。

必要な共変量がデータにない場合、 回帰分析でセレクションバイアスを取り除くことは基本的に不可能である 。得られた推定値は、投入した共変量による影響は取り除けているものの、データにない変数によるOVBを含んでいる。ドメイン知識があれば、そのOVBの大きさや向きを推察することはできるが、それは主観を含む方法であり、分析者の願望が入り込みやすい。

これら二つの問題は、明確な指標を見ながらモデルを選ぶことができず、モデルの限界も定量的に評価できないことを意味する。したがって分析者は、手持ちのデータにどのようなバイアスがあり、選んだ共変量がどのバイアスをコントロールしているのかをよく認識して判断する必要がある。

なお、より応用的な手法を使えば、これらの問題に対処できる場合がある。操作変数法や固定効果モデルは、データが一定の条件を満たすときに、観測できない交絡があっても効果を推定できる。第4章の差分の差分法も、必要な共変量が手に入らないときの対処法の一つである。

Sensitivity Analysis：結果がどれだけ頑健かを確かめる

手持ちのデータにない変数がセレクションバイアスを起こしている可能性は、常につきまとう。回帰分析はデータにある変数しか使えないので、この問題には直接対処できない。したがって、そのようなバイアスが起きていそうかどうかを確認することが重要になる。

そのための方法が Sensitivity Analysis（感度分析） である。経済学の実証論文では、複数の共変量の組み合わせでの推定結果が並べて報告されることが多い。これは、組み合わせごとに推定された効果が大きく変化するかを見るためのものである。

具体的には、分析者が重要だと認識している共変量以外の共変量をモデルから抜いてみて、効果の推定値が大きく変動しないかを確認する。変動が小さければ、回帰の結果は他の変数の影響を受けにくいことを示しており、データにない変数を含めたとしても大きく変化しないだろうと示唆される。逆に、共変量の出し入れで推定値が大きく動くなら、まだ入れていない変数によっても大きく動く可能性が高く、結果の信頼性は低い。

本ノートの例で言えば、treatment の推定値は共変量なしで0.954、recency と channel を入れて0.860、さらに history を入れて0.819と、追加するごとに動いている。これは「まだ動く余地がある」ことの示唆であり、RCTの0.770との差がまさに残っているバイアスである。感度分析をより定量的に行う方法は、追加手法の節で紹介する。

仕組み：来訪者の「中身」が群によって違う

ポイントは、メールの配信が もともと購買傾向が弱いユーザの来訪を増やしている ことにある。

• メール配信群の来訪者は、「メールがなくても来訪する購買傾向の強いユーザ」と「メールがあったから来訪した購買傾向の弱いユーザ」の混合である。
• メール非配信群の来訪者は、「もともと購買傾向の強いユーザ」だけである。

したがって、来訪したユーザ同士で比較すると、配信群の方が購買傾向の弱い人を多く含むため、平均売上が低く見える。visit をモデルに入れることは、まさに「来訪者同士で比較する」ことに相当するので、効果が過小に推定される。

ワークド例：数値で確かめる

架空の数値で確認する。RCTで1,000人ずつを配信群と非配信群に分けたとする。ユーザは2種類おり、それぞれ半数ずつである。

真の効果を計算する。高意欲ユーザにはメールの効果がなく、低意欲ユーザはメールによって来訪し2購買するので、一人当たりの平均効果は $0.5\times0+0.5\times2=1$ である。全員で平均を比較すれば、配信群 $(500\times10+500\times2)/1000=6$ 、非配信群 $(500\times10+500\times0)/1000=5$ で、差は正しく $+1$ になる。

ところが来訪者に絞ると、次のようになる。

差は $6-10=-4$ である。真の効果は $+1$ なのに、来訪者同士の比較では $-4$ と、符号まで逆になった。 介入によって「誰が来訪者になるか」が変わってしまったため、来訪者という条件で揃えることが、比較可能性を壊してしまった のである。

群ごとの平均購買額を「全員で比べた場合」と「来訪者だけで比べた場合」で並べると、同じデータから正反対の結論が出ることが一目で分かる。

実際のRCTデータでも同じ構造が見える。来訪率は非配信群10.6%、配信群18.3%であり、メールは来訪を大きく増やしている。そして来訪者に絞った平均購買額は非配信群6.15、配信群7.78である。この差1.63は、全員で比較した効果0.77とは別の量であり、「来訪者の中身が違う」ことの影響を含んでいる。

このように、介入によって影響を受けた変数を分析に入れることで起きるバイアスを Post Treatment Bias と呼ぶ。重要なのは、 これは介入がどのように割り当てられていたとしても発生する ということである。上のワークド例はRCTの設定であった。つまり、メールをランダムに配信していても、そのデータで visit を投入した回帰を行えば、同じ問題に直面する。

対処：介入より前に値が決まる変数だけを使う

この問題を避けるには、 介入よりも後のタイミングで値が決まる変数を分析から除外する 必要がある。介入が割り振られる前に値が決まっている変数は、介入の影響を受けることがないので、この問題は起きない。ただし、どの変数がこれに該当するかの判断は分析者の知識に依存する。

Huluの推薦施策では、この落とし穴に特に注意が必要である。訴求の効果を測るときに、「訴求後のセッション数」や「トレイをクリックしたか」を共変量に入れると、まさに visit と同じ問題が起きる。訴求は「もともとあまりアプリを開かない人」の起動やクリックを増やすので、これらで条件付けると効果が過小に見える。

データの取得範囲によるPost Treatment Bias

Post Treatment Biasのような問題は、変数の選択以外でも発生する。実際のビジネスでは、 サービスを利用したユーザのログしか残っていない ことがよくある。たとえば、どのユーザが広告を目にしたかは分からないが、広告をクリックして流入したユーザのログは残っている、という状況である。

このとき、サービスへのアクセスを誘発する広告施策を行うと、広告配信グループには購買意向の弱いユーザが多く流入し、非配信グループはもともと購買意向の強いユーザだけがアクセスしている、という状態になる。このデータで分析すれば、先ほどのメールの例と同じく「介入グループの方が平均売上が低い」という結果を得る。割り振りがランダムであっても同じである。

自社サービス内の施策を評価する場合は、すべてのユーザのデータが得られることが多いので、この問題は起きにくい。しかし、自社サービスにユーザを誘引する施策（プッシュ通知、外部広告、メール）を評価する場合、分析するデータの取得は施策を実行するツールに依存するため、大きな問題になる。

Huluで言えば、プッシュ通知の効果を「通知後にアプリを起動した人」のログだけで評価すると、この問題に直面する。通知を受けた群には「通知がなければ起動しなかった人」が混ざり、受けなかった群には「自発的に起動する人」だけが残るからである。 施策を受けたユーザだけでなく、施策の対象になり得たユーザ全体のデータを用意すること が、分析において決定的に重要である。第1章で「視聴ログから始めず、割当台帳を基準に結合する」と述べたのは、この理由による。

私立学校への入学と奨学金の利用

図2.5は、上が女子生徒、下が男子生徒の結果である。

PRSCHA_1 では、男子生徒は約9ポイントの有意な効果があり、当選によって私立学校へ通い始める傾向が強くなっている。一方、女子生徒は2ポイントで有意ではないので、「私立学校へ通わせる効果はない」という可能性を否定できない。仮に有意だったとしても効果量は男子の4分の1以下であり、性別によって効果が大きく異なる。

USNGSCH では、男子生徒は約47%、女子生徒は約54%の効果であり、女子生徒の方が当選した際に奨学金を利用する傾向が強い。

通学の継続と留年

次に私立学校へ通い続ける効果と留年への効果を見る（図2.6）。

PRSCH_C の効果は男子13%、女子17%であり、私立学校への通学を続けさせる効果は女子の方が高い。 INSCHL は男女ともに有意な傾向がない。

REPT6 では、男子は当選グループで留年が8.7ポイント低く有意だが、女子は3.6ポイントの減少で有意ではない。男子には留年を減らす効果があるが、女子については効果がないことを否定できず、あったとしても男子より小さい。

FINISH6〜8 では、男子は8年生の修了を10ポイント増やす効果があり、6、7年生では有意差がない。女子は6年生の修了が3ポイント増え、8年生の修了は10ポイント増加している。女子も男子と同様に、より高い学年を修了できるようになっている。

男子の結果から描かれるストーリーは、全体のデータでの分析と変わらない。しかし女子の結果は大きく違う。 私立学校への通学を維持する効果は女子の方が高いのに、留年に対する効果はほぼない 。これは、女子生徒が私立学校へ通い続けられない原因が、学力や留年以外の要因にあることを示唆している。

一つの仮説は、女子生徒には家庭内の経済的リソースが分配されにくく、経済的なショックの影響を受けやすいというものである。急な所得減少が起きると、家庭内の男子よりも女子の方から公立への転校を迫られたり、労働によって家計を助けることを求められる可能性がある。この仮説が正しければ、割引券によって私立学校への通学が維持されるとともに、 労働時間の減少 が起きると予想される。

労働時間

労働時間 HOURSUM に対する結果を確認する（図2.7）。

男子生徒では0.6時間の減少で有意ではない。一方、女子生徒では2.1時間の減少であり、効果量が大きいうえに有意である。割引券によって経済的な問題がいくらか解決され、その結果として女子生徒が私立学校へ継続して通える状況になった、という解釈と整合している。

図2.5〜2.7の内容を、本ノートの推定値で一枚にまとめると次のようになる。女子（赤）と男子（青）を横に並べているので、どのアウトカムで男女の差が大きいかを見比べやすい。

Huluに置き換える：推薦施策の「しくみ」を追う

たとえば「新シーズン配信に合わせた上段トレイでの訴求」の効果を評価するとき、最終的なKPIである「7日以内の視聴」だけを見るのではなく、割引券の分析と同じ順序で複数のアウトカムを見ることができる。

「訴求で当該作品の視聴は増えたが、総視聴時間は変わらなかった」という結果が出れば、それは他作品からの置き換えが起きたことを示唆する。「新規会員では効果が大きく、継続会員では小さい」という結果が出れば、対象の絞り込みに使える。このように、 介入がどの経路で効いたかを複数のアウトカムから推理する ことが、探索的な効果検証の目的である。

ワークド例：近似の精度を確かめる

$X_2=100$ 、 $X_1=110$ （10%の増加）のとき、 $\log(110)-\log(100)=\log(1.1)\approx0.0953$ であり、変化率0.10とほぼ等しい。一方、 $X_1=150$ （50%の増加）では $\log(1.5)\approx0.405$ となり、0.50からずれが大きくなる。 対数の差による近似は、変化率が小さいときに正確である 。介入効果が数%〜十数%の範囲なら問題ないが、効果が数倍になるような場合は $\exp(\beta)-1$ で正確な変化率に戻すべきである。

次のモデルで解釈を確認する。

• $\beta_Z=0.1$ なら、介入された場合に目的変数が平均的に約10%高い。
• $\beta_1=0.01$ なら、 $X_1$ が1%増加したときに目的変数が約0.01%増加する。
• $\beta_2=0.1$ なら、 $X_2$ が1増加したときに目的変数が約10%高くなる。

本書の分析で対数を取る理由は二つに絞られる。一つは、 目的変数に対する介入の効果が比率で扱われるべきとき である。たとえばHuluで「訴求が総視聴時間を増やすか」を見る場合、月100時間見るヘビーユーザと月2時間のライトユーザに同じ「+30分」を期待するのは不自然で、「+5%」のように比率で考える方が自然なことが多い。もう一つは、 共変量と目的変数の関係が比率で扱われるべきとき である。本来は比率で影響する共変量を線形のまま入れると、共変量の関数の形が誤っているために、セレクションバイアスが残る可能性がある。

なお、視聴時間のように0を多く含む変数に対数を取るときは、 $\log(0)$ が定義できないため工夫が必要である。 $\log(1+Y)$ とすることが慣習的に行われるが、この変換は比率の解釈をゆがめる。0/1の視聴有無をアウトカムにするか、0を含めた線形モデルで効果を「分」で解釈する方が、無理がないことも多い。

ワークド例：分散の膨らみ方

他の変数と無相関（ $R_k^2=0$ ）の場合を基準にする。

$R_k^2=0.9$ なら分散は10倍、標準誤差は $\sqrt{10}\approx3.16$ 倍になる。信頼区間の幅も3倍以上に広がるので、本当は効果があっても有意にならなくなる。

この倍率を $R_k^2$ の関数として描くと、0.9を超えたあたりから急激に立ち上がることが分かる。

重要なのは、 多重共線性は関連の強い変数以外には影響を与えない という点である。したがって、介入変数以外の共変量同士で多重共線性が起きていても、興味のある $\beta_Z$ には影響がない。共変量の係数は解釈しないので、その標準誤差が大きくても問題にならない。一方、 介入変数と強く相関する共変量がある場合は $\beta_Z$ の標準誤差が大きくなる ので、注意が必要である。

Huluでこれが起きるのは、たとえば訴求のルールが「過去90日のアニメ視聴日数が30日以上」というしきい値だけで決まっていて、その視聴日数を共変量にも入れる場合である。 $Z$ が $X$ からほぼ決まってしまうと、 $Z$ と $X$ の相関が非常に高くなり、 $\beta_Z$ の標準誤差が大きくなる。この状況は「 $X$ が同じ人の中で $Z$ が異なる人がほとんどいない」ということでもあり、次節の分散重み付けの話とつながる。そして、しきい値で介入が決まる状況を積極的に利用する手法が、第5章の回帰不連続デザインである。

残差を使った表現（FWL定理）

重回帰における介入変数の係数 $\beta_Z$ は、次のように書ける。

ここで $\epsilon_i$ は、介入変数 $Z_i$ を他の共変量で回帰したときの残差である。

つまり $\epsilon_i$ は、 $Z_i$ のうち 共変量では説明できない部分 だけを残したものである。この残差と $Y$ の共分散を、残差の分散で割った値が重回帰の $\beta_Z$ になる。この事実はFrisch–Waugh–Lovell（FWL）定理と呼ばれ、重回帰の最小化問題を解くことで導出できる。

直観的には、重回帰は次の2段階を一度に行っている。まず「共変量から予測できる分の介入」を取り除き、次に「予測できなかった分の介入（ $\epsilon$ ）」と $Y$ の関係を単回帰で見る。共変量で完全に決まってしまう人（ $\epsilon_i=0$ ）は、2段階目に何も寄与しない。

Pythonで確かめる。

first_stage = ols(biased_df, "treatment", ["recency", "channel", "history"])
eps = first_stage.residuals
y = biased_df.select("spend").to_numpy().ravel()
beta_fwl = np.cov(y, eps, bias=True)[0, 1] / np.var(eps)
print(beta_fwl)   # 0.818992 = 重回帰の treatment 係数

結果は0.818992であり、重回帰 spend ~ treatment + recency + channel + history の treatment の係数と完全に一致する。

飽和モデルでの重み付け平均

共変量がすべて0/1のダミー変数で表され、その組み合わせごとに別々の平均を許すモデルを 飽和モデル（saturated model） と呼ぶ。このとき、上の式は次の形に変形できる。

この式の読み方は次の通りである。 $E[Y\mid Z=1,X]-E[Y\mid Z=0,X]$ は「共変量が $X$ の人たちの中での、介入ありとなしの平均の差」、つまり層ごとの効果である。 $\beta_Z$ は、 層ごとの効果を、その層の介入変数の分散 $\operatorname{Var}(Z\mid X)$ で重み付けして平均したもの に等しい。

介入変数の分散は $P(Z=1\mid X)(1-P(Z=1\mid X))$ なので、 $P(Z=1\mid X)=0.5$ のとき最大の0.25になり、0または1のときは0になる。

重みの形と、それが回帰係数にどう影響するかは、次のワークド例のあとに図で確認する。

したがって、 介入グループと非介入グループの割合が50:50に近い層が、回帰の結果にもっとも影響を与える 。一方、どちらか一方しか含まれない層は推定結果にほぼ影響しない。直観的には、回帰分析は「与えられたデータの中で、より適正な比較が行われていそうな層」の差を効果としてとらえている。

ワークド例：層ごとに効果が違うと、回帰係数はATEとずれる

Huluの架空データで確かめる。共変量 $X$ は「ヘビーユーザか（1/0）」、 $Z$ は訴求の有無、 $Y$ は7日間の視聴分数とする。

ヘビーユーザは訴求がなくてもたくさん見るので効果は小さく（+10分）、ライトユーザは訴求で初めて見るので効果が大きい（+30分）とする。人数は同じなので、 ATE（全員の平均効果）は $(10+30)/2=20$ 分 である。

一方、回帰 minutes ~ recommended + heavy_user の係数は、分散で重み付けした平均になる。

Pythonでこのデータを作って回帰すると、係数は15.2941となり、手計算と一致する。ATEの20より小さい。 訴求が半々のヘビーユーザ層の効果（+10分）が強く反映され、訴求がほとんどないライトユーザ層の効果（+30分）は軽く扱われた からである。

これが、回帰分析の結果とRCTの結果が同一にならない一つの理由である。RCTでは、どの $X$ においても介入割合の期待値は0.5なので、すべての層の重みが等しくなり、回帰はATEを推定する。しかしRCTでないデータでは、介入割合が層によって異なるので、回帰は必ずしもATEを推定していない。ただし重要なのは、 これはセレクションバイアスがうまく取り除けていないという話ではない という点である。回帰分析が「妥当に比較できそうなサンプルに絞り込んで比較する」という特徴を持っていることが原因である。効果が層間で大きく変わらないなら、この問題は起きず、RCTと同じ結果が得られる。

手法ごとに「推定している効果」は違う

因果推論の手法の多くは、手持ちのデータに何らかの重みを付けて妥当な比較を行う構造になっている。どのような重みを使うかによって、どのような効果が推定されているかが異なる。回帰は $\operatorname{Var}(Z\mid X)$ で重み付けしたもの、第3章の傾向スコアマッチングは介入群の分布に合わせたもの（ATT）、逆確率重み付けは母集団全体に合わせたもの（ATE）というように、狙う効果が違う。したがって、複数の因果推論の方法を使った場合、それらが同じ値を返すことはほとんどない。値が違うからといって、どれかが間違っているとは限らない。

Huluの施策評価で言えば、「ヘビーユーザとライトユーザで訴求の効果が違う」ことが予想されるなら、回帰係数を一つ報告するだけでは不十分である。層別に効果を推定するか、次の追加手法の節で述べる交互作用を入れたモデルを使い、「誰に対する効果か」を明示する必要がある。