対数尤度比統計量(逸脱度)を使ってモデルの評価をしよう
モデルの評価をするために何をすればよい？
→フルモデルと興味のあるモデルを比較
ただし, 興味のあるモデルとフルモデルは同じ確率分布およびリンク関数をもつ一般化線形モデルとする
N個の線形線分\(x_i^T\beta\)に対してN個の観測地\(Y_i,i=1,...,N\)が得られているとする.このときフルモデルはN個のパラメータで示される
\(L(b_{max}; y)\) はフルモデルの尤度関数,
\(L(b; y)\)は関心のあるモデルの尤度関数の最大値とするとき,
尤度比
\[\lambda=\frac{L(b_{max}; y)}{L(b; y)}\] がモデルの適合度を測る方法を与える
実際は扱いやすい尤度比の対数,すなわち対数尤度関数の差を用いる
\[log\lambda=l(b_{max}; y)-l(b; y)\]
\(log\lambda\)が大きい,すなわちフルモデルの対数尤度が関心のあるモデルの対数尤度より大きいとき, 関心のあるモデルはフルモデルに比べてあまりよくデータを記述していないことを意味する
\(log\lambda\)の棄却域を決めるには標本分布が必要であるが, 一般に\(log\lambda\)に2をかけた\(2log\lambda\)カイ二乗分布に従うことを利用する
\(2log\lambda\)は逸脱度(Deviance)とよばれる
逸脱度がカイ二乗分布に従うことを確認しよう
逸脱度（deviance）または対数尤度比統計量（rog likelihood ratio statistics）は以下のように定義される：
\[D=2[l(b_{max}; y)-l(b; y)]\]
\(b\)がパラメータ\(\beta\)の最尤推定量ならば式（5.4）より(U(b)=0より)近似的に以下が成り立つ
\[l(\beta)-l(b)=-\frac{1}{2}(\beta-b)\Im(b)(\beta-b)\] したがって
\[2[l(b; y)-l(\beta; y)]=(\beta-b)^{T}\Im(b)(\beta-b) \]
となり, 式（5.7）より自由度pのカイ二乗分布\(\chi^2(p)\)に従う
p:関心のあるモデルのパラメータ数
以上から,逸脱度 \[
\begin{align}
D &= 2[l(b_{max}; y)-l(b; y)] \\
&= 2[l(b_{max}; y)-l(\beta_{max}; y)]-2[l(b; y)-l(\beta; y)]+2[l(\beta_{max}; y)-l(\beta; y)]
\end{align}
\] の標本分布が導かれる
うち
第1項は\(\chi^2(m)\)に従う(m: フルモデルのパラメータ数) 第2項は\(\chi^2(p)\)に従う(pは関心のあるモデルのパラメータ数)
第3項はフルモデルと関心のあるモデルのあてはまりの程度が近いとゼロに近づく
ゆえに
\[D～\chi^2(m-p, v)\] ただしvは非心パラメータ
応答変数が正規分布に従うとき, Dは正確にカイ二乗分布に従うが, \(\sigma^2\)が未知であるため, 直接利用することはできない
応答変数が他の分布に従う場合はDの標本分布は近似的にカイ二乗分布となるが, 2項分布やポアソン分布なのではDは実際に計算可能である
応答変数\(Y_1 ,..., Y_N\)が互いに独立に2項分布に従う\(（Y_i ～ binomial (n_i ,\pi_i )）\)とき，対数尤度関数は
\[l(\beta; y)=\sum_{i=1}^{N}\left[y_ilog\pi_i-y_ilog(1-\pi_i)+n_ilog(1-\pi_i)+log\begin{pmatrix}n_i\\p_i\end{pmatrix} \right]\]
と表される
フルモデルでは\(\pi_i\)がすべて異なるので, \(\beta=\left[\pi_1 ,..., \pi_N\right]\)
最尤推定値は\(\hat{\pi_i}=y_i/n_i\)となり, 対数尤度関数の最大値は以下の通り
\[l(b_{max}; y)=\sum\left[y_ilog\frac{y_i}{n_i}-y_ilog\frac{(n_i-y_i)}{n_i}+n_ilog\frac{(N_i-y_i)}{n_i}+log\begin{pmatrix}n_i\\p_i\end{pmatrix} \right]\] ここで,\(p<N\)となるような他のモデル(=関心のあるモデルのこと？)をあてはめ, \(\hat{\pi_i}\)を確率の最尤推定値, \(\hat{y_i}=n_i\hat{pi_i}\)を当てはめ値とする.
これらの値から求められた対数尤度関数は, \[l(b; y)=\sum\left[y_ilog\frac{\hat{y_i}}{\hat{n_i}}-y_ilog\frac{(n_i-\hat{y_i})}{n_i}+n_ilog\frac{(N_i-\hat{y_i})}{n_i}+log\begin{pmatrix}n_i\\p_i\end{pmatrix} \right]\] したがって逸脱度は,
\[ \begin{align} D &= 2[l(b_{max}; y)-l(b; y)] \\ &=2\sum_{i=1}^{N}\left[ y_ilog(\frac{y_i}{\hat{y_i}})+(n_i-y_i)log(\frac{n_i-y_i}{n_i-\hat{y_i}})\right] \end{align} \]
以下のモデルについて考える： \[E(Y_i)=\mu_i=x_i^T\beta ; Y_i～N(\mu_i, \sigma^2), i=1,...,N\] 対数尤度関数は
\[l(\beta; y)=-\frac{1}{2\sigma^2}\sum_{i=1}^{N}(y_i-\mu_i)^2-\frac{1}{2}Nlog(2\pi\sigma^2)\] となる
フルモデルではすべての\(\mu_i\)が異なりうるので, \(\beta\)はN個の要素\(\mu_1,...,\mu_N\)からなる.対数尤度比関数をおのおのの\(\mu_i\)について微分し推定方程式を解くと\(\hat{\mu_i}=y_i\)が得られる
ゆえにフルモデルに対する対数尤度関数の最大値は:
\[l(b_{max}; y)=-\frac{1}{2}Nlog(2\pi\sigma^2)\] パラメータ数が\(p<N\)である様な他のモデル(＝関心のあるモデル？)に対して最尤推定量は式（5.11）より \[b=(X^TX)^{-1}X^Ty\] となり,
対応する対数尤度関数の最大値は
\[l(b; y)=-\frac{1}{2\sigma^2}\sum(y_i-x_i^Tb)^2-\frac{1}{2}Nlog(2\pi\sigma^2)\] したがって逸脱度は
\[ \begin{align} D &= 2[l(b_{max}; y)-l(b; y)] \\ &=\frac{1}{\sigma^2}\sum_{i=1}^{N}(y_i-x_i^Tb)^2\\&=\frac{1}{\sigma^2}\sum_{i=1}^{N}(y_i-\hat{\mu_i})^2 \end{align} \] ここで\(\hat{\mu_i}\)は当てはめ値\(x_i^Tb\)である
パラメータが１つしかないという特別な場合 例えばすべてのiについて\(E(Y_i)=\mu\)が成り立つとき（＝一定モデル？） において, $XはN個の１からなるベクトルとなり, \(b=\hat{\mu}=\sum_{i=1}^{N}y_i/N=\bar{y_i}\)かつ, すべてのiについて\(\hat{\mu_i}=\bar{y}\)が成り立つ
よって
\[D=\frac{1}{\sigma^2}\sum_{i=1}^{N}(y_i-\bar{y})^2\] この統計量は標本分散\(S^2\)との間に \[S^2=\frac{1}{N-1}\sum_{i=1}^{N}(y_i-\bar{y})^2=\frac{\sigma^2D}{N-1}\]
のような関係があり, \((N-1)S^2/\sigma^2～\chi^2(N-1)\),
したがって,
\(D～\chi^2(N-1)\)は正確である
一般の場合について考えると,
\[\begin{align}
D &=\frac{1}{\sigma^2}\sum_{i=1}^{N}(y_i-x_i^Tb)^2 \\&=\frac{1}{\sigma^2}(y-Xb)^T(y-Xb)
\end{align}\] ここで, デザイン行列Xの行は\(x_i\)である
\((y-Xb)\)は, \[\begin{align}
y-Xb &=y-X(X^TX)~{-1}X^T\\&=[I-X(X^TX)^{-1}]y=[I-H]y\end{align}\] ここで,
\(H=X(X^TX)^{-1}X^T\)であり, ハット行列(‘hat’ matrix)とよばれる
うんぬんかんぬんで, (Dの２次形式の部分はHがべき等(すなわち\(H=H^T\)かつ\(HH=H\))だから)
\((y-Xb)^T(y-Xb)={[I-H]y}^T[I-H]y=y^T[I-H]y\)
うんぬんかんぬんで, （詳細はGraybill 1976）
Dは中心カイ二乗分布\(\chi^2(N-p)\)に従うといえる
未知の\(\sigma^2\)を含むDに対して
両辺に\(\sigma^2\)をかけて
\[\sigma^2D=\sum(y_i-\hat{\mu_i})^2\] を尺度つき逸脱度(scaled deviance)と呼ぶことがある
モデルがよく当てはまっているとき,Dは漸近的に\(\chi^2(N-p)\)にしたがい,その期待値はN-pである
したがって尺度パラメータ\(\sigma^2\)の１つの不偏推定量として以下が得られる
\[\tilde{\sigma^2}=\frac{\sum(y_i-\hat{y_i})^2}{N-p}\]
Dが直接計算できる場合
応答変数\(Y_1,...,Y_N\)が互いに独立で\(Y_i～Poisson(\lambda_i)\)のとき, 対数尤度関数は
\[l(\beta; y)=\sum y_ilog\lambda_i-\sum \lambda_i -\sum log y_i!\] フルモデルでは\(\lambda_i\)は全て異なるので\(\beta=[\lambda_1,...,\lambda_N]^T\)である
最尤推定値は\(\hat{\lambda_i}=y_i\)である, 対数尤度関数の最大値は
\[l(b_max; y)=\sum y_ilog y_i-\sum y_i -\sum log y_i!\] となる
関心のあるモデルのパラメータ数を\(p<N\)とすると 対数尤度関数の最大値は
\[l(b; y)=\sum y_ilog \hat{y_i}-\sum \hat{y_i} -\sum log y_i!\] となる
したがって逸脱度は
\[
\begin{align}
D &= 2[l(b_{max}; y)-l(b; y)] \\
&=2[\sum y_i log(y_i/\hat{y_i})-\sum (y_i-\hat{y_i})])\end{align}
\] と表される \(\sum y_i=\sum \hat{y_i}\)(演習問題9.1参照)から,
\[D=2\sum o_i log(o_i/e_i)\]
ただし観測値\(y_i\)を\(o_i\), 推定された期待値\(\hat{y_i}\)を\(e_i\)とする
このとき,Dはデータから計算可能である
yi <- c(2,3,6,7,8,9,10,12,15)
xi <- c(-1,-1,0,0,0,0,1,1,1)
pois <- data.frame(x=xi,y=yi)
fit <- glm(y~x, data=pois, family="poisson")
summary(fit)##
## Call:
## glm(formula = y ~ x, family = "poisson", data = pois)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8472 -0.2601 -0.2137 0.5214 0.8788
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.8893 0.1421 13.294 < 2e-16 ***
## x 0.6698 0.1787 3.748 0.000178 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 18.4206 on 8 degrees of freedom
## Residual deviance: 2.9387 on 7 degrees of freedom
## AIC: 41.052
##
## Number of Fisher Scoring iterations: 4長さpのパラメータベクトル\(\beta\)に関する仮説はワルド統計量\((\hat{\beta}-\beta)^T\Im(\hat{\beta})～\chi^2(p)\)を用いて検定できるワルド検定
まれに, スコア統計量\(U^T\Im^{-1}U～\chi^2(p)\)が使われることもあるスコア検定
別のアプローチ：２つのモデルのあてはまりの良さの比較
条件：２つのモデルは入れ子状（nested）あるいは階層的（hierarchical）になっている必要あり
モデル\(M_0\)に対応する帰無仮説
\[H_o: \beta=\beta_0=\begin{bmatrix} \beta_1 \\ \vdots\\ \beta_{q} \end{bmatrix} \] と, モデル\(M_1\)に対するより一般的な仮説 \[H_1: \beta=\beta_1=\begin{bmatrix} \beta_1 \\ \vdots\\ \beta_{p} \end{bmatrix} \] を考える(\(q<p<N\))
対立仮説\(H_1\)に対する帰無仮説\(H_0\)の検定は以下に示す逸脱度統計量の差を用いて行うことができる
\[\begin{align} \Delta D &= D_0-D_1 \\&=2[l(b_{max}; y)-l(b_0; y)]-2[l(b_{max}; y)-l(b_1 ; y)] \\ &=2[l(b_1; y)-l(b_0; y)] \end{align}\]
どちらのモデルもデータをよく記述しているとき, \(D_o～\chi^2(N-q)\)および
\(D_1～\chi^2(N-p)\)が成り立ち, したがって, \(\Delta D～\chi^2(p-q)\)となる
観測値から計算された \(\Delta D\) の値が\(\chi^2(p-q)\)分布に矛盾しないなら, より単純であるという理由から\(H_0\)に対応するモデル\(M_0\)を選ぶ
\(\Delta D\)の値が棄却域に入っていれば,モデル\(M_1\)のほうが\(M_0\)よりもデータを有意によく記述しているとして, \(H_0\)を棄却して\(H_1\)を採択する
データから逸脱度を計算できるとき：カイ二乗近似を利用
データから逸脱度を直接計算できないとき：どうする？
独立な確率変数\(Y_1,....,Y_N\)が正規線形モデル
\[E(Y_i)=\mu_i=x_i^T\beta ; Y_i ～ N(\mu_i, \sigma^2)\] に従うとき, 逸脱度は
\[D=\frac{1}{\sigma^2}\sum_{i=1}^{N}(y_i-\hat{\mu_i})^2\] である
\(\hat{\mu_i}(0)\) : \(M_0\)のあてはめ値
\(\hat{\mu_i}(1)\) : \(M_1\)のあてはめ値 とする
このとき
\[D_0=\frac{1}{\sigma^2}\sum_{i=1}^{N}[y_i-\hat{\mu_i(0)}]^2\]
\[D_1=\frac{1}{\sigma^2}\sum_{i=1}^{N}[y_i-\hat{\mu_i(1)}]^2\]
である
未知の\(\sigma^2\)を取り除くために,比を用いると：
\[\begin{align} F &=\frac{D_0-D_1}{p^q}/\frac{D_1}{N-p}\\ &=\frac {\frac{{\sum[y_i-\hat{\mu_i(0)}]^2-\sum[y_i-\hat{\mu_i(1)}]^2}} {(p-q)}} {\sum[y_i-\hat{\mu_i(1)}]^2/(N-p)} \end{align}\] \(H_0\)が正しくないとき, Fは分布\(F(p-q,N-p)\)から期待される値よりも大きくなる
d<-read.csv("2.2_weight_duration.csv")
summary(d)## duration weight sex
## Min. :35.00 Min. :2412 female:12
## 1st Qu.:37.00 1st Qu.:2785 male :12
## Median :38.50 Median :2952
## Mean :38.54 Mean :2968
## 3rd Qu.:40.00 3rd Qu.:3184
## Max. :42.00 Max. :3473head(d)## duration weight sex
## 1 40 2968 male
## 2 38 2795 male
## 3 40 3163 male
## 4 35 2925 male
## 5 36 2625 male
## 6 37 2847 malem0 <- lm(weight~duration+sex,data=d)
m1 <-lm(weight~duration*sex,data=d)
anova(m0,m1,test="F")## Analysis of Variance Table
##
## Model 1: weight ~ duration + sex
## Model 2: weight ~ duration * sex
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 658771
## 2 20 652425 1 6346.2 0.1945 0.6639