因果関係の推定
名古屋大学 山本俊行

2021/09/20-21

第20回 行動モデル夏の学校

1

報告内容
• 導入：日本心理学会YouTubeライブ• ランダム化比較試験と傾向スコアマッチング• 動学モデルの適用

2021/09/17-19

第20回 行動モデル夏の学校

2日本心理学会での議論
日本心理学会YouTubeライブ・話題の論文について著者と語るシリーズ 第1回 吉田・村井(2021)「心理学的研究における重回帰分析の適用に関わる諸問題」（2021/7/2配信）• 観察研究のデータのみに基づく因果関係の確率は不可能である(岩崎, 2015)• 説明のためには変数間の因果関係についてのモデルが重要になる
3

議論の内容著者らが指摘する問題点：• 交絡変数を組み入れない問題• 逆方向の因果関係の存在に関する問題– 断面データでは因果関係の検証は不可能• 媒介変数を独立変数として組み入れてしまう問題(昨年度資料参照)

私なりの勝手な議論のまとめ：• 結果や行動の予測としては，汎化性能の評価や正則化が必須（過適合，推定に用いないデータに対する予測性能）– R2 ≥ .87 でないと予測の標準誤差が従属変数の標準偏差の半分以下にならない（吉田，2018）• 因果関係の探求としては，理論的背景が重要• 個人間差異だけでなく個人内変動が重要

4

因果関係の探求• ランダム化比較試験（デュフロ他，2019）– 母集団からのランダムな抽出や，処置群(with)と対照群(without)のランダムな割り当てを行う– 現実社会ではそんな実験は出来ないことも多い• マッチング– 処置群の各被験者に対して，対照群の中から同様の属性Xを持つ被験者を 1 人または複数見つける– マッチング手法として傾向スコア・マッチング等（Rosenbaum & Rubin, 1983）

2021/09/17-19

第20回 行動モデル夏の学校

5

傾向スコアマッチング
手順：(Wikipedia)
1. 処置群(y = 1)と対照群(y = 0)を被説明変数とする2項選択モデルを推定し，Pr(Y = 1|X)を傾向スコアとする
2. 傾向スコアで層別化し，層内の処置群と対照群の間で共変量Xのバランスが取れていることを確認
3. 傾向スコアに基づいて処置群の各被験者を対照群の被験者 1 人または複数とマッチング– 最近傍探索，キャリパーマッチング，マハラノビス距離マッチング，層化抽出法，差分マッチング，完全一致，等• King & Nielsen （2019）による批判– モデルの不均衡，非効率，モデル依存，バイアスを増加させる

2021/09/17-19

第20回 行動モデル夏の学校

6

ヘックマン補正との関係• ヘックマン補正（Heckman, 1974）– Pr(Y = 1|X) = Φ(γX)– Z = βX + ε，Z: 調べたいこと（処置による変化）– E[Z|X, Y = 1] = βX + E[ε|X, Y = 1]
= βX + ρσελ(γX)ρ: 相関係数，λ: 逆ミルズ比• ヘックマン補正は非観測異質性の補正⇔ マッチングは観測異質性の補正？
2021/09/17-19

第20回 行動モデル夏の学校

7

動学化へのアプローチ1：変動について（北村，2003）
交通は移ろいやすく，交通現象に変動はつきもの• 差異：一断面における個体間の変動• 変化：特定の個体の系時的変動• 変動：確率過程的変動
2007/09/20-21

第6回 行動モデル夏の学校

8

交通行動分析データの進展• 断面データ– 環境の異なる個人間差異• パネルデータ– 時点間の個人の行動変化• 長期観測データ– 個人内変動の考慮
2007/09/20-21

第6回 行動モデル夏の学校

9長期観測の利点• 断面データ（1時点）
y = f(x) + e(差異, 変化, 変動)• パネルデータ
yｔ = f(xｔ) + e(差異) + eｔ(変化, 変動)
yｔ’ = f(xｔ’) + e(差異) + eｔ’(変化, 変動)
yｔ’ - yｔ = f(xｔ’) - f(xｔ) + eｔ’(変化, 変動) - eｔ(変化, 変動)
2007/09/20-21

第6回 行動モデル夏の学校

10長期観測の利点• 長期観測yｔ1 = f(xｔ1) + e(差異) + eｔ(変化) + eｔ1(変動)yｔ2 = f(xｔ2) + e(差異) + eｔ(変化) + eｔ2(変動)yｔ’1 = f(xｔ’1) + e(差異) + eｔ’(変化) + eｔ’1(変動)yｔ’2 = f(xｔ’2) + e(差異) + eｔ’(変化) + eｔ’2(変動)E(yｔ’N) - E(yｔN) = E(f(xｔ’)) – E(f(xｔ))+ eｔ’(変化) - eｔ(変化)2007/09/20-21

第6回 行動モデル夏の学校

11長期観測データの活用• 変動のモデル化– 変動を明示的に導入して日々の行動をモデル化– 行動の基にある，より安定的な意思決定原理のモデル化• 変動の除去– 一定期間に集計した行動のモデル化– 一週間の時間配分モデル等への multiple discretecontinuous choice model や多変量頻度モデル等の適用

2007/09/20-21

第6回 行動モデル夏の学校

12

Causal impactの推定• 事前の時系列から時系列モデルを構築し，事後の時系列の予測値を作成• 事後の時系列の予測値（対照群とみなす）と事後の時系列の観測値(処置群)の偏差を(Brodersen et al., 2015;処置による変化量と捉える
Zhang & Fricker, 2021)（適用例）• コロナ禍が神戸市の商店街の客足に及ぼした影響の分析に適用（瀬谷・大森，2021）
2021/09/20-21

第20回 行動モデル夏の学校

13神戸市内の全商店街の客足：状態空間モデル• 𝑦𝑦𝑖𝑖,𝑡𝑡 = 𝜇𝜇𝑖𝑖,𝑡𝑡 + 𝑐𝑐𝑖𝑖,𝑡𝑡 + ℎ𝑖𝑖,𝑡𝑡 +𝜀𝜀𝑖𝑖,𝑡𝑡 , 𝑡𝑡 = 1, … , 𝑇𝑇

トリップ数

𝜀𝜀𝑖𝑖,𝑡𝑡 ~𝑖𝑖. 𝑖𝑖. 𝑑𝑑. 𝑁𝑁(0, 𝑉𝑉𝜀𝜀,𝑖𝑖 )𝜇𝜇𝑖𝑖,𝑡𝑡 ：トレンド（ローカルレベル）成分

𝜇𝜇𝑖𝑖,𝑡𝑡+1 = 𝜇𝜇𝑖𝑖,𝑡𝑡 + 𝜉𝜉𝑖𝑖,𝑡𝑡 , 𝜉𝜉𝑖𝑖,𝑡𝑡 ~𝑖𝑖. 𝑖𝑖. 𝑑𝑑. 𝑁𝑁(0, 𝑉𝑉𝜉𝜉,𝑖𝑖 )𝑐𝑐𝑖𝑖,𝑡𝑡 : 週周期成分（曜日成分）ℎ𝑖𝑖,𝑡𝑡 : 祝日成分：時間不変の回帰項ℎ𝑖𝑖,𝑡𝑡 = 𝒙𝒙′𝑖𝑖,𝑡𝑡 𝜷𝜷𝜀𝜀𝑖𝑖,𝑡𝑡 ：誤差成分（ガウス分布）

出典：瀬谷・大森(2021)

14兵庫県の感染者数

CASE4

事前

CASE1: 第1波

CASE2: 第2波
CASE3: 第3波

出典：瀬谷・大森(2021)
2021/09/17-19

第20回 行動モデル夏の学校

15

結果：CASE4，垂水区の商店街(メッシュID:513570543)「予測値（対照群）」と「観測値（処置群）」

予測値からの偏差偏差の累積
2021/09/17-19

第20回 行動モデル夏の学校 出典：瀬谷・大森(2021) 16

県間移動自粛と外出自粛が新型コロナの感染に及ぼす影響の分析(Liu & Yamamoto, 2021)
外出トリップ数

県間トリップ数• モバイル空間統計データ：2020/1/15-12/31，首都圏4県• 最小移動仮定での近似トリップ数(Hara and Yamaguchi, 2021)
2021/09/17-19

第20回 行動モデル夏の学校

17

DELPHI model (Li et al., 2020)
SEIRモデルの発展：無症状や自宅療養を考慮

2021/09/17-19

第20回 行動モデル夏の学校

18

外出と県間移動自粛の考慮• 県間移動による影響– 来訪者からの感染– 来訪地での感染
𝐼𝐼𝑘𝑘 𝑡𝑡 + ∑𝑔𝑔≠𝑘𝑘,𝑔𝑔∈𝑅𝑅 𝑄𝑄𝑔𝑔𝑔𝑔 𝑡𝑡 𝐼𝐼𝑔𝑔 𝑡𝑡 ⁄𝑁𝑁𝑔𝑔 𝑡𝑡
𝑑𝑑𝑆𝑆𝑘𝑘= −𝑅𝑅1,𝑘𝑘 𝑡𝑡 𝑆𝑆𝑘𝑘 𝑡𝑡
𝑁𝑁𝑘𝑘𝑑𝑑𝑑𝑑𝐼𝐼𝑔𝑔 𝑡𝑡 + ∑𝑗𝑗≠𝑔𝑔,𝑗𝑗∈𝑅𝑅 𝑄𝑄𝑗𝑗𝑗𝑗 𝑡𝑡 𝐼𝐼𝑗𝑗 𝑡𝑡 ⁄𝑁𝑁𝑗𝑗 𝑡𝑡−�𝑅𝑅1,𝑔𝑔 𝑡𝑡 𝑄𝑄𝑘𝑘𝑘𝑘 𝑡𝑡 𝑆𝑆𝑔𝑔 𝑡𝑡 ⁄𝑁𝑁𝑔𝑔𝑁𝑁𝑔𝑔𝑔𝑔≠𝑘𝑘,𝑔𝑔∈𝑅𝑅

R1,k (t ) = λk2021/09/17-19

Qk (t )
QkN (t )外出トリップの影響
第20回 行動モデル夏の学校

19

推定結果
MAPE

東京

感染者数

19.9% 16.7% 17.8% 17.7%

死者数

7.9%

千葉

神奈川

埼玉

33.1% 18.3% 19.6%

東京都

千葉県

神奈川県

埼玉県

感
染者
数

死者
数

感染状況

次頁で影響分析

2021/09/17-19

感
染者
数

死者
数
第20回 行動モデル夏の学校

20

県間移動自粛と外出自粛の効果

千葉
県

東
京都

神奈
川県埼玉県
2021/09/17-19

第20回 行動モデル夏の学校

21

付録• 次頁からは2007年第6回行動モデル夏の学校の講義資料（抜粋）の再掲です

2021/09/17-19

第20回 行動モデル夏の学校

22

動学化へのアプローチ2• 通常のモデルはある時点のサービス水準が同じ時点の交通手段選択を決定すると仮定している• サービス水準が変化したら直ぐに行動も変化するのか？

動的モデル構築の必要性
2007/09/20-21

第6回 行動モデル夏の学校

23

複数時点での状態のモデル化
Xt-1

Xt

Xt+1

時点間効果

慣性遅れ効果
先取り効果
Yt-1

Yt

Yt+1
時点T-1
2007/09/20-21

T
第6回 行動モデル夏の学校

T+1
24

複数時点の状態モデルの問題点
Yt = aYt-1 + bXt + cXt-1 + dXt+1 + εt
• 誤差項は時点間で相関を持つ（非観測異質性は時点間で同一となる）• 上記のようなモデル構造だと，Xtの値が増加した場合も減少した場合もYtに与える影響の大きさは同じで正負対象（変化の対象性）• 自動車保有台数等では，変化の対象性が成り立たないケースが多く見られる
2007/09/20-21

第6回 行動モデル夏の学校

25

複数時点での変化のモデル化
∆Xt-1

∆ Xt

∆Xt+1

∆Yt-1

∆ Yt

∆Yt+1

時点
T-1
2007/09/20-21

T
第6回 行動モデル夏の学校

T+1
26

複数時点の変化モデルの問題点∆Yt = a∆Yt-1 + b ∆Xt + c ∆Xt-1 + d ∆Xt+1 + εt
• 説明変数は ∆Xt の他に， Xt も含まれる• 変化の非対称性も考慮可能• 時間は連続なのに，離散的に扱っている• 意思決定の時点とモデル化が整合していないケースも生じる
2007/09/20-21

第6回 行動モデル夏の学校

27

連続時間軸上での変化のモデル化
生存時間モデルの適用• 連続時間を扱える• 時点間の誤差項の相関等で苦労する必要がない• 意思決定（選択）行動について，明示的には効用理論の枠組みで表現していない– 効用理論との整合性について示した論文もあり関連付けることは可能（小林ら, 1997; 佐々木,1997）
2007/09/20-21

第6回 行動モデル夏の学校

28

生存時間モデル
ハザード関数： t 以前に事象が発生していないという条件下で t に事象が発生する条件付き確率密度h(t ) = lim∆t →0Pr (t + ∆t > T ≥ t T ≥ t )∆tf (t )=1 − F (t )f(t): 確率密度関数
F(t): 累積分布関数
2007/09/20-21

第6回 行動モデル夏の学校

29

観測打ち切り• 事象が最後まで発生しないケースも解析可能– データから削除すると推定パラメータのバイアスにつながる

被験者2
被験者1
Time
観測開始時点2007/09/20-21

観測終了時点S ( t ) = Pr( T ≥ t ) = 1 − F ( t )第6回 行動モデル夏の学校

30

尤度関数

{S (t ) }( )= ∏ {h(t )S (t )}

L = ∏ f (t )(1−δ i )δii

1−δ ii

δi: 事象が観測できなかったことを表すダミー変数
以下のように，S(t) が h（t）で表せることから上記

は h（t）の関数
2007/09/20-21

S (t ) = exp − ∫ h(s )ds  s =0t

第6回 行動モデル夏の学校

31

ハザード関数の構造化• 比例ハザードモデル：基準ハザード関数に異質性が影響を与える．時点間で一定

h(t ) = h0 (t )r (θX )• 加速故障モデル：時間の進み方に異質性が影響を与える

h(t ) = h0 {r (θX )t}
2007/09/20-21

第6回 行動モデル夏の学校

32

構造化の方法によるモデルの分類
モデル

基準ハザード異質性

ノンパラメトリックモデル ノンパラメトリック

ノンパラメトリック

セミパラメトリックモデル ノンパラメトリック

パラメトリック

パラメトリックモデル

パラメトリック

2007/09/20-21

パラメトリック

第6回 行動モデル夏の学校

33

観測精度との整合性• 通常の調査では，変化の時点を連続的に観測することは稀であり，通常は1年に1度等のパネル調査が用いられる• 生存時間モデルでは「変化」がある時間内に生起したことを表現することも生起確率を t について積分することで表現可能（積分区間は任意に設定可能）
Pr (t1 < s ≤ t2 ) = ∫t2s = t1f (s )ds= F (t2 ) − F (t1 )2007/09/20-21= S (t1 ) − S (t2 )第6回 行動モデル夏の学校

34

参考文献
岩崎学（2015）統計的因果推論，朝倉書店
北村隆一（2003）変動についての試行的考察，土木計画学研究・論文集，20, 1-15.
小林潔司，喜多秀行，後藤忠博（1997）ランダム効用理論に基づく滞在時間モデルに関する理論的研究，土木学会論文集，576/IV-37, 43-54.
佐々木邦明（1997）行動論的背景を持たせた生存時間モデルによるパネル
アトリションバイアスを修正した選択モデル，土木計画学研究・講演集，20(2), 181-184.
瀬谷創・大森優希（2021）COVID-19が神戸市の商店街の滞在人口にもたらした影響の分析，「モバイルビッグデータを活用したコロナ禍における移動活動実態分析」成果報告会，2021年8月17日，オンライン
エステル・デュフロ，レイチェル・グレナスター，マイケル・クレーマー（2019）
政策評価のための因果関係の見つけ方 ランダム化比較試験入門，日本評論社
吉田寿夫（2018）本当にわかりやすいすごく大切なことが書いてあるちょっと進んだ心に関わる統計的研究法の本Ⅲ，北大路書房
吉田寿夫，村井潤一郎（2021）心理学的研究における重回帰分析の適用に関わる諸問題，心理学研究，https://doi.org/10.4992/jjpsy.92.19226
2021/09/17-19

第20回 行動モデル夏の学校

35

Brodersen, K. H., Gallusser, F., Koehler, J., Remy, N., & Scott, S. L. (2015)
Inferring causal impact using Bayesian structural time-series models. The
Annals of Applied Statistics, 9(1), 247-274.
Hara, Y., & Yamaguchi, H. (2021) Japanese travel behavior trends and change
under COVID-19 state-of-emergency declaration: Nationwide observation
by mobile phone location data. Transportation Research Interdisciplinary
Perspectives, 9, 100288
Heckman, J. (1974) Shadow prices, market wages, and labor supply.
Econometrica, 42(4), 679–694.
King, G., & Nielsen, R. (2019). Why propensity scores should not be used formatching. Political Analysis, 27(4), 435-454.
Li, M.L., Bouardi, H.T., Lami, O.S., Trikalinos, T.A., Trichakis, N.K., & Bertsimas,
D. (2020) Forecasting COVID-19 and analyzing the effect of government
interventions. medRxiv.
Liu, S., & Yamamoto, T. (2021) Role of travel restrictions and stay-at-homerequests in preventing the spread of COVID-19 in Japan. (under review)
Rosenbaum, P. R., & Rubin, D. B. (1983). The central role of the propensity
score in observational studies for causal effects. Biometrika, 70 (1), 41–55
Zhang, Y., & Fricker, J. D. (2021). Quantifying the impact of COVID-19 on nonmotorized transportation: A Bayesian structural time seriesmodel. Transport Policy, 103, 11-20.
2021/09/17-19

第20回 行動モデル夏の学校

36