[統計] GLMでの区間推定

はじめに

前回の記事では、正規分布の平均を区間推定する方法を説明した。

cyanatlas.hatenablog.com

同様に、一般化線形モデル(GLM)でもパラメータを区間推定することがあるが、この区間推定をどうやっているのか気になったので調べてみた。

GLM

パラメータを p 個持つGLMを考える。

 Y_i \sim \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p

正確に書けば次のようになる

 g(E[Y_i]) = \boldsymbol{x}^{T} \boldsymbol{\beta}
 \boldsymbol{x}^{T} = (x_{i1}, x_{i2}, \dots, x_{ip})
 \boldsymbol{\beta}^{T} = (\beta_1, \beta_2, \dots, \beta_p)

ここで g はリンク関数、\boldsymbol{x} は説明変数の行列、\boldsymbol{\beta} はパラメータのベクトルである。i 番目の応答変数  Y_i の期待値は、リンク関数を介して、説明変数とパラメータの線形結合によって説明される。

GLMのパラメータの最尤推定量が従う分布

GLMのパラメータ \boldsymbol{\beta}最尤推定\boldsymbol{b} は、近似的に正規分布に従うことが分かっている (Dobson, 2008, p.87)。

 \boldsymbol{b} \sim N(\boldsymbol{\beta}, \mathcal{I}^{-1})

ここで  \mathcal{I}^{-1} は、フィッシャー情報行列 \mathcal{I}逆行列である。フィッシャー情報行列は、対数尤度から計算される量であり、対数尤度関数が最尤推定量の周辺でどれだけ凸型になっているかを表している(詳しくはフィッシャー情報量 - Wikipedia)。

したがって、最尤推定量は、真の値を平均とし、フィッシャー情報行列の逆行列を分散共分散行列として持つ正規分布に従うということである。上の式を推定量と真の値の差で書き直せば、次のようになる。

 \boldsymbol{b} - \boldsymbol{\beta} \sim N(0, \mathcal{I}^{-1})

この式から、最尤推定量が真の値とどれだけ離れている可能性があるかを評価することができるので、信頼区間を求めることができる。

例:パラメータ2つのロジスティック回帰

例として次のデータのようなデータ(  N=40 )を考える。

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

f:id:cyanatlas:20201018014744j:plain

パラメータを2つ持つロジスティック回帰を考える。

 logit(q_i) = \beta_0 + \beta_1 x_i

 i 番目のサンプル  Y_i が1となる確率  q_i は、切片  \beta_0 とパラメータ  x_i の係数  \beta_1 によって決まるというモデルである。

ちなみにデータを生成するときに使ったパラメータ(真のパラメータ)は \beta_0 = -2.0, \beta_1 = 0.8 である。

では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メソッドを使うことで、フィッシャー情報行列を得ることができる(おそらく符号が逆になっていると思われたので、ソースコードの中では符号を反転させている)。

以下は  \beta_1区間推定を行う例。フィッシャー情報行列の逆行列を求め、その  (1, 1) 成分(  \beta_1 - b_1 の分散)から標準偏差を求める。±標準偏差×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%信頼区間になっているか、データの生成と信頼区間の推定を繰り返して検証してみる。

先程と同じく  N=40 のデータにロジスティック回帰を行って \beta_1 の信頼区間を求める。これを M=1000 回繰り返して、信頼区間の的中率を計算する。

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()

結果の一例を示す。

f:id:cyanatlas:20201018022108j:plain

この例では、的中率は96.5%であり、おおよそ95%になっている。

正規分布近似とWald統計量

以上の方法では、GLMのパラメータの最尤推定量は真の値を平均とする正規分布に従うという法則を用いていた。しかし、これは近似であるということに注意が必要である。データの数が少ないときには近似の精度が悪くなるので、正規分布に近似する代わりに、パラメトリックブートストラップを行う方が良いかもしれない。

最尤推定量と真の値の差が正規分布に従うのであれば、その二乗和はカイ二乗分布に従う(パラメータは独立でないので、分散共分散行列の逆行列を挟む必要がある。詳しくはDobson(2008, p.87)など)。この考え方で求まるのが、Wald統計量であり、このWald統計量を用いてWald検定が行われる。したがって、Wald検定は上で説明した信頼区間と類似した考え方だろうと思うのだが、詳細を知らないので確実なことは言えない。信頼区間の場合は分布の両側から2.5%を削って求めたが、Wald検定の場合はパラメータ  \beta が0から離れているかがポイントになるので、0に近い方から分布を5%削って有意性を調べるのだろうか。この辺りは久保(2012, p.51)に簡単な解説がある。とは言え、多くの場合はWald検定ではなく尤度比検定を使うだろうと思うので、この辺りの詳細はあまり気にする必要はないかもしれない。

参考文献

  • Annette J, Dobson (2008)『一般化線形モデル入門』(訳:田中豊・森川敏彦・山中竹春・冨田誠、原著第二版)共立出版
  • 久保拓弥 (2012)『データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC』(シリーズ 確率と情報の科学)岩波書店