7.5 適合度統計量

この章を理解する上で参考になるページ

  • パラメータを推定するために, 最尤推定を用いる代わりに,以下のような重み付き平方和の最小化を考える場合もある

\[S_w=\sum_{i=1}^{N}\frac{(y_i-n_i\pi_i)^2}{n_i\pi_i(1-\pi_i)}\]

  • ここに
    • \(E(Y_i)=n_i\pi_i\),
    • \(var(Y_i)=n_i\pi_i(1-\pi_i)\)

\(X^2\)(ピアソンカイ２乗統計量)

上記の方法は, ピアソンカイ２乗統計量（Pearson chi-squared statistic）を最小化することに同値

\[X^2=\sum_{i=1}^{n} \frac{(o-e)^2}{e}\]

証明

  • カイ二乗適合度統計量(\(X^2\))において,
    • \(o\): 観測度数（\(y_i\), \(n_i-y_i\)に該当）
    • \(e\): 期待度数（\(n_i\pi_i\), \(n_i(1-\pi_i)\)に該当)
  • \(X^2\)値を展開[右辺第１項は”成功”(\(pi_i\)), 第２項は”失敗”(\(1-\pi_i\))]

\[ \begin{align} X^2 &= \sum_{i=1}^{N}\frac{(y_i-n_i\pi_i)^2}{n_i\pi_i}+\sum_{i=1}^{N}\frac{[(n_i-y_i)-n_i(1-\pi_i)]^2}{n_i(1-\pi_i)} \\ &= \sum_{i=1}^{N}\begin{Bmatrix} \frac{(y_i-n_i\pi_i)^2(1-\pi_i)}{n_i\pi_i(1-\pi_i)}+\frac{[(n_i-y_i)-n_i(1-\pi_i)]^2\pi_i}{n_i\pi_i(1-\pi_i)} \end{Bmatrix} \\ &= \sum_{i=1}^{N}\begin{Bmatrix} \frac{(y_i-n_i\pi_i)^2(1-\pi_i)+(n_i\pi_i-y_i)^2\pi_i}{n_i\pi_i(1-\pi_i)} \end{Bmatrix} \\ &= \sum_{i=1}^{N}\begin{Bmatrix} \frac{(y_i-n_i\pi_i)^2(1-\pi_i)-(y_i-n_i\pi_i)^2\pi_i}{n_i\pi_i(1-\pi_i)} \end{Bmatrix} \\ &= \sum_{i=1}^{N}\frac{(y_i-n_i\pi_i)^2}{n_i\pi_i(1-\pi_i)} (1-\pi_i+\pi_i) \\&=S_w \end{align} \]

\(D\)(逸脱度適合度統計量)

推定期待値度数（ロジスティック回帰による推定値\(\hat{\pi_i}\)）を用いて\(X^2\)を評価したときの統計量 \[X^2=\sum_{i=1}^{N}\frac{(y_i-n_i\hat{\pi_i})^2}{n_i\hat{\pi_i}(1-\hat{\pi_i})}\tag{7.6}\] は, 式7.5の逸脱度 \[D=2\sum_{i=1}^{N}[y_ilog(\frac{y_i}{n_i\hat{\pi_i}})+(n_i-y_i)log(\frac{n_i-y_i}{n_i-n_i\hat{\pi_i}})]\] に漸近的に等しい (\(D=2\sum　olog\frac{o}{e}\))

証明

  • \(f(s)=slog(s/t)\)を, tを中心にテイラー展開した式
    \[slog\frac{s}{t}=(s-t)+\frac{1}{2}\frac{(s-t)^2}{t}+…\] を用いる

  • この式から \[\begin{align} D &= 2\sum_{i=1}^{N} \begin{Bmatrix} (y_i-n_i\hat{\pi_i})+ \frac{1}{2}\frac{(y_i-n_i\hat{\pi_i})^2}{n_i\hat{\pi_i}}+[(n_i-y_i)-(n_i-n_i\hat{\pi_i})]+ \frac{1}{2}\frac{[(n_i-y_i)-(n_i-n_i\hat{\pi_i})]^2}{n_i-n_i\hat{\pi_i}}+… \end{Bmatrix} \\ &\cong\sum_{i=1}^{N}\frac{(y_i-n_i\hat{\pi_i})^2}{n_i\hat{\pi_i}(1-\hat{\pi_i})} \\ &=X^2 \end{align}\] となり, 両者が漸近的に一致することがわかる

  • 仮定したモデルが正しいとき,Dは漸近的に\(D～\chi^2(N-p)\)となる
    上式から近似的に\(X^2～\chi^2(N-p)\)

  • \(D\)と\(X^2\)どちらを使うべき？
    　データのフォーマットに依る
    • 事象・試行データ：特定の共変量パターンにおける試行回数が多いため,応答変数は離散的な値をとる　　 　
    • ２値データ：試行回数が極めて少ない（1回）のため,応答変数は0か1の２値をとる
      　試行回数が多いとき（事象・試行データ）では \(X^2\)のほうが良い近似を与える(Cressie & Read 1989)
      試行回数が極めて少ないとき（2値データ）では\(D\)や\(X^2\)を用いることはできなくなる
      この場合,後述のホズマー・レメショウ統計量(かAIC)を用いるのが良い

ホズマー・レメショウ統計量(Hosmer & Lemeshow 1980)

2値データのロジスティック回帰において適用可能な適合度統計量

がん疫学の分野で使われるらしい(Wikiより) 詳細は7.8にて

手順
まず, 予測確率に従いデータを適当なグループに分割(g=10程度)
各グループにおける２値カテゴリ（成功,失敗など）の各度数を表にまとめる　この分割表に対してピアソンの\(X^2\)を計算し, 適合度の指標とする
→ホズマー・レメショウ統計量\(X^2_{HL}\)
\(X^2_{HL}～\chi^2(g-2)\)となることが推知実験により明らか
\(H^0\): 当てはめたモデルは正しい(p値が大きければ\(H^0\)の棄却を保留)

\[X^2_{HL}=\sum_{j=1}^{k}\frac{(O_j-E_j)^2}{N_j\hat{\pi_j}(1-\hat{\pi_j})}=\sum_{j=1}^{k}\frac{(O_j-N_j\hat{\pi_j})^2}{N_j\hat{\pi_j}(1-\hat{\pi_j})}～\chi^2(g-2)\]

library(ResourceSelection)

  set.seed(123)
n <- 500
x <- rnorm(n)
y <- rbinom(n, 1, plogis(0.1 + 0.5*x))
m <- glm(y ~ x, family=binomial)
hoslem.test(m$y, fitted(m))#p>0.05のとき適合しているとみなす
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  m$y, fitted(m)
## X-squared = 4.5227, df = 8, p-value = 0.8072

当てはめモデルの対数尤度関数と一定モデルの対数尤度関数の比較

  • いわゆる“カイ２乗近似の尤度比検定？”

  • 一定モデルとは, 全ての\(\pi_i\)が等しいようなモデル
    唯一のパラメータは\(\tilde{\pi}=(\sum y_i)/(\sum n_i)\)と推定される
    いま, p個のパラメータを含むモデルを当てはめたときに,\(Y_i\)の予測値が\(\hat{\pi_i}\)と推定された とする（当てはめ値は\(\hat{y_i}=\pi_i\hat{\pi_i}\)となる）
    このとき, 対数尤度関数lの差で定義された統計量 \[C=2[l(\hat{\pi; y})-l(\tilde{\pi};y)]\] は,
    \[c=2\sum[y_ilog(\frac{\hat{y_i}}{n_i\tilde{\pi}})+(n_i-y_i)log(\frac{n_i-\hat{y_i}}{n_i-n_i\tilde{\pi}})]\] となる

  • 5.5節の結果より,Cは帰無仮説\(\beta_2=…=\beta_p=0\)のもとで近似的に\(C～\chi^2(p-1)\)となり,対立仮説のもとでは非心カイ二乗分布に従う
    したがって, Cを「一定モデルが正しい」という帰無仮説の検定に利用できる

m0 <- glm(y~1,family=binomial)
anova(m0,m,test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: y ~ 1
## Model 2: y ~ x
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       499     692.95                          
## 2       498     651.54  1   41.406 1.236e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Cは尤度比カイ2乗統計量（likelihood ratio chi-squared statistic）とよばれる

擬似\(R^2\)

重回帰における決定係数（\(R^2\)）からの類推で, 擬似\(R^2\)が定義される \[pseudoR^2=\frac{l(\tilde{\pi};y)-l(\hat{\pi};y)}{l(\tilde{\pi};y)}=1-\frac{(\hat{\pi};y)}{(\tilde{\pi};y)}\] これは一定モデルに比べ当てはめモデルが対数尤度をどれだけ改善しているかをはかっている

library(BaylorEdPsych)

PseudoR2(m)
##         McFadden     Adj.McFadden        Cox.Snell       Nagelkerke 
##       0.05975417       0.05109550       0.07947670       0.10598307 
## McKelvey.Zavoina           Effron            Count        Adj.Count 
##       0.10232192       0.08094056       0.61800000       0.22040816 
##              AIC    Corrected.AIC 
##     655.54068258     655.56482744

この出力結果に表示されている 擬似^2について

7.6残差統計量

  • ロジスティック回帰における適合度指標Dや\(\chi^2\)に対応して用いられる２つの残差統計量について

    • \(n_k\)を第k共変量パターンにおける試行数,
    • \(Y_k\)を第k共変量パターンにおける成功数,
    • \(\hat{\pi_k}\)を第k共変量パターンにおける推定成功確率とする

  • ピアソン残差（Peason residuals）は, 以下のように定義される

\[X_k=\frac{(y_k-n_k\hat{\pi_k})}{\sqrt{n_k\hat{\pi_k}(1-\hat{\pi_k})}}, k=1,...,m\tag{7.7}\]

  • 式7.6より, \(X^2=\sum_{i=1}^{N}\frac{(y_i-n_i\hat{\pi_i})^2}{n_i\hat{\pi_i}(1-\hat{\pi_i})}\)
    であるから,
    \(\sum_{k=1}^{m}X^2_k=X^2\) である

  • さらに, 標準化ピアソン残差（standardized Pearson residuals）は,
    \[rP_k=\frac{X_k}{\sqrt{1-h_k}}\] で定義される. ただし\(h_k\)はハット行列から得られるてこ比（6.2.6参照）

  • 逸脱度残差（deviance residuals）も以前と同様に,

\[d_k=sign(y_k-n_k\hat{\pi_k}) \begin{Bmatrix} 2[ y_klog(\frac{y_k}{n_k\hat{\pi_k}})+(n_k-y_k)log(\frac{n_k-y_k}{n_k-n_k\hat{\pi_k}}) ] \end{Bmatrix}^{1/2}\tag{7.8} \]

式7.5より
\(D=2\sum_{i=1}^{N}[y_ilog(\frac{y_i}{\hat{y_i}})+(n_i-y_i)log(\frac{n_i-y_i}{n_i-\hat{y_i}})]\)であるから,

\(\sum_{k=1}^{m}d_k^2=D\)である

  • 標準化逸脱度残差（standardized deviance residuals）は

\[rD_k=\frac{d_k}{\sqrt{1-h_k}}\] で定義される

これらの残差統計量はモデルの妥当性をチェックするために利用される

  • 横軸に説明変数（連続）, 縦軸に残差をプロット：線形性の仮定が成立しているかを検討可能
  • 測定順に残差をプロットし, 系列相関の有無を確認

共変量パターンのほとんどでデータが2値の場合, 残差プロットはあまり意味がない
→互いに異なる値をもつ残差の数が少ないため

\(X^2\)や\(D\)などのモデル全体としての適合度統計量などを用いて判断すべし