OLSをbootstrapでやるためのコードを書いてみた（Python）

Pythonで回帰分析（OLS）を行うときに、bootstrap法を指定するコードが見つからなかったので、「OLSをbootstrapで行うための関数」bootstrapOLSを自分で作ってみました。

#バージョン
python 3.8.10
statsmodels 0.12.2
scipy 1.6.2
tqdm 4.59.0

#簡単な説明
###やり方
この記事で紹介するbootstrap法では、サブサンプルをオリジナルのサンプルから復元抽出（母集団から取り出した標をは母集団に戻してからつぎの標本を取りだすやり方）でK回取り出します。
そして、各サブサンプルについてOLSで各変数の係数を推定して、その結果を1000回分箱に入れます。
最後に、それぞれの変数の係数について平均値と不偏標準偏差（解釈としては標準誤差）を計算して、result箱に入れます。
###R2
ここでは、オリジナルのサンプルのadjR2をresult箱に入れておきました。おそらく、STATAと同じのはずです。
###t値
p値ではなくt値をresult箱に入れたので、統計的な有意水準を判定するために、dataframeのサンプルサイズと自由度を使って調整したt値をresult箱に入れておきました。

#引数の説明
K: 復元抽出を繰り返す回数
→1000回くらいやれば良いと思います。5000,10000回でも試してみましたが、私の使ったデータ（サンプルサイズ4000）では、特に結果が変わらなくなりました。

sapmplesize: 復元抽出で構成するサブサンプルのサイズ
→復元抽出でサブサンプルをK回構成してOLSで係数を推定します。もしかすると、オリジナルのdataframeと違うサイズを指定する方法もあるのかな？と思ったので、蛇足かもしれませんが、サンプルサイズは指定することにしました。

Formula: 回帰式
→statsmodels.formula.apiを内部で使ってOLS推定してるので、
Formula = "y~x1+x2+C(year)"
みたいな感じで回帰式を指定してください。

df: dataframe

#コード

#　モジュール
import pandas as pd
import numpy as np
from tqdm import tqdm_notebook as tqdm
import statsmodels.formula.api as smf
from scipy.stats import t

def bootstrapOLS(K, samplesize, Formula, df):
    for r in tqdm(range(K)):
        #　dtから復元抽出（r回繰り返して変数を定義しなおさないと、毎回同じ標本が抽出される）
        ds = df.sample(n=samplesize, random_state=r, replace=True)

        #　回帰
        reg = smf.ols(formula=Formula, data=ds)
        res = reg.fit()

        #　係数を箱に入れる
        if r==0:
            #　係数を入れる箱
            coef = pd.DataFrame(res.params).T   #　res.paramsは辞書型っぽいので、そのままdataframeにして転置
        else:
            coef = pd.concat([coef, pd.DataFrame(res.params).T])

    #　係数の平均
    result = pd.DataFrame(coef.mean()).T

    #　係数の不偏標準偏差（解釈としては標準誤差）
    bootde = coef - coef.mean()   #　偏差
    bootcov = bootde.T.dot(bootde) / (K-1)   #　共分散行列（対角成分は不偏分散）
    bootse = pd.DataFrame(np.sqrt(np.diag(bootcov)), index=coef.columns).T   #　対角成分のルート（標準偏差）
    result = pd.concat([result, bootse])

    #　平均と標準誤差の列名をつける
    result = result.T.set_axis(["mean", "se"], axis="columns")

    #　t値
    result["t"] = result["mean"] / result["se"]
    result = result.reset_index()

    #　adjR2とかその他諸々を（オリジナルのdataframeから計算）
    reg = smf.ols(formula=Formula, data=df)
    res = reg.fit()
    #　係数を箱に入れる
    result = result.append({"index":"adjR2(original)", "mean":res.rsquared_adj, "se":np.nan, "t":np.nan}, ignore_index=True)
    result = result.append({"index":"N", "mean":res.nobs, "se":np.nan, "t":np.nan}, ignore_index=True)
    result = result.append({"index":"t_99", "mean":t.ppf(0.995, len(df)-(len(result)-3)-1), "se":np.nan, "t":np.nan}, ignore_index=True)
    result = result.append({"index":"t_95", "mean":t.ppf(0.975, len(df)-(len(result)-3)-1), "se":np.nan, "t":np.nan}, ignore_index=True)
    result = result.append({"index":"t_90", "mean":t.ppf(0.95, len(df)-(len(result)-3)-1), "se":np.nan, "t":np.nan}, ignore_index=True)
    
    #　Indexを設定
    result = result.set_index("index")
    
    return result
    

#結果の出力
結果はdataframeとして帰ってきます。
普通に回帰分析するみたいに、res = bootstrapOLS(~~)　とすれば、resの中に結果のdataframeが格納されます。
なお、1000回OLSを繰り返すので、それなりに時間がかかります。

手元の適当なデータで「bootstrapOLS(1000, len(df), Formula, df)」を回した結果、↓こんな感じのdataframeで帰ってきました。

#bootstrap法は何をやってるのか？
※以下の内容は、僕の指導教員の先生から教えてもらったことを僕なりにかみ砕いたものです。

bootstrap法では、復元抽出でサブサンプルを作りまくって係数を推定するので、実質的にはユニークな観測を使って係数の推定を行っていることになります。
なお、観測の重複部分は係数の推定に影響を与えないので、言い換えると、ブートストラップは復元抽出によりランダムに一部の観測をサンプルから落とすことによって、係数の分布をサンプルを変えまくって計測していることになります。

#不均一分散に対応するためにbootstrap法を使うことついて
※以下の内容は、僕の指導教員の先生から教えてもらったことを僕なりにかみ砕いたものです。

一般的に、普通にOLSで係数の標準誤差を推定すると、均一分散の仮定に基づいた推計方法が採用されます。
そうではないやり方を採用したいなら、ロバスト標準誤差（僕の解釈：残差の共分散行列を標準誤差を推定する式に挟み込む）やクラスターロバスト標準誤差（できればこれを使うのが良さそう）を使えばいいですが、サンプルサイズやクラスターの数が小さい・少ない時は標準誤差の推定方法としていまいち。

一方、bootstrapは普通にOLSを行う際に必要な仮定を置かない代わりに，手元にあるサンプルから集計量の分布をサンプルを変えまくって導出します。
したがって、ブートストラップは誤差項の不均一分散に対応するための手法というわけではなく、誤差項の分布に対する仮定全般をなくした状態で係数の標準誤差を推測していることになります。
その結果が、不均一分散時に誤って均一分散を仮定して計算した標準誤差よりも正確なものになるかどうかは、手元のサンプルが母集団をどの程度正確に反映しているか、また、誤差項に対する他の仮定はどの程度成立しているのかに依存します。