[統計] GLMでの区間推定
はじめに
前回の記事では、正規分布の平均を区間推定する方法を説明した。
同様に、一般化線形モデル(GLM)でもパラメータを区間推定することがあるが、この区間推定をどうやっているのか気になったので調べてみた。
GLM
パラメータを 個持つGLMを考える。
正確に書けば次のようになる
ここで はリンク関数、 は説明変数の行列、 はパラメータのベクトルである。 番目の応答変数 の期待値は、リンク関数を介して、説明変数とパラメータの線形結合によって説明される。
GLMのパラメータの最尤推定量が従う分布
GLMのパラメータ の最尤推定量 は、近似的に正規分布に従うことが分かっている (Dobson, 2008, p.87)。
ここで は、フィッシャー情報行列 の逆行列である。フィッシャー情報行列は、対数尤度から計算される量であり、対数尤度関数が最尤推定量の周辺でどれだけ凸型になっているかを表している(詳しくはフィッシャー情報量 - Wikipedia)。
したがって、最尤推定量は、真の値を平均とし、フィッシャー情報行列の逆行列を分散共分散行列として持つ正規分布に従うということである。上の式を推定量と真の値の差で書き直せば、次のようになる。
この式から、最尤推定量が真の値とどれだけ離れている可能性があるかを評価することができるので、信頼区間を求めることができる。
例:パラメータ2つのロジスティック回帰
例として次のデータのようなデータ( )を考える。
x | y |
---|---|
2.614672698355135 | 1 |
0.2311712788458009 | 0 |
4.130713626870605 | 1 |
3.4980589040126513 | 0 |
1.2242786762488267 | 1 |
0.5373907435482667 | 0 |
1.8970854531195753 | 1 |
0.3183052275680498 | 0 |
3.9466709572816496 | 0 |
3.5309785781654632 | 1 |
3.3761561081893414 | 0 |
0.720001443910317 | 0 |
2.064750554734893 | 0 |
2.009122195532624 | 0 |
2.2860103892990615 | 0 |
4.7927043238872375 | 1 |
1.1305070181779313 | 0 |
4.027287985243035 | 0 |
2.3250037833200765 | 0 |
4.554878548309989 | 0 |
0.6777507705415364 | 0 |
1.0177827317855659 | 0 |
4.042219354035978 | 1 |
1.5374167164460428 | 1 |
0.419902823072304 | 0 |
3.2341776629722996 | 0 |
1.2807428825075295 | 0 |
1.1185882637563382 | 0 |
0.37086337116461865 | 1 |
0.6935115140150583 | 1 |
4.775067677568351 | 1 |
3.5673200749785567 | 1 |
3.4595515158432804 | 1 |
3.195628563216748 | 0 |
0.3674178782304305 | 1 |
0.4163850367959582 | 0 |
0.23154185542532957 | 0 |
0.11134517325856519 | 0 |
4.287761228825091 | 1 |
2.2963424298511153 | 0 |
パラメータを2つ持つロジスティック回帰を考える。
番目のサンプル が1となる確率 は、切片 とパラメータ の係数 によって決まるというモデルである。
ちなみにデータを生成するときに使ったパラメータ(真のパラメータ)は である。
ではPythonのstatsmodelsライブラリを使って、ロジスティック回帰を行ってみよう。0の数もカウントする必要があるので、0の数をカウントした列を追加して、ロジスティック回帰を行う。
import pandas as pd import statsmodels.api as sm import statsmodels.formula.api as smf df = pd.read_csv('test.csv') df['z'] = 1 - df['y'] glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial()) fit = glm.fit() print(fit.summary())
先に示したデータに対してロジスティック回帰を行い、print(fit.summary())
で出力される結果を以下に示す。
Generalized Linear Model Regression Results ============================================================================== Dep. Variable: ['y', 'z'] No. Observations: 40 Model: GLM Df Residuals: 38 Model Family: Binomial Df Model: 1 Link Function: logit Scale: 1.0000 Method: IRLS Log-Likelihood: -24.577 Date: Sun, 18 Oct 2020 Deviance: 49.153 Time: 01:57:24 Pearson chi2: 40.2 No. Iterations: 4 Covariance Type: nonrobust ============================================================================== coef std err z P>|z| [0.025 0.975] ------------------------------------------------------------------------------ Intercept -1.4937 0.648 -2.306 0.021 -2.763 -0.224 x 0.4327 0.231 1.871 0.061 -0.021 0.886 ==============================================================================
サマリーの右下に2つのパラメータの信頼区間が表示されている。今回はせっかくなので、この値をフィッシャー情報行列から求めてみよう。
statsmodelsでは、GLMオブジェクトのinformation
メソッドを使うことで、フィッシャー情報行列を得ることができる(おそらく符号が逆になっていると思われたので、ソースコードの中では符号を反転させている)。
以下は の区間推定を行う例。フィッシャー情報行列の逆行列を求め、その 成分( の分散)から標準偏差を求める。±標準偏差×1.96の範囲に95%が含まれるので、ここから信頼区間を得ることができる。
import numpy as np import pandas as pd import statsmodels.api as sm import statsmodels.formula.api as smf df = pd.read_csv('GLM_test/test.csv') df['z'] = 1 - df['y'] glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial()) fit = glm.fit() x_est = fit.params['x'] # Maximum likelihood estimation of beta1 fim = - glm.information(params=fit.params) # Fisher Information Matrix (FIM) inv_fim = np.linalg.inv(fim) # Inverse of FIM std1 = np.sqrt(inv_fim[1, 1]) # Standard error of beta1 lower = x_est - 1.96 * std1 upper = x_est + 1.96 * std1 print('[{:.3f}, {:.3f}]'.format(lower, upper))
出力結果は
[-0.021, 0.886]
となって、先にサマリーで確認した信頼区間と一致する。
信頼区間の検証
GLMのパラメータの信頼区間の求め方が分かったところで、95%信頼区間が本当に95%信頼区間になっているか、データの生成と信頼区間の推定を繰り返して検証してみる。
先程と同じく のデータにロジスティック回帰を行って の信頼区間を求める。これを 回繰り返して、信頼区間の的中率を計算する。
import pandas as pd import numpy as np import statsmodels.api as sm import statsmodels.formula.api as smf import matplotlib.pyplot as plt from tqdm import tqdm # Parameters of a population beta0 = -2.0 beta1 = 0.8 # The number of samples drawn from the population at once N = 40 # The number of resample and estimation M = 1000 # The confidence coefficient is (1 - alpha) alpha = 0.05 def logistic(x): return 1 / (1 + np.exp(-x)) def make_data(n): x_arr = np.random.uniform(low=0, high=5, size=n) err_arr = np.random.normal(loc=0, scale=0.3, size=n) probs = logistic(beta0 + beta1 * x_arr + err_arr) y_arr = np.array([np.random.binomial(n=1, p=p) for p in probs]) df = pd.DataFrame({'x': x_arr, 'y': y_arr}) return df def calc_confint(df): df['z'] = 1 - df['y'] glm = smf.glm(formula='y + z ~ x', data=df, family=sm.families.Binomial()) fit = glm.fit() x_est = fit.params['x'] # Maximum likelihood estimation of beta1 fim = - glm.information(params=fit.params) # Fisher Information Matrix (FIM) inv_fim = np.linalg.inv(fim) # Inverse of FIM std1 = np.sqrt(inv_fim[1, 1]) # Standard error of beta1 lower = x_est - 1.96 * std1 upper = x_est + 1.96 * std1 print(fit.summary()) return (lower, upper) def estimation(): df = [] for i in tqdm(list(range(M))): samples = make_data(n=N) confint = calc_confint(samples) df.append(confint) df = pd.DataFrame(df, columns=['lower', 'upper']) df.to_csv('estimation.csv', sep=',') return def n_hit(): df = pd.read_csv('estimation.csv') n_hit = 0 n_miss = 0 for lower, upper in zip(df['lower'], df['upper']): if (lower <= beta1) and (beta1 <= upper): n_hit += 1 else: n_miss += 1 return n_hit def plot(): df = pd.read_csv('estimation.csv') fig, ax = plt.subplots() for i, (lower, upper) in enumerate(zip(df['lower'], df['upper'])): ax.vlines(i + 1, ymin=lower, ymax=upper, color='C0') ax.hlines(beta1, xmin=0, xmax=M + 1, color='red') ax.text(x=0.95, y=0.95, s='{:d}/{:d}'.format(n_hit(), len(df)), transform=ax.transAxes, ha='right') fig.savefig('estimation.jpg', dpi=300) return if __name__ == '__main__': estimation() plot()
結果の一例を示す。
この例では、的中率は96.5%であり、おおよそ95%になっている。
正規分布近似とWald統計量
以上の方法では、GLMのパラメータの最尤推定量は真の値を平均とする正規分布に従うという法則を用いていた。しかし、これは近似であるということに注意が必要である。データの数が少ないときには近似の精度が悪くなるので、正規分布に近似する代わりに、パラメトリックブートストラップを行う方が良いかもしれない。
最尤推定量と真の値の差が正規分布に従うのであれば、その二乗和はカイ二乗分布に従う(パラメータは独立でないので、分散共分散行列の逆行列を挟む必要がある。詳しくはDobson(2008, p.87)など)。この考え方で求まるのが、Wald統計量であり、このWald統計量を用いてWald検定が行われる。したがって、Wald検定は上で説明した信頼区間と類似した考え方だろうと思うのだが、詳細を知らないので確実なことは言えない。信頼区間の場合は分布の両側から2.5%を削って求めたが、Wald検定の場合はパラメータ が0から離れているかがポイントになるので、0に近い方から分布を5%削って有意性を調べるのだろうか。この辺りは久保(2012, p.51)に簡単な解説がある。とは言え、多くの場合はWald検定ではなく尤度比検定を使うだろうと思うので、この辺りの詳細はあまり気にする必要はないかもしれない。