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法を使うことついて