【因果と矢印】DAGの初歩を解説する【統計的因果推論】

はじめに

ここ最近、統計的因果推論の文脈で特殊なモデルとしてITSを扱ってきましたが、よくよく考えてみると初歩的なことを全く書いていないので少し整理した記事を書きます。今日はDAGという、因果推論の初歩、矢印で因果関係を表現する方法を初心者向けに解説します。

目次

DAGとはなにか

有向非巡回グラフ（DAG）の定義と特徴

DAG：directed acyclic graph 有向非巡回グラフと言います。
DAGは因果推論の文脈でのみ使われるという訳ではなく、一般的なグラフ理論用語でもあります。因果推論のグラフィカルモデルはグラフ理論から派生した理論なので、例えば競技プログラミングなどで、DAGという言葉を聞いたことがある方は多いのではないでしょうか。
以下、本記事では「一般的なグラフ理論用語で使われるDAG」という文脈で「DAG」を用いる時は必ずこの旨を明記するようにし、断りなく「DAG」という用語を使うときは「因果推論の文脈でのDAG」という意味で書きます。DAGには以下のような特徴があります。

DAGの特徴と書き方

1. 頂点を矢印で結んで書く
2. DAGは巡回しない
3. 結ばれた矢印は変数同士の直接の因果関係を表す

1. 頂点を矢印で結んで書く
因果推論の文脈では、変数のことを頂点（あるいはnode）と言います。
そして一般的なグラフ理論では頂点を辺(edge)で結んでいきます。辺には有向辺と無向辺がありますが、グラフ理論において、DAGは、無向辺（線分）ではなく、有向辺（矢印）で結んでいきます。

2. DAGは循環、巡回しない
グラフ理論においてDAGは循環しないという原則があります。例えば、下図のような巡回するグラフはDAGではありません。

DAGではない有向グラフ

3. 結ばれた矢印は変数同士の直接の因果関係を表す
今までの2つの定義は一般的なグラフ理論におけるDAGと同様です。これを因果推論に拡張します。
因果推論におけるDAGでは、頂点は変数を表し、X→Yという矢印があったとき、それは変数Xから変数Yに直接の影響があることを意味します。
「風が吹けば桶屋が儲かる」という言葉がありますが、あれはたくさんの間接的な因果関係を介在しているので風が吹く→桶屋が儲かるというDAGを書いてはいけません。

因果推論における因果グラフの役割

DAGを書くことで、検証したい因果関係に対して、いろんな変数が介在しているとき、統計モデルを作る上で、考慮すべき変数と、無視すべき変数を機械的に整理することができます。

経路が開いている、閉じているとは？

このようなDAGをまず見てみましょう。このDAGを使って、XとYの直接の因果関係を見たいとします。ただ、上で見た通り、矢印で繋がっているところはそれぞれ、直接の関連がありますから、例えば、「AとX」「AとY」がこの図では直接の関連があるということになります。そうすると、X→Yという直接の経路以外にも、Aを介した、「X←A→Y」という経路がXとYの関係に影響を与えている可能性があります。ただし、直接の関連があるもの同士が繋がっているからと言って、その相方同士に影響を及ぼしているとは限りません。あなたの知り合い2人を連想したとき、その2人がそれぞれ知り合いであることもあれば、全く繋がりがないこともあるでしょう。そんな感じです（？）

DAGにおいて、経路が開いているということは、その経路で関連が存在することを意味します。反対に、閉じている経路に関してはその経路で関連が存在しないということを意味します。

因果推論においては、知りたい因果関係以外の関連をすべて無くして、知りたい因果関係があるかを調べたいというモチベーションがあります。知りたい因果関係以外の開いている経路を把握することがまずモチベーションになります。つまり、この図では、「X→Y」以外に「X←A→Y」、「X→C→Y」、「 X→B←Y」という経路があるが、まずこれらのうちどれが開いている経路なのかを特定することがやりたいことになります。

矢印の向きで関連があるかを見る

実は、DAGを利用することで、矢印の向きを見るだけで、数学的な処理をすることなく、視覚的に開いている、閉じているが判断できます。これがDAGを利用するメリットです。ただ、最初は説明があった方が良いと思うので具体例を引きながら説明していきます。

共通の原因があるとそれを介した経路は開いている

共通原因(common cause)としてよくある例は、毛髪量と年収の例です。毛髪量が少ないほど年収が高いというデータがあったとします。直感的には散髪をして毛髪量を減らしたら年収が増える、植毛をしたら年収が減る、ということはなさそうです。ということはこれらには因果関係はなさそうです。これの「妥当な」説明としては、裏に「年齢」という共通原因があり、これがそれぞれに影響を与えているというようなことができます。これをDAGで書くとこのようになります。

共通原因を介した経路、すなわち「毛髪量←年齢→年収」という経路がありますが、毛髪量と年齢は関連があるので、この経路は「開いている」という言い方をすることができます。
このような年齢のような、共通の原因をもたらす変数を「交絡因子」と言います。

共通の結果があるとそれを介した経路は閉じている

共通効果(common effect)の具体的な例として、肥満を挙げます。肥満は遺伝的要因と、運動不足などの環境要因の相乗効果で生じていると考えられます。

しかし、一般的には遺伝的要因と環境要因は関係がないので、例えば、より運動量が減ったからと言って、遺伝的要因に影響を及ぼす訳ではありません。
共通の結果を介した経路、すなわち「遺伝的要因←肥満→環境的要因」という経路がありますが、遺伝的要因と環境要因は関連がないので、この経路は「閉じている」という言い方をすることができます。
このような肥満のような、共通の結果となるような変数を「合流点(collider)」と言います。

中間因子があるとその経路は開いている

一般に、糖尿病があると死亡リスクは高くなると考えられます。糖尿病と死亡リスクの関係を知りたいです。糖尿病というのは糖の代謝異常などにより高血糖が生じるような病態ですが、高血糖だけで死ぬことは（糖尿病性ケトアシドーシスのような一部の例外疾患を除けば）死ぬことはありません。糖尿病で死亡リスクが上がるのは、糖尿病による高血糖が血管を損傷し、心筋梗塞のような心血管疾患リスクが上昇することで、死亡率の上昇につながります。
このDAGを書くと次のようになります。

このDAGでは、媒介因子である心血管疾患を介する因果経路、すなわち「糖尿病→心血管リスク→死亡率」という経路がありますが、これは糖尿病と死亡率には関連があるので、この経路は「開いている」という言い方をすることができます。
このような心血管リスクのような、因果関係を媒介するような変数を「媒介因子(mediator)」と言います。

上で見た経路3つを整理するとこのようになります。

XからYへの因果関係を見たい時、X→Y以外に「X←A→Y」、「X→C→Y」、「 X→B←Y」という経路があります。上で見たように、「開いている」経路としては「X←A→Y」、「X→C→Y」が当てはまり、「閉じている」経路としては「 X→B←Y」があることがわかります。

フロントドアパスとバックドアパス

ここで、フロントドアパスとバックドアパスという概念を導入します。フロントドア、バックドアというのは野球の変化球でも使われる概念でパスはサッカーとかバスケで使う概念でなんのスポーツなのかよくわかりませんね。
フロントドア、バックドアというのは表道、裏道という意味、パスは経路のことですから、直訳すると表道経路、裏道経路ということになりますね。
フロントドアパス、バックドアパスというのは開いている経路を区別するための概念になります。本来は厳密な定義がありますが、初歩記事なので概念のとしての導入をします。厳密な定義はバックドア基準、フロントドア基準などでググってみてください。余力があれば記事を書きます。
先ほどから多用しているDAGから開いている経路だけを抽出します。

さて、どっちがフロントドアパス（表道）っぽいでしょうか。そう、「X→C→Y」ですね。XからYへの影響を考えたいとき、この流れの向きに矢印が向いている経路を「フロントドアパス」と言い、そうでないものを表道とは反対方向に行く裏道経路、「バックドアパス」と言います。

経路は開くこともできるし閉じることもできる

「ある操作」をすると経路は開くこともできるし、閉じることもできます。経路が「開いている」というのはその経路を介した関連がある、ということでしたから、経路を閉じることができれば、その経路を介した関連を無視することができます。因果推論のモチベーションは、ある変数から、ある変数に及ぼす直接影響を知りたいということでした。ですから、この経路を閉じるという操作が、因果推論のモチベーションと親和性があることはお分かりでしょうか。余計な影響を与える影響を閉じることで、直接的な影響を知ることができます。
しかし、余計な経路を閉じてしまうと、本来知りたかった関連も見れないことになってしまうし、余計な経路を開いてしまうと余計な影響が入ってしまうことになります。この「ある操作」を適切な経路に対して行っていくことが大事です。

「ある操作」って？

この経路上にある変数に対し、統計学的に処理を行うことです。例えば、その変数で層別化することや多変量回帰分析に変数を入れることです。例えば、交絡因子を介するパスは交絡因子で層別化することにより交絡因子の影響を取り除いた解析ができます。これはすなわち、交絡因子を介した関連がなくなるということであり、つまり「交絡因子を統計学的に処理したことにより開いていたパスが閉じた」、ということになります。また、合流点で層別化をしてしまうと、全く関係がなかった2つの変数があたかも関連があるように見えることがあります（合流点バイアスと言います。知らない方は検索してみてください）。有名な選択バイアスなどもこれに該当することがあります。つまり、合流点で層別化をしてしまうことは閉じていたパスを介して、関連が生じてしまうということになります。これはすなわち、「合流点での統計学的処理によって閉じていたパスが開いた」ということに他なりません。これらの例からパス上にある変数に統計学的処理を行うと閉じていたパスは開き、開いていたパスは閉じるということがわかります。

結局どの経路を閉じたらいいの？

結論から言うとバックドアパス経路を閉じたら良いです。余談ですが、ドア、閉じるというと某映画を思い出しますね。あの映画は出現する全部の開いているドアを閉じればいいので簡単ですが、DAG上ではバックドアだけしか閉じてはならない意味でさらに難しいです。
バックドアを閉じるのは交絡因子の影響を取り除きたいと言うことを考えるとわかりやすいと思うんですが、フロントドアを閉じてはいけないのは非直感的な気がします。結局知りたいのはXからYへの影響なので、中間因子を介した影響も加味したい、と言うのが直感的な説明かなと思います。

交絡因子の古典的定義

・交絡因子は原因と関連する
・交絡因子は原因がなくとも結果に関連する
・交絡因子は原因や結果の上流にある

例えばこの例では、年齢は毛髪量と言う原因と関連し、毛髪量という要素がなくても年収と関連し、毛髪量と年収の上流にあることから交絡因子ということができます。
これだけでは対処できない交絡因子が存在する事がわかってきました。

交絡の定義を拡張

→因果推論の文脈では、交絡因子は統計学的に処理する事でバックドアパスを閉じられるかどうか？と言う解釈で捉えられています。バックドアパスを形成していて、その変数を処理することでバックドアパスを閉じられるものが交絡因子となります。

例えば、X→Yの影響を知りたい上のモデルで、CはXの上流にないので、古典的な交絡因子の定義を満たしませんが、Cは「X←A→C→Y」というバックドアパスを形成しており、Cは交絡因子ということができます。直感的には交絡因子とは捉えにくいものでもDAGを書くことで交絡因子を特定しやすくなることがDAGを使う良いところです。

モデルによって何が交絡因子なのかが違う

このDAGを見てみましょう。X→Yの因果を見たい時はAは交絡因子ですが、A→Yの因果を見たい時はXは中間因子です。A→Xの因果を見たい時はYはYは合流点です。
見たい因果関係によって構築するモデルは違うし、採用しなければいけない変数も変わってくることがわかります。

予測と因果推論は違う

因果推論においては選択的に因果を見たいので統計学的処理をする、しないなどがありますが、予測においては、よっぽどそれがノイズにならない限り変数を採用します。全然頭の使い方が違うことは少し頭の片隅に置いておいても良いでしょう。
Deep Learningの予測タスクなんかはデータ量、変数量こそ正義みたいなところもあり、いい感じの特徴量を持つ変数を追加で生み出すみたいな乱暴なこともするので、そういうものとは少し違うということは知っておいても良いと思います。

おわりに

本日は因果推論の初歩、DAGで因果関係を表現することを扱いました。今後はモデル化も扱っていきたいところです。

references

入門 統計的因果推論

Pearl流の因果推論の基礎の本です。

因果推論の科学　「なぜ？」の問いにどう答えるか (文春e-book)

これは数式を使わずにPearlが因果推論のお気持ちを読み物として書いてくれています。いい本ですが、めちゃくちゃ分厚いので、数式がないわりに気合を入れて読まなきゃいけないのがややしんどいです笑

つくりながら学ぶ! Pythonによる因果分析: 因果推論・因果探索の実践入門 Compass Booksシリーズ

初心者向けに因果推論のやりたいことを易しく書きながらPythonでの実装も書いてくれています。やや例が冗長化かつわかりにくい気もしますが・・・
Pythonが得意な機械学習、Deep Leaningの文脈で因果推論、因果探索について触れてくれていて良いです。

統計的因果推論の理論と実装 Wonderful R

Rでの実装を書いてくれています。時代の変遷や根拠となる論文を引きながら説明してくれている良い本です。細かいモデルごとの実装も説明してくれています。僕は間違えて2冊買ってしまいました。

Modern Epidemiology

英語が得意な方はこちらでしょうか。体系的に書いてくれています。なかなか他の章まで手が回っていません。

【統計的因果推論】周期性・季節性のある分割時系列解析をPythonで実装する【ITS・DeterministicProcess・Fourier級数】

1. 前回までのまとめ

少しシリーズが長くなってきたので、過去記事をどのようなモチベーションで書いてきたかを整理します。
syleir.hatenablog.com
分割時系列解析とはなにか？どのようなモチベーションでこの分析がやりたいのか？モデリングを提示し、基本的な考え方や求めたいものについてまとめました。
syleir.hatenablog.com
今までITSはRやStataでばかり実装されてきており、Pythonでの実装は確認できなかったので、Statmodelsというライブラリを用いて、折れ線回帰、セグメント回帰をPythonで実装し、グラフの描画をしました。
syleir.hatenablog.com
介入前後での傾きを求めるモデルについて、式変形を交えながら解説し、モデリングと結果の解釈について解説しました。

2. 今回の目標

時系列データには少なからず周期性が反映されます。季節性や周期性をもつ時系列データに対し、今までのような直線での折れ線回帰を用いた場合、当てはまりが悪いことがあります。そのような場合にどのような処理をしてITSを実装していくのかをPythonを使いながら説明していきます。

3. 目次

4. 使用データ

データが今回長いので折りたたんでいます。

ここをクリックで展開できます

date	y
2020-01-01	2424
2020-02-01	2596
2020-03-01	2789
2020-04-01	3770
2020-05-01	3955
2020-06-01	4203
2020-07-01	4025
2020-08-01	3933
2020-09-01	3610
2020-10-01	2889
2020-11-01	2274
2020-12-01	2376
2021-01-01	2574
2021-02-01	1966
2021-03-01	2249
2021-04-01	3290
2021-05-01	2725
2021-06-01	3543
2021-07-01	3485
2021-08-01	3183
2021-09-01	2290
2021-10-01	1959
2021-11-01	1464
2021-12-01	1716
2022-01-01	5314
2022-02-01	5106
2022-03-01	5699
2022-04-01	5700
2022-05-01	6015
2022-06-01	6363
2022-07-01	6555
2022-08-01	6143
2022-09-01	6340
2022-10-01	6229
2022-11-01	5414
2022-12-01	5496
2023-01-01	5914
2023-02-01	5736
2023-03-01	6629
2023-04-01	6510
2023-05-01	7185
2023-06-01	6843
2023-07-01	7485
2023-08-01	7073
2023-09-01	6760
2023-10-01	6919
2023-11-01	6314
2023-12-01	6306

5. 周期回帰分析について

別名をFourier analysisとか、Harmonic analysisとか言います。今回の記事の肝です。
フーリエ級数展開を利用して周期関数を回帰分析しようというモチベーションです。
周期性とフーリエ級数、フーリエ変換の相性が良いのは過去のこの辺の記事でも触れています。

フーリエ級数展開というのは
周期の関数がフーリエ級数展開できるとき、（実務上大体の関数はできます。できるための条件は専門書を参照ください。）

フーリエ解析 理工系の数学入門コース　新装版

これとかには比較的詳しく書いてあります。

と無限級数展開することを言います。
周期回帰分析では、この無限級数展開を途中で打ち切り、残りを誤差項と考え、実際のデータからを推定するものになります。
例えば、n=1で切るとすると、を推測することができれば、１次のフーリエ級数で近似した予測関数が得られることになります。何次の項まで考えるかはハイパーパラメータで、分析の時にデータを見ながら設計することになります。この係数も通常の回帰分析と同じように解析的に計算することができるので、周期関数を回帰分析で推測できるようになっています。

6. 今回のモデリング

ただしh(T)はフーリエ級数展開できる周期関数：つまり上で定義した関数とします。
今までのモデリングに加えて、h(T)が追加されています。この周期関数でのモデリングを追加することで周期性に対してrobustなモデルになります。

7. Rでのharmonic関数

Rは比較的簡単にharmonic関数という関数を用いることで簡単にフーリエ級数を用意できます。
harmonic(x, nfreq, period, intercept)というような引数を取りますが、
x:データ間隔を表す引数、月なのか日なのか
nfreq: １期あたりに含まれる周期の数（上の式でいう）
period: １期あたりのデータ点（上の式でいう）
intercept: 切片の有無(True: 上の式で、False: )
とすれば良いです。例えば年単位の周期性を考えたく、データが1月に1つの今回のデータ（12データで1年）において、2次までのフーリエ級数展開の式を用意したければ、harmonic(month,2, 12,False) とすれば良いです。

8. PythonにおけるRでのharmonic関数に相当する実装

Pythonでは関数一つでこれを解決する実装は（少なくとも自分の実力・調査範囲では）難しいです。より直接的にのデータ列を用意していくしかありません。ただ、DeterministicProcessというライブラリを使うことで少し楽ができます。これはデータ列を生成するためのライブラリで、季節ごとのダミー変数とか、月毎のダミー変数とか、今回利用するフーリエ級数列を用意することができる便利なライブラリです。
これを2次までのフーリエ級数を求める方針で実装してみると下のようになります。

#ライブラリの読み込み
from statsmodels.tsa.deterministic import DeterministicProcess, Fourier

#周期12, 2次までのフーリエ級数を用意する
fourier = Fourier(period=12, order = 2)

#データをDeterministicProcessにより生成
dp = DeterministicProcess(
    index=data.index, constant=False, order=0, drop=True, additional_terms=[fourier])

H = dp.in_sample()

pythonでのFourier関数はFourier(period, order)という引数を取ります。
period： 1周期の長さ()に相当します。Rのharmonic()におけるperiodと全く同じです。
order: フーリエ成分の数() Rのnfreqと全く同じです。
xはありませんが、前処理で調節するか、DeterministicProcessにおいてindexを目的のdatetime型に設定することで実現可能です。おそらくコードを読んでいても難しいと思うので結果を載せます。Hを最初の方だけ出力するとこうなります。

試しにH['sin(1,12)']を描画してみると次のようになります。

plt.figure(figsize=(16,8))
plt.plot(H['sin(1,12)'])
plt.xlabel('time',fontsize=18)
plt.xticks(np.arange(0, 48, 12))
plt.ylabel('y', fontsize=18)

ちゃんととなっていることがわかります。
これらのデータ列の線形結合がharmonic関数に相当するので、これらを変数に追加して今までのモデリングに適応してみれば良いです。

9. 変数の抽出

## 研究開始時点から経過した時間T
t = np.arange(0, (len(data)))
data["t"] = t

## 処置前後を表す二値変数Xt
x = np.where(data['date']>="2022-01-01", 1, 0)
data["xt"]=x

### 処置タイミングTiを取得、介入からの経過時間Xt(T-Ti)
ti = 24
data['xt(t-ti)'] = data['xt']*(data['t']-ti)

#変数の抽出
X1 = data[['t','xt','xt(t-ti)']]
H = dp.in_sample()[["sin(1,12)", "cos(1,12)","sin(2,12)","cos(2,12)"]]
y = data[['y']]

#X1とHを統合
X = X1.join(H)

これでデータの準備ができました。現在Xはこのようになっています。

モデルの作成と学習

ここからは前回記事とほぼ同じです。適宜省略します。

# 回帰モデル作成
mod = sm.OLS(y, sm.add_constant(X))
# 訓練
result = mod.fit()
# 統計サマリを表示
print(result.summary())

このように分析結果を出力することができます。原著論文では2次まで採用していたので2次で実装していますが、結果を見てみると、sin(2, 12) , cos(2, 12) は有意な変数ではないので1次までのモデリングでよかったかもしれません。

10. 推測を行う

#Xの0列目に'cep'というXと同じ長さの[1, 1, 1, ...]という列を挿入
X.insert(0, 'cep', [1]*len(X))
#学習したモデルから推測を行う
y_predict = result.predict(X)

ここで、推測してみたデータを描画してみます。

#グラフを描画
plt.figure(figsize=(16,8))
plt.title('ITS_test')
plt.plot(y, 'o')
plt.plot(y_predict)
plt.xlabel('time',fontsize=18)
plt.ylabel('event', fontsize=18)

赤実線が予測した線、青点が実際のデータですが、いい感じに学習できていることがわかります。

11. 反事実データの作成、推測

#反事実の作成
cf_X = X.copy(deep=True) 
cf_X['xt'] = 0
cf_X['xt(t-ti)'] = 0

#反事実を予測する
X_cf = cf_X[['t','xt','xt(t-ti)','sin(1,12)','cos(1,12)','sin(2,12)','cos(2,12)']]
X_cf.insert(0, 'cep', [1]*len(X_cf))
y_cf_predict = result.predict(X_cf)

#事実データに対する信頼区間を返す
predict = result.get_prediction(X)
frame = predict.summary_frame(alpha=0.05)
y_predict_l = frame['mean_ci_lower']
y_predict_u = frame['mean_ci_upper']

#反事実データに対する信頼区間を返す
predict_cf = result.get_prediction(X_cf)
frame_cf = predict_cf.summary_frame(alpha=0.05)

y_predict_cf_l = frame_cf['mean_ci_lower']
y_predict_cf_u = frame_cf['mean_ci_upper']

12. 描画

###描画
plt.style.use('fivethirtyeight')

#itv（intervention time）介入したところが何番目のデータか
#tiは上で定義済み（介入した時間）
itv = 24

plt.figure(figsize=(16,8))
plt.title('ITS_test')
plt.plot(y_predict)
plt.plot(y_cf_predict[(itv-1):], 'r--')
plt.plot(y, 'o')
plt.fill_between(data["date"], y_predict_l, y_predict_u, where= data["date"] >= "2022-01-01", color='b', alpha=0.2)
plt.fill_between(data["date"], y_predict_cf_l, y_predict_cf_u, where= data["date"] >= "2022-01-01",  color='r', alpha=0.2)
plt.xlabel('date',fontsize=18)
plt.xticks(np.arange(0, 48, 12))
plt.ylabel('event', fontsize=18)

出力がこれです。今までのモデルと比較してみましょう。

今までの折れ線回帰、セグメント回帰での出力はこれですから、相当当てはまりが改善していることがわかります。

pre-postの傾きの変化のモデルも前回の式変形の両辺にharmonic項を足すだけなので前回と同様に変数を用意して作ることができます。やってみてください。

おわりに

Rでのharmonic関数に相当する実装をPythonで行い、折れ線回帰分析に周期回帰分析を組み合わせたモデリングを行うことで周期性・季節性を考慮したITSがPythonで実装できました。ある程度は実務に耐えられるレベルまで昇華できた気がしますがいかがでしょうか。

少し難しかったでしょうか。過去の記事も見ながらじっくり取り組んでみてください。過去にはあまりフーリエ解析の記事は書けていないので基礎があやふやな方はこの辺の本を読んでみてください。統計検定準1級レベルでもたまに出るのでこの辺は押さえておくと良いです。

マンガでわかるフーリエ解析

フーリエ解析の初歩の初歩を教えてくれる良書です。漫画だからといってナメてたんですけど、初学者にはわかりやすかったのでおすすめです。

道具としてのフーリエ解析

次にフーリエ解析について学ぶならこの本がおすすめです。どう実務に生かしていくかを含めてわかります。

全実装

#ライブラリのインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.deterministic import DeterministicProcess, Fourier

#ファイルをDataFrameとして読み込む
data = pd.read_csv("dummy_data.csv")

#周期12, 2次までのフーリエ級数を用意する
fourier = Fourier(period=12, order = 2)

#データをDeterministicProcessにより生成
dp = DeterministicProcess(
    index=data.index, constant=False, order=0, drop=True, additional_terms=[fourier])

H = dp.in_sample()

## 研究開始時点から経過した時間T
t = np.arange(0, (len(data)))
data["t"] = t

## 処置前後を表す二値変数Xt
x = np.where(data['date']>="2022-01-01", 1, 0)
data["xt"]=x

### 処置タイミングTiを取得、介入からの経過時間Xt(T-Ti)
ti = 24
data['xt(t-ti)'] = data['xt']*(data['t']-ti)

#変数の抽出
X1 = data[['t','xt','xt(t-ti)']]
H = dp.in_sample()[["sin(1,12)", "cos(1,12)","sin(2,12)","cos(2,12)"]]
y = data[['y']]

#X1とHを統合
X = X1.join(H)

# 回帰モデル作成
mod = sm.OLS(y, sm.add_constant(X))
# 訓練
result = mod.fit()
# 統計サマリを表示
print(result.summary())

#Xの0列目に'cep'というXと同じ長さの[1, 1, 1, ...]という列を挿入
X.insert(0, 'cep', [1]*len(X))
#学習したモデルから推測を行う
y_predict = result.predict(X)

#反事実の作成
cf_X = X.copy(deep=True) 
cf_X['xt'] = 0
cf_X['xt(t-ti)'] = 0

#反事実を予測する
X_cf = cf_X[['t','xt','xt(t-ti)','sin(1,12)','cos(1,12)','sin(2,12)','cos(2,12)']]
X_cf.insert(0, 'cep', [1]*len(X_cf))
y_cf_predict = result.predict(X_cf)

#事実データに対する信頼区間を返す
predict = result.get_prediction(X)
frame = predict.summary_frame(alpha=0.05)
y_predict_l = frame['mean_ci_lower']
y_predict_u = frame['mean_ci_upper']

#反事実データに対する信頼区間を返す
predict_cf = result.get_prediction(X_cf)
frame_cf = predict_cf.summary_frame(alpha=0.05)

y_predict_cf_l = frame_cf['mean_ci_lower']
y_predict_cf_u = frame_cf['mean_ci_upper']

###描画
plt.style.use('fivethirtyeight')

#itv（intervention time）介入したところが何番目のデータか
#tiは上で定義済み（介入した時間）
itv = 24

plt.figure(figsize=(16,8))
plt.title('ITS_test')
plt.plot(y_predict)
plt.plot(y_cf_predict[(itv-1):], 'r--')
plt.plot(y, 'o')
plt.fill_between(data["date"], y_predict_l, y_predict_u, where= data["date"] >= "2022-01-01", color='b', alpha=0.2)
plt.fill_between(data["date"], y_predict_cf_l, y_predict_cf_u, where= data["date"] >= "2022-01-01",  color='r', alpha=0.2)
plt.xlabel('date',fontsize=18)
plt.xticks(np.arange(0, 48, 12))
plt.ylabel('event', fontsize=18)

【統計的因果推論】分割時系列解析をPythonで結果の解析・解釈をする【ITS】

1. はじめに

前回の記事ではITSをPythonで実装しながら、視覚的にわかりやすいグラフを作るところまでを目的に実装してみました。今回の記事では視覚情報だけでなく、解析的にもわかりやすいデータを出力するコードを提供し、実務に耐えるような表を実装することを目標にしています。

記事①　分割時系列解析の基本
syleir.hatenablog.com
記事②　Pythonでの分割時系列解析の実装
syleir.hatenablog.com
過去記事2作です。適宜こちらの記事に立ち戻ります。前回記事を読んでからのご参加を推奨させていただきます。

2. 求めたいものを確認する

記事①より、変数の定義と分割時系列解析によって知りたいものを確認します。

　
：切片 intercept
： 介入前の傾き pre slope
：介入後の傾き post slope
：介入前後の切片の差 level change
：介入前後の傾きの差 slope change

でした。これが求めたいものになります。論文などでITSを使う際にはこれらのパラメータの点推定値、区間推定値、(p値)などが欲しいです。
前回記事②で実装した、サマリーによって出力が得られました。

これを見てみると、ということを示唆すると捉えることができたのでした。
上述の5変数のうち、4つはこのモデルで求まりますが、はこれだけでは求まりません。
正確にはと計算することで求まりますが、これらの変数は点推定値だけではなく、95％信頼区間も欲しいです。どれくらいの精度で推定できているのかも興味がある事項になります。
信頼区間については、サマリーでこの部分を確認すれば良いことになります。

これにより、については点推定、区間推定ができているので、あとはが区間推定できれば良いことになります。

3. の推定の準備をする

これについては少しコツが必要です。
まずは、記事①で定義した式、

を変形します。
　(展開する)

( を追加し、の形に合わせる)

(の形を作る←この式変形の目的)


(変数ごとにまとめる)
とします。
式変形しただけですが、介入後の傾きが求まるモデルに様変わりしていることがわかります。
変数として採用するのが、です。元々も3変数だったので自由度も変わらず、式変形しただけなのでこのモデルで予測される折れ線回帰も前回のモデルで実装したものと同じ結果が得られます。

4. 実装していく

前回の全実装部分を実行した後にデータの前処理を追加しの項を用意していきます。

data["t(1-xt)"] = data["t"] *(1 - data["xt"])
data["t*xt"] = data["t"] * data["xt"] 

これで変数を追加しモデルを再構築します。

X = data[["t(1-xt)",'t*xt','xt']]
y = data[['y']]

# 回帰モデル作成
mod2 = sm.OLS(y, sm.add_constant(X))
# 訓練
result2 = mod2.fit() 
# 統計サマリを表示
print(result2.summary())

出力された結果が以下です。

式変形前のモデルと式変形後のモデルの学習結果を比較してみます。
（当然ですが、）偏回帰係数のパラメータ以外は全く同じものができています。

これらのパラメータを表にいい感じにまとめれば論文にも耐える形での実装が可能です。

5. おわりに

簡単ではありますが、PythonでのITSの実装例を提供しました。

小児科診療 2021年 06 月号 [雑誌]

この本で触れられている初歩的なものをさらに解説し、Pythonistaに向けた実装ができたと思います。この本では実際の論文を用いながら解説してくれているので併用しながら読んでみるとさらなる理解が深まると思います。

次は周期性のあるデータをITSで分析する記事を出そうかと考えています（いつになるかはわかりません）。
時系列記事のPart4あたりが予習に良いかと思います。ここまではただの折れ線回帰をお話していましたが、この辺から時系列ぽさが出てきます。
syleir.hatenablog.com

6. 全実装

前回実装と今回実装をまとめたものを出しておきます。

#ライブラリのインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

#ファイルをDataFrameとして読み込む
data = pd.read_csv("dummy_data.csv")

## 研究開始時点から経過した時間T
t = np.arange(0, (len(data)))
data["t"] = t

## 処置前後を表す二値変数Xt
x = np.where(data['date']>="2021-10-01", 1, 0)
data["xt"]=x

### 処置タイミングTiを取得、介入からの経過時間Xt(T-Ti)
ti = data[data['date'] == "2021-10-01"].index
data['xt(t-ti)'] = data['xt']*(data['t']-ti)

#変数の抽出
X = data[['t','xt','xt(t-ti)']]
y = data[['y']]

# 回帰モデル作成
mod = sm.OLS(y, sm.add_constant(X))
# 訓練
result = mod.fit() 
# 統計サマリを表示
print(result.summary())

#Xの0列目に'cep'というXと同じ長さの[1, 1, 1, ...]という列を挿入
X.insert(0, 'cep', [1]*len(X))
#学習したモデルから推測を行う
y_predict = result.predict(X)

#反事実の作成
cf_data = data.copy(deep=True) 
cf_data['xt'] = 0
cf_data['xt(t-ti)'] = 0

#反事実を予測する
X_cf = cf_data[['t','xt','xt(t-ti)']]
X_cf.insert(0, 'cep', [1]*len(X_cf))
y_cf_predict = result.predict(X_cf)

#事実データに対する信頼区間を返す
predict = result.get_prediction(X)
frame = predict.summary_frame(alpha=0.05)
y_predict_l = frame['mean_ci_lower']
y_predict_u = frame['mean_ci_upper']

#反事実データに対する信頼区間を返す
predict_cf = result.get_prediction(X_cf)
frame_cf = predict_cf.summary_frame(alpha=0.05)

y_predict_cf_l = frame_cf['mean_ci_lower']
y_predict_cf_u = frame_cf['mean_ci_upper']

###描画
plt.style.use('fivethirtyeight')

#itv（intervention time）介入したところが何番目のデータか
#tiは上で定義済み（介入した時間）
itv = ti[0]
#tot(total time)総観測期間
tot = len(y)

plt.figure(figsize=(16,8))
plt.title('ITS_test')
plt.plot(y_predict)
plt.plot(y_cf_predict[(itv-1):], 'r--')
plt.plot(y, 'o')
plt.fill_between(np.arange(itv,tot),y_predict_cf_l[itv:], y_predict_cf_u[itv:], color='r', alpha=0.2)
plt.fill_between(np.arange(itv,tot),y_predict_l[itv:], y_predict_u[itv:], color='b', alpha=0.2)
plt.xlabel('time',fontsize=18)
plt.ylabel('event', fontsize=18)

###slope changeを求めるモデルを作る
#データの前処理
data["t(1-xt)"] = data["t"] *(1 - data["xt"])
data["t*xt"] = data["t"] * data["xt"] 

#変数を準備する
X = data[["t(1-xt)",'t*xt','xt']]
y = data[['y']]

# 回帰モデル作成
mod２ = sm.OLS(y, sm.add_constant(X))
# 訓練
result2 = mod2.fit() 
# 統計サマリを表示
print(result2.summary())

【統計的因果推論】分割時系列解析をPythonで実装する【ITS 回帰不連続デザイン】

はじめに

さて、前々回記事で分割時系列解析、分割時系列デザイン(Interrupted Time Series Analysis : ITS, ITSA)について解説したばかり（5ヶ月前）ですね。当該記事で2週間後には書くと言った気もしますが。。。
syleir.hatenablog.com
今日はこれをPythonで実装してみたいと思います。

この図を書くのが今日の目標です。
調べてみたところ該当する記事は見当たらなかったので参考にしてください。独力で実装したので改善点、誤りなどあれば教えてください。全実装はページの一番下にあります。目次からもジャンプできます。

もくじ

0.使用データ

date	y
2020-10-01	2030
2020-11-01	2201
2020-12-01	2100
2021-01-01	2000
2021-02-01	2035
2021-03-01	2343
2021-04-01	2854
2021-05-01	2400
2021-06-01	2748
2021-07-01	2443
2021-08-01	2600
2021-09-01	3194
2021-10-01	1506
2021-11-01	1340
2021-12-01	1456
2022-01-01	1256
2022-02-01	1345
2022-03-01	1100
2022-04-01	1304
2022-05-01	1405
2022-06-01	1245
2022-07-01	852
2022-08-01	1200
2022-09-01	840

今回はこのデータを.csvで持ちます。これをpythonで読み込んでデータ処理していきます。2022年9月に公開しようと思ってたのでデータは2022年9月までです。

1.データの読み込み

import pandas as pd
#ファイルをDataFrameとして読み込む
data = pd.read_csv("dummy_data.csv")

2.データの前処理

読み込んだデータを処理していきます。まずは今回作るモデルの
を処理していきます。この辺の話は前回記事を参考にしてください。
今回のデータでは2021-10-01に処置開始のイベントが生じたとします。

は実験開始からの時間
は処置前後を表す二値変数（0なら処置前、1なら処置後）
は処置が行われた時間 を表します。

import numpy as np

## 研究開始時点から経過した時間T
t = np.arange(0, (len(data)))
data["t"] = t

## 処置前後を表す二値変数Xt
x = np.where(data['date']>="2021-10-01", 1, 0)
data["xt"]=x

### 処置タイミングTiを取得、介入からの経過時間Xt(T-Ti)
ti = data[data['date'] == "2021-10-01"].index
data['xt(t-ti)'] = data['xt']*(data['t']-ti)

ここでdataを出力してみるとこのようになります。

ほしい変数がdataframeに追加されていることがわかります。

3.変数の抽出

説明変数Xから目的変数yを学習するモデルを作ります。余談ですが、慣習的に機械学習分野では(行列は大文字、ベクトルは小文字で表すことに倣って)説明変数は大文字、目的変数は小文字で表します。

#変数の抽出
X = data[['t','xt','xt(t-ti)']]
y = data[['y']]

このように抽出するとXとyは下のようになります。

4.モデルの作成と学習

pythonにおいては回帰分析はscikit-learnなどでもできますが、今回のような詳細なパラメータがほしい場合はstatsmodelsというライブラリを使用します。pythonにおける回帰分析はめちゃくちゃ簡単に実装でき、3行で結果が表示できます。良い時代ですね。add_constantという関数はではなくのように定数項も含めて評価してねという関数になります。今回はこれが必要です。

import statsmodels.api as sm

# 回帰モデル作成
mod = sm.OLS(y, sm.add_constant(X))
# 訓練
result = mod.fit() 
# 統計サマリを表示
print(result.summary())

最終行はGoogle ColaboratoryやJupyter notebookのような環境を使っている方であれば単純に

result.summary()

でも良いです。

自分はこちらの方が体裁が整っていて好きです。

5.学習結果を見る

今回の回帰分析で見るのはココです。

ちなみに、この部分を取り出すときはこうです。

coef = np.array(result.params)

としてcoefを出力すると、[ 1948.06410256, 84.41258741, -1483.502331 , -128.06643357]と出てきます。
これは学習結果が、

においてということを示しています。

6.推測を行う

実際に、この学習した式に対してデータを入れてみたときどのような挙動をするのかを検証してみたいと思います。
これは1行で実装ができて、

#学習したモデルから推測を行う
y_predict = result.fittedvalues

これで終わりです。y_predictを出力してみると

このように出てきます。

ただ、今回はこの式に介入をしていなかったときを想定した反事実の予測もしてみたいので、別実装を用意します。

別実装

#Xの0列目に'cep'というXと同じ長さの[1, 1, 1, ...]という列を挿入
X.insert(0, 'cep', [1]*len(X))

こうすると

このようにXの0列目に大量の1が挿入されます。この操作は行列のサイズを合わせるために必要で、定数項＝1*定数項とみなしています。
便宜上やるお作法です。

この行列Xを利用して推測をする関数が以下で、実質的には上と変わらない挙動をしていますが、Xを他の行列に置換することで違うXに対しても予測ができます。

#学習したモデルから推測を行う
y_predict = result.predict(X)

この時のy_predictは以下で、上の実装での結果と一致していることがわかります（当然ですが。）

ちょっとここまでの結果を描画してみます。

plt.style.use('fivethirtyeight')

plt.figure(figsize=(16,8))
plt.title('ITS_test')
plt.xlabel('time',fontsize=18)
plt.ylabel('event', fontsize=18)
plt.plot(y, 'o')

plt.plot(y_predict)

まともに学習できていることがわかります。

7.反事実データを作成

さて、統計的因果推論のキモは事実と反事実の違いを描出することでした。
この場合の反事実というのは介入を行っていない、ということに相当します。

元の変数に立ち戻ってみると、
は実験開始からの時間
は処置前後を表す二値変数（0なら処置前、1なら処置後）
は処置が行われた時間 を表します。

でした。
介入を行わない場合、
はいついかなる時も0です。
または定義できません。（無限大という考え方もできます。）
ということで、も介入前は全て0として扱っていましたから、これも0になります。

これを実装すると、

#反事実の作成
cf_data = data.copy(deep=True) 
cf_data['xt'] = 0
cf_data['xt(t-ti)'] = 0

とすれば良いです。ちなみにcfというのは反事実的な（counter-fuctual）の略です。

8.反事実を予測する

先ほどと同じように使用する変数を取り出し、1を挿入してから予測してみます。

#反事実を予測する
X_cf = cf_data[['t','xt','xt(t-ti)']]
X_cf.insert(0, 'cep', [1]*len(X_cf))
y_cf_predict = result.predict(X_cf)

9.反事実を描画する

plt.style.use('fivethirtyeight')

plt.figure(figsize=(16,8))
plt.title('ITS_test')
plt.xlabel('time',fontsize=18)
plt.ylabel('event', fontsize=18)
plt.plot(y, 'o')

plt.plot(y_cf_predict)

介入がなかった場合に時系列データがこのような推移をするという結果がよくわかると思います。
で後ろ2項が常に0の場合に相当し、
となります。というのは前回記事で書いたように処置前の傾きの推定量ですから、
このような直線形になるわけです。

10.信頼区間を表示する

www.statsmodels.org

さて、今回の記事の肝はこれです。
Pythonの問題点は折れ線回帰（セグメント回帰）における信頼区間と予測区間の描出が難しいことにあります。

信頼区間と予測区間

結構信頼区間と予測区間は混同されがちではありますが、
信頼区間：実際の直線がこの範囲に入っている確率が1-α以上
予測区間：サンプルデータがこの範囲に入っている確率が1-α以上
というように信頼区間と予測区間は意味が違います。基本的には信頼区間の方が予測区間より狭いです。
今回欲しいのは信頼区間になりますが、予測区間についても一応実装を残しておきます。

さて、一般的に「回帰分析　信頼区間」などでググると出てくるのはwls_prediction_stdを用いた以下のようなコードではないでしょうか。

from statsmodels.sandbox.regression.predstd import wls_prediction_std
prstd, lower, upper  = wls_prediction_std(result)

これで、lower, upperを出力してみます。

出てくるのは予測区間です。これを混同した記事がかなりあるので注意してください。そして今回欲しいのは信頼区間です。これは今回のwls_prediction_stdでは出てきません。predictionなので。
tedboy.github.io
この公式ドキュメントを暇な時に読んでいただければ良いですが、これ以上の機能はないようです。困りましたね。

一般的なセグメント回帰、折れ線回帰に対して予測区間を表示するための実装としてはこのような実装がありますが、（ちなみに信頼区間もできます）

#信頼区間を抽出
from statsmodels.stats.outliers_influence import summary_table
st, data, ss2 = summary_table(result, alpha=0.05)

y_predict_l, y_predict_u = data[:, 4:6].T

plt.plot(y_predict_l, 'r--', lw=2)
plt.plot(y_predict_u, 'r--', lw=2)

#描画
plt.style.use('fivethirtyeight')

plt.figure(figsize=(16,8))
plt.title('ITS_test')
plt.plot(y_predict)
plt.plot(y_predict_cf[11:], 'r--')
plt.plot(y, 'o')


plt.xlabel('time',fontsize=18)
plt.ylabel('event', fontsize=18)

plt.plot(y_predict_l, 'g--', lw=2)
plt.plot(y_predict_u, 'g--', lw=2)

これではITSにおいて肝となる反事実のグラフにおける描画ができません。（とはいえ反事実はただの単回帰分析＋αなので実装はそこまで難しくないですが）

今回使うのは
www.statsmodels.org
これのsummary_frameを使います。ドキュメント見ても難しいと思うので実装していきます。

#事実データに対する信頼区間を返す
predict = result.get_prediction(X)
frame = predict.summary_frame(alpha=0.05)
y_predict_l = frame['mean_ci_lower']
y_predict_u = frame['mean_ci_upper']

#反事実データに対する信頼区間を返す
predict_cf = result.get_prediction(X_cf)
frame_cf = predict_cf.summary_frame(alpha=0.05)

y_predict_cf_l = frame_cf['mean_ci_lower']
y_predict_cf_u = frame_cf['mean_ci_upper']

frame というのはこのような表で、信頼区間、予測区間の上下を示してくれている表になります。mean_ci_lower/upper が信頼区間で、obs_ci_lower/upper が予測区間です。ちなみにobsというのはobservationの略で、新しくデータが「観測」されるとするならこの区間にきますよ、という予測区間の意味を反映しています。

この実装の良いところはほぼ同一の実装で事実、反事実ともに信頼区間を返してくれ、重実装にならないところです。
あとは信頼区間の間をいい感じに塗れば完成します。

塗るときはfill_betweenを使いましょう。

plt.style.use('fivethirtyeight')

#itv（intervention time）介入したところが何番目のデータか
#tiは上で定義済み（介入した時間）
itv = ti[0]
#tot(total time)総観測期間
tot = len(y)

plt.figure(figsize=(16,8))
plt.title('ITS_test')
plt.plot(y_predict)
plt.plot(y_cf_predict[(itv-1):], 'r--')
plt.plot(y, 'o')
plt.fill_between(np.arange(itv,tot),y_predict_cf_l[itv:], y_predict_cf_u[itv:], color='r', alpha=0.2)
plt.fill_between(np.arange(itv,tot),y_predict_l[itv:], y_predict_u[itv:], color='b', alpha=0.2)
plt.xlabel('time',fontsize=18)
plt.ylabel('event', fontsize=18)

これでいい感じのグラフができました。

おわりに

一般的にこういう統計的手法はRやSASなどを使いますが、Pythonでも可能です。
matplotlibやseabornがある分描画はかなりやりやすく、ライブラリの使い方が難しい程度で自由度はかなり高いです。
なんらかのものに引用する場合は一声かけるか引用を明示してくれるとありがたいです。なんらかの間違いはあると思うので。。

10000字書きました。疲れました。またいつか気が向いたらなにかを書きます。

→続編です。slope change（）を求めて結果の解析を行う記事です。
syleir.hatenablog.com

おすすめの書籍

小児科診療 2021年 06 月号 [雑誌]

おそらく現状で唯一（他にあったらぜひ教えてください）分割時系列解析に触れた成書です。ぜひ。医学雑誌ですけど統計回で普通に統計の勉強になります。

syleir.hatenablog.com
上にも貼ってますが分割時系列解析の解説記事です。

syleir.hatenablog.com
統計的因果推論の初歩の初歩の解説記事です。

syleir.hatenablog.com
統計検定準1級レベルの時系列解析をまとめています。

余談

最近流行りのChatGPTに聞いてみたら大嘘吐かれました。

全実装

説明のための余計なコードを省略したものです。上からStep by stepで実行していただくとわかりやすいかもです。

#ライブラリのインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

#ファイルをDataFrameとして読み込む
data = pd.read_csv("dummy_data.csv")

## 研究開始時点から経過した時間T
t = np.arange(0, (len(data)))
data["t"] = t

## 処置前後を表す二値変数Xt
x = np.where(data['date']>="2021-10-01", 1, 0)
data["xt"]=x

### 処置タイミングTiを取得、介入からの経過時間Xt(T-Ti)
ti = data[data['date'] == "2021-10-01"].index
data['xt(t-ti)'] = data['xt']*(data['t']-ti)

#変数の抽出
X = data[['t','xt','xt(t-ti)']]
y = data[['y']]

# 回帰モデル作成
mod = sm.OLS(y, sm.add_constant(X))
# 訓練
result = mod.fit() 
# 統計サマリを表示
print(result.summary())

#Xの0列目に'cep'というXと同じ長さの[1, 1, 1, ...]という列を挿入
X.insert(0, 'cep', [1]*len(X))
#学習したモデルから推測を行う
y_predict = result.predict(X)

#反事実の作成
cf_data = data.copy(deep=True) 
cf_data['xt'] = 0
cf_data['xt(t-ti)'] = 0

#反事実を予測する
X_cf = cf_data[['t','xt','xt(t-ti)']]
X_cf.insert(0, 'cep', [1]*len(X_cf))
y_cf_predict = result.predict(X_cf)

#事実データに対する信頼区間を返す
predict = result.get_prediction(X)
frame = predict.summary_frame(alpha=0.05)
y_predict_l = frame['mean_ci_lower']
y_predict_u = frame['mean_ci_upper']

#反事実データに対する信頼区間を返す
predict_cf = result.get_prediction(X_cf)
frame_cf = predict_cf.summary_frame(alpha=0.05)

y_predict_cf_l = frame_cf['mean_ci_lower']
y_predict_cf_u = frame_cf['mean_ci_upper']

###描画
plt.style.use('fivethirtyeight')

#itv（intervention time）介入したところが何番目のデータか
#tiは上で定義済み（介入した時間）
itv = ti[0]
#tot(total time)総観測期間
tot = len(y)

plt.figure(figsize=(16,8))
plt.title('ITS_test')
plt.plot(y_predict)
plt.plot(y_cf_predict[(itv-1):], 'r--')
plt.plot(y, 'o')
plt.fill_between(np.arange(itv,tot),y_predict_cf_l[itv:], y_predict_cf_u[itv:], color='r', alpha=0.2)
plt.fill_between(np.arange(itv,tot),y_predict_l[itv:], y_predict_u[itv:], color='b', alpha=0.2)
plt.xlabel('time',fontsize=18)
plt.ylabel('event', fontsize=18)

Google chrome で はてなブログ内のTeX の数式が乱れる方へ

はじめに

最近ようやく重い腰を上げてmacOSをVentura 13.0にアップデートしChromeを最新にしたんですが、この頃からはてなブログにおけるTeXの表示がこのように見辛くなってしまいました。

Chromeで閲覧している時に数式が乱れる

TwitterやGoogleで検索してもなかなか該当する事象が見られなかったので自分だけの再現性のないものなのかもしれないですが、解決したので備忘録がてら書いておきます。

テスト環境

M1 Macbook pro 13 inch (2020)
Apple M1, RAM 16 GB, SSD 1TB
macOS Ventura 13.0

Google Chrome Version : 107.0.5304.110（Official Build） （arm64）

対策

対処方法

結論から言うと、数式の上で右クリック（or control + クリック）してMath Settings → Math Renderer → Common HTMLと設定変更すれば良いです。

どのアップデートが原因かわからないんですが、Chrome内でのデフォルトの設定が変わってしまったのかもしれません。

この設定変更により元の美しい数式へと戻すことができます。

再現されなかった環境

手持ちの

・Windows PC＋chrome
・Windows ＋ Edge
・Linux + Firefox

・iPad or iPhone + chrome
・iPad or iPhone + safari

・M1mac + safari
・M1mac + Windows OS (dual boot) + Chrome

・M1mac + Chrome + 非はてなブログのTeX環境（note）などでは今回の事象を再現できませんでした。

arm64 の Chrome でサイト内でTeXのレンダーを行うときだけ、デフォルトの設定がどこかでズレてしまったのかもしれません。まだ同じようなことで困っている方は見つけられませんでしたが、今後も困った時の備忘録として残しておきます。

余談：使った画像の引用元です。ぜひ。

syleir.hatenablog.com

【統計的因果推論】分割時系列解析の初歩を解説する【ITS 回帰不連続デザイン】

1. はじめに

たまには統計検定に関係ない統計記事を書きたいなと思って筆（キーボードでは？）を取りました。需要が迷子ですが、今後HOTになってくると思っているのでよければぜひ。このシリーズは統計的因果推論の文脈で成り立っているので初心者の方はこの記事をご参照ください。

syleir.hatenablog.com

目次

2. 今分割時系列解析がアツい！

ITSを用いた論文数の推移

これが今回の記事のモチベーションになります。分割時系列解析を用いた論文数は年々増加傾向にあり、人口に膾炙しつつあります。その割には解説記事が（特に日本語のものは）少ないです。原稿案は2年前くらいからあったんですが、そろそろ解放しても良いかと思ったので書いている次第です。

3. 分割時系列解析をなぜ導入するのか？

統計的因果推論において、介入の評価を行いたいときは一般的にランダム化比較試験：RCTを行いますよね。一般的にはランダム化によってさまざまな恩恵が得られ、介入によってアウトカムが変化しているという因果関係を示すことができます。しかしながら対照群を取れないような以下の研究デザインについてはいかがでしょうか。
・ある時間を境に法改正があった前後での事件件数の比較
・治療に対する推奨が変わった後の治療成績の比較
このような例では、介入（法改正や推奨変更）があった世界と、介入がない世界を用意して比較するという、ランダム化比較試験のようなことができません。現実に起きてしまったことは変えられないからです。You cannot change others or the past.ってよく言いますね。それです。
そのような例ではどうしたら良いでしょうか。一つの案としてはその介入が起きていない世界を頑張って探すということが考えられます。法改正に関しては例えばアメリカの特定の州で行われたとすれば、他の州を対照にとる。推奨変更は特定の病院で行われたのであれば他の病院を対照にとるなどが考えられます。この方法の欠点はなんでしょうか？
それは、母集団の違いによって信頼度が落ちるということです。母集団が違うということは、考慮に入れている変数、あるいは入れていない変数すらも調整可能というランダム化比較試験の利点を完全に消しているということです。ある程度の考察は可能ではあるものの、因果推論に用いることができるほどの信頼度は担保できなくなってしまいます。
そこで導入される概念が分割時系列解析です。

4. 分割時系列解析とは？

ということで導入が遅くなってしまいました。この分野はまだ未発展なので英語でのInterrupted Time Series：ITSに対応する日本語表記が曖昧で表記揺れがあったりします。分割時系列デザインとか回帰不連続デザイン、分割時系列解析、中断時系列解析とかって言われてるものは全部これです。この記事では分割時系列解析と呼びます。なんとか統一して欲しいんですが。


なおGoogle的民主主義では分割時系列デザインの方が検索ヒット数が多いですが反抗することとします。理由はカッコいいからです。ただし、時系列解析とは言っても、時系列解析というよりは普通の回帰分析よりの手法なのでそこを注意すると良いと思います。あまり時系列解析要素は出てきません。

4.1 分割時系列解析の概念

分割時系列解析はあるイベントの前と後に分けて、イベントが起きる前のトレンドが続いていた場合を反事実、イベントが実際に起きたあとのトレンドを事実として下図のようにモデリングします。反事実、事実という言葉がわからない場合は冒頭のリンクの記事に立ち戻ってもらえると気持ちがわかるかと思います。

非常に直感的な概念で反事実と事実のグラフが「離れている」ほど介入（イベント）の効果があることがわかります。「離れている」というものの評価も大まかに2つあるのですが、ここでは省略し後ほど説明します。

5. どう実装するか？

反事実部分のモデリングというのは簡単で、ただの回帰分析です。
イベントが起きる前のデータを用いて単回帰すれば良いです。統計検定2級くらいの知識でしょうか。それを直線を伸ばせば良いだけなので中学2年生くらいでもできてしまうかもしれません。問題は事実部分の実装です。

5.1 折れ線回帰モデル（Segmented regression model）

見た目のままですね。このような例では折れ線回帰モデルを使います。

というモデリングをします。
ただし
は実験開始からの時間（横軸）
は処置前後を表す二値変数（0なら処置前、1なら処置後）
は処置が行われた時間 を表します。

これと実際のデータを用いて回帰分析を行い、を決定します。

5.2 パラメータの解釈

さて、このモデルにおける、はそれぞれ何を表すのでしょうか。数式をいじって検討してみます。絶対値を含む式やこのような折れ線というのは場合分けをするのが鉄板です。

5.2.1 処置前

処置前はですから、モデルの式は

となりますね。これは解釈が簡単なのではないでしょうか。
はにおける回帰直線の切片であり、は処置前の傾きです。

5.2.2 処置後

処置後はですから、モデルの式は

となります。これを変形していきます。


となります。 が処置後の傾きということになります。

5.2.3 ：なにを表す？

処置前の関数で、としてみます。つまり介入直前の値を計算してみます。
となります。
同様に、処置後の関数でとしてみます。つまり介入直後の値を計算してみます。

となります。この差がとなっていることがわかります。
つまり介入前後のにおける切片の差がになります。平易に言い換えると「介入によって値がどれだけズレたか？」になります。

5.3 パラメータの意味のまとめ

上記をまとめると、
介入前の切片
介入前の傾き
介入前後の切片の差
介入後の傾き

ということになるでしょうか。

これを図解すると下図のようになります。

パラメータの図示

5.4 level change と slope change

先ほどの章では、反事実と事実のグラフが「離れている」ほど介入（イベント）の効果があると述べました。ここまでをみて、「離れている」の定量評価を考えてみましょう。
介入の効果とは一体なんでしょうか。介入と言えば先日の2022年10月、日銀が為替介入を行いました。自分は為替とかそういうのはサッパリわかりませんが、その結果を見てみましょう。

2022/10/24 為替介入チャート

このチャートからわかるのは
・介入前は上昇トレンドを形成していた。
・介入により即時的な効果はかなり強そう。
・介入後はむしろ介入前より強く上昇している。
ということかなと思います。為替介入が円高方向（グラフの下の方向）に持っていきたい操作であるとすればこの介入は成功なのでしょうか？

少し考えてみると、介入前後でのトレンドの変化とトレンドの即時効果という2点で介入の効果を捉えてみることができそうです。

この図を見てみましょう。
介入前後のトレンドの変化というのは緑の矢印（1本目と3本目）の傾きの変化で表されます。これを、分割時系列解析ではslope changeと言います。slope(スロープ)とは傾きのことです。そのままの意味でわかりやすいですね。この介入では介入したにもかかわらず、上昇の勢いが続いていること、さらに上昇の勢いが強くなっていることがわかります。

続いて、介入の即時的な効果について考えます。
介入の即時的効果というのは青の矢印（2本目）の矢印の大きさで表されます。これを、分割時系列解析ではlevel changeと言います。level(レベル)は水平、水準といった意味です。RPGにおけるレベルは水準から転じて階級と言った意味で使われるでしょうか。ここでいうlevel changeは水平面の変化、水準変化くらいの認識で考えれば良いです。この介入は（この図からはその結論を導くのは尚早ですが）かなりの効果がありそうと見てとれます。

ということで、為替相場については素人ですが、「介入によってかなりの円高に向けた即時効果があるが、その前後では円安への上昇傾向は加速している」というように介入の効果をまとめることができます。

5.5 level change と slope changeの定量評価

さて、これを先ほどのパラメータで考えると、介入前の傾きが、介入後の傾きが であったので、slope changeは、介入前後の切片の差が であったから、level changeは と表すことができます。極論、分割時系列解析では、 と を求めるものと言うこともできますね（これは本当に極論ではあるのですが）。

6.まとめ

統計的因果推論の文脈では、ランダム化比較試験などのデザインがよく用いられるが、過去の自分のデータを対照群とする研究デザインもあり、そのひとつが分割時系列解析（ITS）である。ITSにおいては介入前のトレンドが持続するものとして反事実を作り、実際のデータと離れ具合を比較する。その指標がslope change と level change であり、統計的に定量評価が可能になっている。

7. おわりに

長くなってきてしまったので導入はこの辺で終わりにします。次の記事では実際に実装しながら概念の理解を深めていきたいと思います。なお、適用にあたっては細かい前提や実装にあたっての過分散や周期性などの注意点を大きく省略したので各自で勉強されてください。分割時系列解析についてはこの本がおすすめです。

統計的因果推論の理論と実装 Wonderful R

この分野は論文をつまみ食いしながら勉強するしかなかったんですが、やっと最近になって分割時系列解析も成書で扱われるようになってきました。この本は歴史的な背景や論文を挙げながら詳細に解説してくれているのでオススメです。ぜひこの記事を読んで興味が出てきた方は勉強してみてください。Rが読める方はオススメです。読めない方は次の記事を楽しみに待っててください（2週間くらいで頑張ります）。

次の記事です。これをPythonで実装しました。
【統計的因果推論】分割時系列解析をPythonで実装する【ITS 回帰不連続デザイン】 - Syleir’s note

【統計検定準1級】MCMCの初歩のお気持ちをできるだけわかりやすく解説する

MCMC名前はかっこいいですよね。使いこなせるとカッコいいんですが、なかなか初見で理解するのが難しいです。以前自分で勉強したときも段階的に理解できるような記事がなくて苦労した記憶があるので、今回はそのような記事を作るというのが目的です。

本日の内容です。

MCMCとは？

Markov Chain Monte Carlo法（マルコフ連鎖モンテカルロ法） のことです。一般的な記事ではマルコフ連鎖とモンテカルロ法という風に分解され〜から入ってマルコフ連鎖とモンテカルロ法のそれぞれの解説が始まる記事が多いのですが、あまり本質ではないので今回は一旦置いておきます。

そもそもMCMCで何がしたいの？

ベイズモデリングにおいて、事後分布の性質について知りたい時に使う方法です。ベイズモデリングのお気持ちをざっくり話すと、「適当な」事前分布を用意し、情報が与えられることで情報を更新し、事後分布を生成し、この事後分布を使って意思決定や分析、検定などを行う、ということを行います。つまり、事後分布についての情報ができるだけ欲しいということになります。
事後分布の生成にはベイズの定理を使います。

ベイズの定理と共役事前分布

一般的には事後分布を知りたいときにはベイズの定理を用います。
を事前分布とすると、データ が与えられた時の事後分布 は

の形で与えられます。細かいところはいつか記事になります。今理解しておくべきことは、このベイズの定理によって事前分布から事後分布が出ること、この計算は分母の積分が難しく、適切にモデリングしなければ解くことができない、ということの2点です。
一般的に事前分布に共役事前分布と呼ばれる特殊な形を使うと解析的に解けますが一般的には難しいです。一般的に積分が難しかったり、一般的な漸化式が解けなかったりとかと同じようなものです（実際ほぼ同じ）

事後分布が解析不能な形な時の方法について

事後分布が解析不能な時は、サンプリングしてなんとか情報が取れないかを検討します。これがMCMCの大事な発想につながっていきます。

サンプリングって何？

標本抽出を行うことです。標本抽出とは、母集団がどのような分布をしているかわからない時に一部の集団をとってくることです。例えば選挙とかの例が有名ですね。選挙に行ったうち一部の人にアンケートを実施して選挙行動の動向を聞くことで人口全体の選挙行動を推測して事前の趨勢を把握したり、開票が始まったタイミングで一部しかまだ開票が進んでいない段階でも全体の結果を推測して当確を出したりということができます。

今回は説明のために母集団が分かっている分布を用いて色々試してみます。

サイコロの出目のサンプリング

・例えば、サイコロを初めてみた人がいて、サイコロが、1/6ずつ等しく同じ目が出ると分かっていなかった場合、10000回サイコロを振るという試行をしてみます。
実際に振って結果をメモするのはめんどくさいのでPythonを使って試行をシミュレーションしてみます。
結果は以下の通りです。

60000回サイコロを振った結果

サイコロが等しく出目が出る一様分布に従うということを知らなかったとしてもこの結果を見るとなんとなく「出目が等しく出そう」という結論を導くことができるのではないでしょうか？

正規分布からのサンプリング

サイコロの出目のような1, 2, 3, 4, 5, 6 といった離散でなくとも、連続分布からもサンプリングができます。
ただし連続分布において1.49359203といった特定の値を取る確率は0なので、サンプリングの結果は1以上2未満に入ったなどと「幅」を持たせて集計します。今回のサンプリングは標準正規分布からのサンプリングを100万回行い-4から4までの区間を100分割して、それぞれの区間に入った個数を集計します。

標準正規分布からの100万回サンプリング

これを元々の母集団である標準正規分布のグラフと比較すると次のようになります。

標準正規分布の確率密度関数との比較

サンプリングによりかなりの確度で元の分布の性質を反映できていることがわかるのではないでしょうか。

現在の情報を使ったサンプリング

2状態、確率1/2の状態遷移

今度は、下図のような状態遷移図に従う、1秒ごとに確率1/2で状態を遷移し、確率1/2でその状態に留まるモデルを考えます。初期状態では状態1にいるものとします。これも理論上解くことができ、高校で習う確率漸化式の要領でn秒後の確率分布を求めることができますし、十分時間が経ったあとの状態1,2それぞれにいる確率は対称性から1/2と求めることができます。

さて、十分賢くない中学3年生に戻り、確率漸化式のことを学んでいない状態に戻ったとして、この十分経ったあと状態Aにいる確率がどれほどかを知りたくなったとします。
そうすると、コイントスをしながら実際にこの移動を10000回シミュレーションして、実際に何回Aにいたかを見てみることでなんとなく確率がわかることが想像できるでしょうか。コインを10000回トスしてメモしたくないので上記の設定でシミュレーションを行い、Pythonくんに結果だけ出してもらいます。

試しに最初の10回の遷移を書き出すと
['1', '2', '2', '1', '2', '1', '2', '1', '1', '2']
というような感じになります。

10000回のシミュレーションは以下のようになります。

10000回試行

ほぼ1/2ずつになっていますね。
このように時系列データに対しても長時間シミュレーションしてサンプリングすることにより収束した際の確率分布の情報を取ることができます。

MCMCのアルゴリズムを見てみる

ついに、本題に入ります。MCMCのアルゴリズムは次のようになっています。ここではMCMCのアルゴリズムのうち最も初歩的なMH（メトロポリス・ヘイスティングス法）について扱います。

①初期値を適当に設定する（定義域内ならなんでもいい）
②ループ回数を決める(多ければなんでもいい)
③遷移確率を適当に定義する（計算しやすければなんでもいい）（この状態移動に関する確率分布を提案分布と言います）

④以下を設定した回数ループする
❶現在位置を確認し、②で決めた遷移状態に従って次に移動する予定値を決める
❷実際にその予定値に移動するか、現在値に止まるかどうかを事前分布 と適当に決めた遷移ルールで決まる値を用いて確率的に決める
❸移動したか移動していないかによらず、今いるところをメモし、サンプリングする

ご覧のようにMCMCというのは①〜③まではめちゃめちゃ適当です。MCMCの本質と言っても過言ではない一番重要なエッセンスは❷です。❷がなければ先ほどの2状態の遷移となんら変わりません。「状態遷移をして、その現在位置をサンプリングする」というのが先ほどの例ですが、MCMCは「確率的に状態遷移をしてその現在位置をサンプリングする」というのに過ぎません。

MCMCで事後分布の情報を取れる理由

MCMCでは❷によって、「事後分布において事後確率が高いところにいるとき、遷移を保留し、事後確率が低いところにいるとき、遷移を積極的に行うように移動する」というようになっています。

これにより事後確率が高いところで頻回のサンプリングができるように、事後確率が低いところではあまりサンプリングが行われないようになります。

これはMCMCでは詳細釣り合いという条件として付加されているのでその辺を参照ください。そしてこのアルゴリズムはループ回数を極限まで増やすことでサンプリングの結果と事後分布が一致することが知られています（ただし収束までのループ回数は保証されていませんが。。。）

実際の挙動を見てみる

事前分布が正規分布、正規分布から得られるデータが得られる時、これらは共役事前分布であり、事後分布は正規分布になることが解析的に求められます。
ここでMCMCを実際に実装してみて挙動を比較してみます。

正規分布のMCMC

折線が移動経路の挙動、右側のグラフが集計した事後分布のサンプリングになります。実際の挙動を見てみると、上に外れても、下に外れてもすぐに真ん中に帰ってきて、真ん中にいる間はしばらく留まっているということがわかります。

マルコフ連鎖モンテカルロ法って結局なんなの？

ここまでMCMCについて、一切の言葉の説明をせずに突っ走ってきました。マルコフ連鎖モンテカルロ法はマルコフ連鎖とモンテカルロ法に分割され、それを合成した言葉なわけで、実際ある意味MCMCを表しているということなのですが、それから解説を始めるとよくわからなくなってしまうということが往々にしてあるんだと思います。MCMCのわからなさって結局ここにあったんだと思います。本稿の最後に後方視的に意味を解説していきたいと思います。

マルコフ連鎖って？

マルコフ連鎖とは、現在の状況によってのみ次の移動条件が決まる確率過程のことです。現在の状態さえ分かれば、未来の状況が予測できる、過去の状態は未来を予測するのに不要ということをマルコフ性と言います。マルコフ性に従う確率変数の列がマルコフ連鎖です。
例えば、2項間漸化式、あれがマルコフ連鎖です。先ほどの1/2の確率で状態を移動する過程もマルコフ連鎖と言えます。MCMCにおいてはシミュレーションにおいて次の移動先を決めるのに現在の位置しか関係ないのでマルコフ性があるということになります。先ほどの「事後分布において事後確率が高いところにいるとき、遷移を保留し、事後確率が低いところにいるとき、遷移を積極的に行うように移動する」とかも現在の位置のみを参照していて、過去の位置については見ていないというのが大事なポイントでしたよね。遷移するかどうかもそうですが、具体的な遷移先の座標についても現在の位置のみを参照して決まっています。

モンテカルロ法って？

モンテカルロ法とは、確率的な試行（乱数を用いた試行）によってシミュレーションを行い推定を行う方法のことです。有名なものだとビュフォンの針などでしょうか。均一幅で引かれた直線群に対して針を投げると円周率が求まると言ったものです。針を投げるという確率的操作で円周率の推定ができることがモンテカルロ法に相当します。今回の説明で言えば❷における、移動するかしないかが確率的に決まることによって事後分布が推定される、というところにモンテカルロ法が使用されています。

マルコフ連鎖モンテカルロ法って？

結局、ベイズ推定において「現在の情報のみを参照しながら遷移先を決定し」、「確率的に遷移するかしないかを決める」ことで事後分布をシミュレーションする、という方法になります。それ以上の発展的アルゴリズムや、具体的な数式の内容については成書をご参照ください。

さらなる学習のために

いくつかベイズ統計のおすすめの本を置いた記事を書いておきます。ご自身のレベルに合わせてチョイスください。
syleir.hatenablog.com