1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

因果推論 Day 28/全30回 空間データ×因果推論、隣が効くとき何が壊れるか

1
Posted at

この連載について

因果推論を「本を読んだ」で終わらせず、自分の言葉で説明でき、コードで再現できる状態まで落とす30日連載です。直前のDay 27では、26日分の道具を出す前の総点検として15項目のチェックリストに畳みました。その3番(SUTVA)と7番(オーバーラップ)が同時に壊れる場所が、私の主戦場である空間データです。

今日はそこへ踏み込みます。私は衛星データと水道インフラの会社で機械学習をやっていて、扱うデータはほぼ全部が「位置」を持ちます。位置が付いた瞬間に、隣り合う単位は似た未観測要因を共有し、施策の効果は境界を越えて漏れます。データはすべて格子状の合成自治体で、真の効果を私たちが握った状態で壊し方と直し方を見ます。


TL;DR(3行)

  • 空間交絡は、未観測交絡が地図の上でつながっている状態です。 真値2.00に対して素朴な差は 4.25、観測共変量Xで調整しても4.16、2次の空間トレンドを足しても4.05。効いたのは空間ブロック固定効果の 2.41 で、勝負を決めるのは「調整するかどうか」ではなく「どのスケールで調整するか」でした。
  • 干渉があると、推定より先に定義が壊れます。 露出(自分の処置+近傍の処置割合)で書き直すと、直接 1.47・波及 1.02・全面展開の総効果 2.48 に分解できます。素朴な差1.97はどれでもない中間の量で、しかもユニット無作為化しても出るのは直接効果1.51だけでした。
  • 未観測自治体への汎化は、ランダムCVが嘘をつきます。 ランダム4分割CVは誤差 0.87 と報告しますが、別地域へ持ち出した真の誤差は 1.92。バッファ付き空間ブロックCVだけが1.92を当てました。

今日の問い

隣の自治体が施策を入れると、うちの数字も動く。効果は境界を越えて漏れ、似た者同士が隣り合って並んでいる。
こういう世界で「この施策の効果は◯◯でした」と言うとき、◯◯とは何の数字なのか?

Day 2で潜在的結果を定義したとき、SUTVA(他人の割付が自分の結果に影響しない)を暗黙に置きました。Day 8ではマッチングの前提として仮定で封じ、「崩れる世界はDay 28で」と先送りしました。今日がその回です。壊れ方は2つあり、片方は推定のズレ(空間交絡)、もう片方は定義の崩壊(空間干渉)です。後者のほうが厄介で、推定を工夫する前に何を推定するのかを書き直すところから始まります。


概念① 空間交絡:滑らかな未観測交絡は、共変量調整では消えない

空間データの出発点は、Toblerの第一法則です。「すべてのものは関連しているが、近いものはより強く関連している」。地価も気温も人口動態も財政力も、隣り合う自治体では似ています。これを空間的自己相関と呼びます。

因果の言葉に翻訳すると、こうなります。未観測交絡 $U$ が位置 $s$ の滑らかな関数になっている状態、つまり

$$
U_s \approx U_{s'} \quad (\text{$s$ と $s'$ が近いとき})
$$

空間交絡です。DAG(Day 10)に描けば、$U \to W$ と $U \to Y$ の分岐が1本あるだけで、構造そのものは Day 11 のバックドアと同じです。違うのは規模と姿です。$U$ は個体ごとにバラバラな乱数ではなく、地図の上でひとつながりの模様をしています。位置が近いというだけで交絡が伝播するので、実質的には全個体を巻き込む1本の巨大な分岐になります。

厄介なのは、$U$ が観測できないことです。ではどうするか。空間交絡の場合だけ使える抜け道があります。$U$ が位置の滑らかな関数なら、位置そのものが $U$ の代理になるからです。手は3段階あります。

  • 大域的な空間トレンド:座標の多項式を回帰に足し、南北や東西のゆるやかな勾配を除く。Day 16 のFWL定理でいう「座標で残差化してから見る」に当たります。
  • 局所的な空間固定効果:地域をブロックに区切り、ブロックダミーを入れる。ブロック内の比較だけで効果を測るので、ブロックより大きなスケールの交絡はすべて落ちます。
  • 空間ラグや空間誤差モデル:隣接構造を明示的にモデル化する。空間計量経済の王道ですが、今日は前2つに絞ります。

ただし、うまい話ではありません。位置で調整するということは、処置の空間的な変動も一緒に削ることを意味します。ブロックを細かくするほど交絡は落ちますが、ブロック内で処置が全員同じになった瞬間、比較相手が消えて効果は識別できなくなります。この綱引きが、手を動かす①の主題です。

概念② 空間干渉:SUTVAが壊れると、個体の処置効果の定義から書き直す

もうひとつの壊れ方が干渉(interference)、いわゆるスピルオーバーです。上流の自治体が漏水対策を進めれば下流の水需給が変わり、隣接自治体が共同発注に加われば単価が下がる。自分が受けていない処置の効果が、自分の結果に乗ります。

ここで壊れるのは推定ではなく定義です。Day 2 で書いた潜在的結果 $Y_i(w_i)$ は、「自分の割付だけで自分の結果が決まる」という前提の上に立っていました。干渉があると、正しい書き方は割付ベクトル全体を引数に取る

$$
Y_i(w_1, w_2, \ldots, w_N)
$$

になります。$N$ 個の自治体に対して、1人あたり $2^N$ 通りの潜在的結果です。自治体が1700あれば $2^{1700}$ 通りで、これは推定以前に定義として使い物になりません。

そこで**露出マッピング(exposure mapping)**を置きます。割付ベクトル全体を、低次元の「露出」に要約する関数 $f$ を仮定するのです。今日は素直に

$$
f(i, \mathbf{w}) = \bigl(w_i,; g_i\bigr), \qquad g_i = \frac{1}{|N(i)|}\sum_{j \in N(i)} w_j
$$

とします。$w_i$ が自分の処置、$g_i$ が隣接自治体のうち処置を受けた割合です。この要約で十分だと仮定する、つまり $Y_i(\mathbf{w}) = Y_i(w_i, g_i)$ と書けると仮定すると、潜在的結果は2引数に戻り、推定対象を3つに分けて定義できます。

  • 直接効果 $Y_i(1, g) - Y_i(0, g)$:近傍の状況を $g$ に固定したまま、自分だけが導入したときの変化。
  • 波及効果(間接効果) $Y_i(0, g) - Y_i(0, 0)$:自分は導入しないまま、近傍が $g$ まで導入したときの変化。
  • 総効果 $Y_i(1, 1) - Y_i(0, 0)$:誰も導入していない世界と、全面展開した世界の差。政策決定でいちばん知りたいのは、たいていこれです。

露出マッピングは仮定であって発見ではありません。「干渉は8近傍まで」「隣の隣は効かない」「効くのは処置された数の割合だけで、誰が処置されたかは関係ない」という強い制約を置いています。SUTVAを捨てたのではなく、SUTVAを露出のレベルで置き直したというのが正確です。この置き直しが妥当かどうかは、水が流れる方向や発注の単位といったドメインの話で、データからは決まりません。


手を動かす①:空間交絡は、共変量調整でも空間トレンドでも消えない

80×80の格子(6400自治体)を作ります。未観測交絡 $u$ は scipy.ndimage.gaussian_filter でホワイトノイズをぼかして作り、空間的に滑らかにします。導入確率は主に $u$ で決まり、$u$ は成果 $y$ も押し上げます。真の効果は2.0に固定し、5通りの調整を比べます。

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy.ndimage import gaussian_filter


def smooth_field(rng, G, sigma):
    """空間的に滑らかな標準化ガウス場(隣り合うセルほど値が似る)"""
    z = gaussian_filter(rng.normal(size=(G, G)), sigma=sigma, mode="wrap")
    return (z - z.mean()) / z.std()


rng = np.random.default_rng(28)
G = 80                                        # 80×80=6400自治体の格子
u = smooth_field(rng, G, sigma=4.0)           # 未観測交絡:空間的に滑らか
x = rng.normal(size=(G, G))                   # 観測共変量:空間構造なし
e = 1 / (1 + np.exp(-(1.5 * u + 0.5 * x)))    # 導入確率は主に U で決まる
w = rng.binomial(1, e)
y = 1.0 + 2.0 * w + 2.0 * u + 0.5 * x + rng.normal(0, 1, (G, G))  # 真の効果 2.0

row, col = np.indices((G, G))
d = pd.DataFrame({
    "y": y.ravel(), "w": w.ravel(), "x": x.ravel(), "u": u.ravel(),
    "r": (row.ravel() / (G - 1)) * 2 - 1,     # 座標を[-1,1]に正規化
    "c": (col.ravel() / (G - 1)) * 2 - 1,
    "blk": (row.ravel() // 5) * 1000 + (col.ravel() // 5),   # 5×5セルのブロック
})

for name, f in [
    ("素朴な差", "y ~ w"),
    ("Xで調整", "y ~ w + x"),
    ("X+2次の空間トレンド", "y ~ w + x + r + c + I(r**2) + I(c**2) + r:c"),
    ("X+空間ブロック固定効果", "y ~ w + x + C(blk)"),
    ("Uを直接調整(オラクル)", "y ~ w + x + u"),
]:
    print("%s = %.2f" % (name, smf.ols(f, d).fit().params["w"]))
素朴な差 = 4.25
Xで調整 = 4.16
X+2次の空間トレンド = 4.05
X+空間ブロック固定効果 = 2.41
Uを直接調整(オラクル) = 2.01

素朴な差 4.25 は真値2.00の倍以上です。ここまでは Day 6 の交絡と同じ話で、驚きはありません。問題はその後です。手元の共変量Xで調整しても4.16、座標の2次多項式でトレンドを抜いても 4.05 と、ほとんど動きません。$u$ の模様は南北や東西の単調な勾配ではなく、数セル規模のまだら模様なので、大域的な多項式ではまるで捕まらないのです。

一方、5×5セルのブロック固定効果を入れると 2.41 まで落ちます。ブロック内なら $u$ はほぼ一定とみなせるので、ブロック内の導入済みと未導入を比べる形になり、交絡の大半が消えました。$u$ を直接調整したオラクルが2.01ですから、残っている0.4はブロック内に残った $u$ の変動です。空間交絡は「空間で調整したかどうか」ではなく、「交絡の空間スケールに合った粒度で調整したかどうか」で決まります。

下の地図が、この数字の正体です。左が未観測交絡、中央が導入の分布、右が成果。中央の地図でオレンジ(導入)のかたまりが、左の明るい領域とほぼ重なっているのが見えます。導入率は $u$ の高い帯で90%、低い帯で11%でした。これは同時に、Day 27 のチェックリスト7番(オーバーラップ)も壊れているという意味です。$u$ の高いかたまりの内部には、比較相手になる未導入自治体がほとんどいません。

Fig 1: 空間交絡の地図。未観測交絡Uの模様に、導入Wも成果Yも従う

Fig 2: 素朴4.25→X調整4.16→空間トレンド4.05→ブロック固定効果2.41→オラクル2.01(真値2.00)


手を動かす②:干渉を、直接効果と波及効果に分解する

次は干渉です。交絡の話と混ざらないよう、ここでは未観測交絡をゼロにします。成果は「自分の処置」と「8近傍の処置割合」の両方で決まる形にし、真の直接効果を1.5、波及の係数を1.0(隣が全部導入したら1.0だけ上がる)と置きます。割付は2通り試します。地理的にかたまる普及(現実の姿)と、1自治体ずつのコイン投げ(RCT)です。

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy.ndimage import convolve, gaussian_filter


def smooth_field(rng, G, sigma):
    z = gaussian_filter(rng.normal(size=(G, G)), sigma=sigma, mode="wrap")
    return (z - z.mean()) / z.std()


rng = np.random.default_rng(280)
G = 60
kern = np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]], float)   # 8近傍


def outcome(w, rng):
    g = convolve(w.astype(float), kern, mode="wrap") / 8    # 近傍の処置割合=露出
    y = 5.0 + 1.5 * w + 1.0 * g + rng.normal(0, 1, (G, G))  # 直接1.5・波及1.0
    return g, y


wc = (smooth_field(rng, G, sigma=1.0) > 0).astype(int)   # 導入が地理的にかたまる
gc, yc = outcome(wc, rng)
wr = rng.binomial(1, 0.5, (G, G))                        # 1自治体ずつコイン投げ
gr, yr = outcome(wr, rng)

for tag, w, g, y in [("空間クラスター割付", wc, gc, yc), ("ユニット無作為化  ", wr, gr, yr)]:
    print("%s: 素朴な差 = %.2f(近傍処置率の差 = %.2f)"
          % (tag, y[w == 1].mean() - y[w == 0].mean(), g[w == 1].mean() - g[w == 0].mean()))

d = pd.DataFrame({"y": yc.ravel(), "w": wc.ravel(), "g": gc.ravel()})
fit = smf.ols("y ~ w + g", d).fit()
print("露出で分解: 直接 = %.2f、波及 = %.2f、総効果 = %.2f"
      % (fit.params["w"], fit.params["g"], fit.params["w"] + fit.params["g"]))
空間クラスター割付: 素朴な差 = 1.97(近傍処置率の差 = 0.50)
ユニット無作為化  : 素朴な差 = 1.51(近傍処置率の差 = 0.02)
露出で分解: 直接 = 1.47、波及 = 1.02、総効果 = 2.48

3つの数字がすべて違います。順に読みます。

地理的にかたまって普及した世界では、素朴な差が 1.97 でした。導入済みの自治体は隣も導入済みであることが多く、近傍処置率の差が0.50あります。つまり素朴な差は、直接効果1.5に「隣が0.5だけ多く導入している分の波及」を混ぜた量です。真の直接効果でも真の総効果でもない、混合比が普及の地理的な形で決まる中間の量が出てきます。

では無作為化すればよいのか。ユニット無作為化の素朴な差は 1.51 で、真の直接効果1.5をきれいに当てました。近傍処置率の差が0.02とほぼゼロだからです。ここが今日いちばん誤解しやすい点だと思います。無作為化は交絡を消しますが、干渉は消しません。1自治体ずつコイン投げをして得られるのは、近傍がおおむね半分導入している状況での直接効果であって、全面展開したときの総効果 2.48 ではありません。差は0.97、率にして4割近くを取りこぼします。パイロットの数字をそのまま全国展開の見積もりに使うと、この分だけ過小評価します。

露出 $g$ を回帰に入れて分解すると、直接1.47・波及1.02となり、真値の1.5と1.0をそれぞれ復元できました。合計2.48が総効果です。ただしこれが成り立つのは、露出マッピングを正しく(8近傍の割合と)指定できたからで、干渉の範囲を読み違えれば分解ごと崩れます。

Fig 3: 直接1.47+波及1.02=総効果2.48。素朴な差はクラスター割付で1.97、無作為化で1.51


手を動かす③:ランダムCVは、未観測自治体への汎化に嘘をつく

3つめは評価の話です。効果が場所によって違う(Day 22のCATE)とき、推定した $\tau(x)$ をまだ観測していない自治体へ持ち出せるかが実務の関心事になります。ここで交差検証の切り方が効きます。

設定はこうです。真の効果は観測共変量 $x_1$ と、観測できない滑らかな地域要因 $z$ で決まります。モデルが使えるのは座標と $x_1$ だけ。割付はコイン投げなので交絡はなく、純粋に汎化の話に絞れます。T-learner(Day 23)で $\tau$ を推定し、真の $\tau$ との誤差(PEHE)をランダム4分割CVと空間ブロック4分割CVで測り、最後にまったく別の地域へ持ち出した真の誤差と突き合わせます。

import numpy as np
from scipy.ndimage import binary_dilation, gaussian_filter
from sklearn.ensemble import RandomForestRegressor


def smooth_field(rng, G, sigma):
    z = gaussian_filter(rng.normal(size=(G, G)), sigma=sigma, mode="wrap")
    return (z - z.mean()) / z.std()


def make_region(rng, G):
    x1 = rng.normal(size=(G, G))
    z = smooth_field(rng, G, sigma=5.0)             # 観測できない滑らかな地域要因
    tau = 1.0 + 0.8 * x1 + 1.2 * z                  # 真のCATE(場所で変わる)
    w = rng.binomial(1, 0.5, (G, G))                # 割付はコイン投げ(交絡なし)
    y = 2.0 + 1.0 * x1 + 1.5 * z + tau * w + rng.normal(0, 1, (G, G))
    row, col = np.indices((G, G))
    feat = np.c_[row.ravel(), col.ravel(), x1.ravel()]   # 使えるのは座標と x1 だけ
    return feat, w.ravel(), y.ravel(), tau.ravel()


def t_learner_tau(ftr, wtr, ytr, fte):
    rf = dict(n_estimators=200, min_samples_leaf=5, random_state=0)
    m1 = RandomForestRegressor(**rf).fit(ftr[wtr == 1], ytr[wtr == 1])
    m0 = RandomForestRegressor(**rf).fit(ftr[wtr == 0], ytr[wtr == 0])
    return m1.predict(fte) - m0.predict(fte)


rng = np.random.default_rng(2828)
G = 48
feat, w, y, tau = make_region(rng, G)
rnd_fold = rng.permutation(len(y)) % 4                          # ランダム4分割
blk = ((feat[:, 0] // 12) * 4 + feat[:, 1] // 12).astype(int)   # 12×12セル=16ブロック
blk_fold = rng.permutation(16)[blk] % 4                         # ブロックごと4分割


def cv_pehe(fold, buffer=0):
    err = []
    for k in range(4):
        te = fold == k
        tr = ~te
        if buffer > 0:                        # 検証ブロックの周囲を訓練から外す
            m = np.zeros((G, G), bool)
            m.ravel()[te] = True
            near = binary_dilation(m, structure=np.ones((2 * buffer + 1,) * 2, bool))
            tr = tr & ~near.ravel()
        err.append((t_learner_tau(feat[tr], w[tr], y[tr], feat[te]) - tau[te]) ** 2)
    return np.sqrt(np.concatenate(err).mean())


f2, _, _, tau2 = make_region(np.random.default_rng(999), G)   # まったく別の地域
print("ランダム4分割CV          : PEHE = %.2f" % cv_pehe(rnd_fold))
print("空間ブロック4分割CV      : PEHE = %.2f" % cv_pehe(blk_fold))
print("ブロック+バッファ5セル  : PEHE = %.2f" % cv_pehe(blk_fold, buffer=5))
print("未観測の別地域(真の誤差): PEHE = %.2f"
      % np.sqrt(((t_learner_tau(feat, w, y, f2) - tau2) ** 2).mean()))
ランダム4分割CV          : PEHE = 0.87
空間ブロック4分割CV      : PEHE = 1.12
ブロック+バッファ5セル  : PEHE = 1.92
未観測の別地域(真の誤差): PEHE = 1.92

ランダム4分割CVは誤差 0.87 と報告します。検証セルのすぐ隣に訓練セルがいるので、モデルは座標から地域要因 $z$ を局所的に埋め合わせられます。実際にやっているのは汎化ではなく空間内挿です。ところが同じモデルを別の地域へ持ち出したときの真の誤差は 1.92、2倍以上でした。ちなみに真の $\tau$ の標準偏差は1.41ですから、この予測は「全員に平均を答える」より悪いことになります。座標に頼って覚えた模様が、新しい土地では役に立たないどころか害になるわけです。

空間ブロックCV(1.12)は方向としては正しいものの、まだ楽観的でした。ブロックを丸ごと外しても、検証ブロックの縁は訓練データと接しているので、そこから染み出します。検証ブロックの周囲5セルを訓練から外したバッファ付きで 1.92、真の誤差とぴたり一致しました。バッファ幅5は、この実験で地域要因を作ったときの平滑化スケール($\sigma = 5$)に合わせています。実務では相関の届く距離を先に見積もる作業(バリオグラムなど)が要ります。

これは Day 24 で「fold の切り方が本丸」と書いた話の実測版です。DMLの交差適合をランダムKFoldでやれば、残差化の段階で同じ漏れが起きます。Day 25 のhonestyも同様で、分割用と推定用をランダムに割るだけでは、空間相関のある2つの標本は独立になりません。

Fig 4: ランダム分割と空間ブロック分割の地図、PEHEはランダム0.87・ブロック1.12・バッファ付き1.92・真の誤差1.92


つまづき・誤解しやすい点

  • 「座標を特徴量に入れたから空間を考慮した」ではありません。 手を動かす①では、2次の空間トレンドを入れても4.25が4.05にしか動きませんでした。逆にブロックを細かくすると2.41(5×5セル)→2.29(4×4)→2.16(3×3)→2.08(2×2)と真値に近づきます。ただしこれは「細かいほど良い」という話ではありません。細かくするほどブロック内に残る処置の変動が減り、極端にはブロック内が全員同じ処置になって比較相手が消えます。今回は格子が大きく処置も半々なので2×2セルでも変動が残り、標準誤差は0.034〜0.036でほぼ横ばいでしたが、現実の自治体データでこの粒度まで刻めば、まず先に比較が枯れます。調整の粒度は、交絡の空間スケールについての仮定そのものです。検証できない仮定という点で、Day 6 の非交絡性と同じ性格を持ちます。
  • 「干渉があるなら無作為化すればいい」ではありません。 ユニット単位で無作為化すると、交絡は消えますが干渉は残ります。手を動かす②の1.51が示すのは、無作為化実験が答えるのは「近傍が現状のままで自分だけが導入したときの効果」だという事実です。全面展開の総効果2.48を知りたいなら、割付の単位を変える(地理クラスター単位で無作為化する)か、露出を明示的にモデル化するしかありません。デザインを変えずに推定量だけ工夫しても届きません。
  • 空間ブロックCVにすれば正直、でもありません。 ブロックだけでは1.12と楽観が残り、バッファ5セルで1.92になりました。バッファを3セルにすると1.68で、まだ足りません。適切なバッファ幅は空間相関の届く距離に依存し、その距離自体がデータから推定する対象です。安全側に倒すならバッファを広く取りますが、今度は訓練データが痩せます。ここにも綱引きがあります。

GISデータ実務での使い方

  • 水道インフラの施策評価は、干渉の宝庫だと最初に宣言します。 上流の更新が下流の水圧や漏水に効き、広域企業団の共同発注は参加していない近隣の単価にも波及します。私は評価設計の1枚目に「効果はどこまで漏れるか」を書く欄を作り、露出変数(隣接自治体の導入割合、上流側の更新率)を最初から定義するようにしています。結果の報告も、直接効果と波及効果を分けた形にします。分けずに1つの数字で出すと、混合比が普及の地理的な形で決まってしまい、後から別の地域に当てはめられません。
  • メッシュの粒度を選ぶことが、交絡の仮定を選ぶことになります。 衛星由来の指標(地表面温度、地盤変動、土地被覆)は空間的に滑らかなので、メッシュ固定効果や地域ダミーでかなり吸えます。逆に言えば、粒度を決める議論は「どのスケールの未観測要因まで落とすか」の議論です。私は粒度を1段変えるたびに推定値がどう動くかを並べて出すことにしています。手を動かす①の2.41→2.29→2.16の並びは、そのまま感度分析(Day 9)の一種として読めます。数字が粒度で大きく動くなら、それは結論が仮定に強く依存している合図です。
  • モデルの精度は、バッファ付き空間ブロックCVの数字だけを提案書に載せます。 漏水リスクや劣化予測のモデルを「未導入の自治体へ展開する」文脈で語るとき、ランダムKFoldの精度は無意味です。手を動かす③のとおり2倍以上の乖離が出ます。特に自治体IDや座標を特徴量に入れたモデルは、新しい地域で平均予測より悪くなりえます。同じ話がDML(Day 24)の交差適合とCausal Forest(Day 25)のhonestyにも効くので、社内では「fold は空間で切る」を既定にしておくのが安全です。
  • パイロットは自治体単位ではなく、地理クラスター単位で割り付けます。 干渉があるとき、ユニット無作為化のパイロットは総効果を過小評価します。処置クラスタと対照クラスタの間に緩衝地帯を置く設計は、手を動かす③のバッファと同じ発想です。現場では「隣り合う自治体で片方だけ実施する」ほうが調整コストは低いのですが、その安さは推定対象のズレという形で後から請求されます。実施の可否を握るのは私ではないことが多いので、設計の相談を受けた段階で「その割り方だと答えられるのは直接効果だけです」と言えるようにしておくのが、この回でいちばん実務に効く準備だと思っています。

参考(本棚)

  • Imbens & Rubin『Causal Inference for Statistics, Social, and Biomedical Sciences』第1章 — SUTVAの定義と、干渉が入ると潜在的結果が $2^N$ 通りに膨らむ話
  • Hudgens & Halloran (2008) "Toward Causal Inference With Interference", JASA — 直接効果・間接効果・総効果の定義。部分干渉(クラスタ内では干渉するが、クラスタ間ではしない)の枠組み
  • Aronow & Samii (2017) "Estimating average causal effects under general interference", Annals of Applied Statistics — 露出マッピングの定式化。今日の概念②の下敷き
  • Roberts et al. (2017) "Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure", Ecography — 空間ブロックCVとバッファの実務
  • Facure『Causal Inference for the Brave and True』 Geo実験の章 — 地理単位で割り付ける実験設計
  • 自分用の空間機械学習サーベイノート(空間的自己相関・空間CV・空間特徴量) — 予測側の整理。今日はその因果版に当たります
  • 本文のコードはすべて合成データです。乱数シードを含めコピペで再現できます

次回予告(Day 29)

今日の3つの実験には、共通する居心地の悪さがありました。ブロックの粒度もバッファ幅も露出マッピングの範囲も、どれも「決め打ち」で、その決め方に自信があるわけではありません。点推定を1つ出して終わりにするより、分からなさを分布として持ち歩くほうが誠実な場面があります。次回Day 29は、ベイズと因果推論です。効果を点ではなく事後分布で語り、「効果がプラスである確率は87%」と直接言えるようにします。PyMCで小さなRCTを解き、事前分布を振って感度を見て、最後に自治体ごとの効果を階層モデルで部分プーリングします。サンプルの少ない自治体の極端な推定値が、どう縮むのかを見ます。

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?