Devienceの和訳には「尤離度」があり、以前金先生も使用しておられたが、あまり一般的ではなく、「逸脱度」が最近は主流である。以下では原語のまま使用する。
尤度を用いる最適化では、何らかの条件で対数尤度\( \log L \)を最大化して、
最尤推定値とする。パラメータ\( \beta_0,\beta_1,...,\beta_p \)を、
観測された data の元で尤度を最大にする値が
\( \hat{\beta}_0,\hat{\beta}_1,...,\hat{\beta}_p \)であるとき、
それを代入した最大値である\( \log L ((\hat{\beta})|data) \)を
\( \log L^* \)と書くことにすれば、そのモデルのdevience D は単純であり
\[ D = -2 \log L^* \]
である。そして、飽和モデル(saturated model)、即ちその設定の元でパラメータを可能なかぎり
（データ数まで）増やしたモデルのdevience を \( D_{sat} \)とすれば、それは
その設定の元での最小のdevienceである（最大対数尤度が最大）。このとき
問題にしているモデルのresidual devience RD を
\[ RD = D - D_{sat} \]
と定義する。そして、residual devienceが最大になる、つまり最も
あてはまりの悪いモデル（回帰ならば定数項のみ）を NULL modelと
よび、その時の devience を \( D_{null} \) とする。Null devience ND はその
redidual devience であり、つまり
\[ ND = D_{null} - D_{sat} \]
によって定義される。R における glm()の出力での
Null Devience: ***
Residual Devience: *** AIC: ***
となっている2つの Devience はその意味である。
また、AIC は、\( D_{sat} \)による調整は行わず
\[ AIC = D + 2 \#(parameters) = -2(\log L^* - \#(parameters) ) \]
と定められている。\( \log L^* \)が大きければ予測は良いのだが、
パラメータの数というバイアスが有り、それを調整してAICとなる。
AICは大小関係をみるのだから、D でなくRDを
もちいても、比較の用途には同じことになる。つまり、次のような
関係となっている。
\[ AIC = RD + D_{sat} + 2 \#(parameters) \]
なお、クラスター分析の距離の補足で、Kalback-Leibler情報量に ふれる予定があるが、そのKLは対数尤度の期待値のようなものである。 （距離の場合は対称化して、Information Distance とする）
GLMにおいて、Residual devience が何か解りにくいと 思う場合は、定義を近似式に表すとよい。 また、\( \chi^2 \)検定と尤度によるG検定の関係をわかっていればなおよい。 式変形によれば確かに２次項の和が主要部になり、通常の ２乗和のResidual と、残差の高次項を除いて同様になる。
まず基本はTaylor展開にある。即ち -1<\( t \)<1として、
\[ \log(1+t) = t - \frac{1}{2}t^2 + \mathcal{O}(t^3), \] \[ \log(1-t) = -t - \frac{1}{2}t^2 + \mathcal{O}(t^3). \]
となる。ここで \( \mathcal{O} \) は Landau symbol.
Pearson の\( \chi^2 \)検定量は、分割表の各マス目の出現頻度を\( O_i \),
独立性の帰無仮説の元での期待値を\( E_i \)とすれば次式であった。ここで
分母の\( E_i \)はPoisson分布での分散（平均に等しい）に該当する。要するに、
各項は標準化残差（調整済み残差ではない）の２乗和である。
\[ \chi^2 = \sum_i \frac{(O_i - E_i)^2}{E_i} \]
そして、G統計量（尤度）は
\[ G = -2 \sum_i O_i \log \frac{E_i}{O_i} \]
である。ここで
\[ \log \frac{E_i}{O_i}=\log\left( 1 - \frac{O_i-E_i}{O_i} \right) \]
となることから、\( | O_i-E_i |/O_i \)が十分小さいとして
\[ G = \sum \left( 2 (O_i-E_i) + \frac{(O_i-E_i)^2}{O_i}
+ \mathcal{O}((O_i-E_i)^3) \right) \]
となる。さらに、\( \sum (O_i-E_i)=0 \) と
\( 1/O_i = (1/E_i)(1 - (O_i-E_i)/O_i) \)
を用いれば、最終的に
\[ G = \sum \left( \frac{(O_i-E_i)^2}{E_i} +
\mathcal{O}((O_i-E_i)^3) \right) \]
となる。この初項はまさに\( \chi^2 \)値である。
つまり、\( \chi^2 \)値は、尤度の近似主要部なのであった。
そして、独立性の帰無仮説（marginal分布からのランダム標本）の
もとで、Gはしかるべき自由度（\( \chi^2 \)と自由度と同じ）の\( \chi^2 \)
分布で近似される。
なお、上記の近似は、\( |O_i-E_i|>E_i \)となる i があると成り立たない。
この場合、\( \chi^2 \)検定よりもG検定を用いるのが望ましいとされる。
ここでは、誤差に二項分布を設定する場合の説明を行う。
B(p,n)の生起確率は下記の通り。
\[ P(y) = {}_n C_y \,p^y(1-p)^{n-y}, \ y=0,1,...,n . \]
そして GLM(family=binomial(logit))は次のモデルである。
\[ logit(p) = \beta_0 + \beta_1 x_1 + ... + \beta_p x_p. \]
\[ logit(p) = \log\frac{p}{1-p}, \ logistic(t)=\frac{1}{1+e^{-t}}.
logit(logistic(t))=t, \ logistic(logit(p))=p. \]
まず個別のDevience Residuals \( dr_i \) を考察する。 \( y_i \) を二項分布\( B(p_i,n_i) \)に従う確率変数の観測度数 、 \( \hat{y}_i \)をGLMモデル(family=binomial)による予測度数とするとき、 \( dr_i \)は符号を残差推定値 \( \hat{e}_i = y_i - \hat{y}_i \) と同じとして
\[ dr_i^2 = -2 \left(y_i \log\frac{\hat{y}_i}{y_i} + (n_i-y_i) \log\frac{n_i-\hat{y}_i}{n_i-y_i} \right) \]
にて定義される。２値の応答変数、即ち\( n_i=1 \)の場合は、\( y_i \)は0,1 のみをとる変数であり、\( \hat{p}_i \)をモデルに依る予測確率とすれば 次式のように単純化される。
\[ dr_i^2 = -2 \left(y_i \log \hat{p}_i + (1-y_i) \log(1-\hat{p}_i) \right) \]
ここで一般の場合に残差\( \hat{e}_i \)が十分小であるとして、 \( \hat{y}_i/y_i=1-\hat{e}_i/y_i \)などの変形を、 logの展開式に入れて整理すれば、\( \hat{e}_i \)の一次項は互いに 相殺し次式を得る。 \[ dr_i^2 = \frac{n_i}{y_i(n_i-y_i)}\hat{e}_i^2 + \mathcal{O}(\hat{e}_i^3). \] 観測値からの生起確率 \( y_i/n_i \) を予測値\( \hat{p}_i \)で置き換えれば、符号を込めて次の通り。 \[ dr_i = \frac{\hat{e}_i}{\sqrt{n_i \hat{p}_i (1-\hat{p}_i)}}+ \mathcal{O}(\hat{e}_i^2). \] ここで、\( n_i \hat{p}_i (1-\hat{p}_i) \)は分散の予測値による 推定値\( \hat{\sigma}^2 \)である。従って、この主要項は通常で言うならば”標準化残差", あるいは分散を標本からの推定値を用いたという意味では “スチューデント化残差”なのである。 （但し、\( \hat{\sigma}^2 \)は母分散推定値であって、\( \hat{e}_i \)そのものの 分散は、線形回帰のときのように簡単には求まらない。） そして、これらの総和が Residual Devience であり、 \[ RD = \sum_{i=1}^n dr_i^2 = \sum \left( \frac{\hat{e}_i^2}{\hat{\sigma}^2} + \mathcal{O}(\hat{e}_i^3) \right). \] この主要項は、”標準化残差（スチューデント化残差）”の２乗和である。
Poisson分布(link=log)を仮定するGLMでのdevience residuals は次式で与えられる。ここで、\( \hat{y}_i \)は、誤差分布Pois(\( \lambda_i \))の母数\( \lambda_i \)の元での GLM回帰モデルによる予測値であり、\( \hat{e}_i \)は残差推定値である。
\[ dr_i^2 = -2(y_i - \hat{y}_i + y_i \log \frac{\hat{y}_i}{y_i}) = \frac{\hat{e}_i^2}{y_i}+\mathcal{O}(\hat{e}_i^3) . \] 第一項の符号付き平方根\( dr_i \)が、標準化残差に該当する。 \[ dr_i = \frac{\hat{e}_i}{\sqrt{y_i}} +\mathcal{O}(\hat{e}_i^2). \] なお、\( y_i=0 \) の対策としては、分母をその予測値\( \hat{y}_i \)で 取り替える方がよいであろう。\( 1/y_i = (1-\hat{e}_i/y_i)/\hat{y}_i \) であるから、近似式には影響しない。 その表記でのResidual devience は次の通り。 \[ RD = \sum_{i=1}^n \left( \frac{\hat{e}_i^2}{\hat{y}_i}+\mathcal{O}(\hat{e}_i^3) \right). \] 念のため。生起確率と対数線形モデルは次の通り。 \[ P(y_i) = \frac{\lambda_i^{y_i} e^{-\lambda_i}}{y_i !} , \ y_i =0,\,1,\,2... \] \[ \log(\lambda) = \beta_0 + \beta_1 x_1 + ... + \beta_p x_p . \]
参考のため、誤差分散が通常の仮定のもとで、通常の線形回帰を行う時に、対数尤度がどうなるかを考える。
定数項を含めた（元の）データベクトルを
\( x_i=(1,x_{i1}, x_{i2}, ..., x_{ip})^T, i=1,...,n \)と置く。係数は\( b=(b_0, b_1,...,b_p) \)である。このすべてをparameterとして採用する。
モデルは下記の通り。
\[ y_i = x_i^T b + \varepsilon_i, \ (i=1,2,...,n). \]
ここで、\( \varepsilon_i \)が等分散（\( \sigma^2 \)）の正規分布に従うとすれば、 \[ f(y_i) = \frac{1}{\sqrt{2 \pi \sigma^2}}\exp \left( - \frac{\varepsilon_i^2}{2 \sigma^2} \right) = \frac{1}{\sqrt{2 \pi \sigma^2}}\exp \left( - \frac{(y_i - x_i^T b)^2}{2 \sigma^2} \right). \]
尤度は\( f(y_i) \)の積となる。よって、対数尤度は
\[ \log L(b, \sigma^2 | x,y) = \sum_{i=1}^n \log f(y_i) = - \frac{n}{2} \log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_{i=1}^n (y_i - x_i^T b)^2. \]
ここで変数は\( b \)と\( \sigma^2 \)の p+2 個であり、 \( x_i, y_i \)は定数である。偏微分して、
\[ \frac{\partial \log L}{\partial (\sigma^2)}= - \frac{n}{2 \sigma^2} + \frac{1}{2(\sigma^2)^2}\sum \varepsilon_i^2, \]
\[ \frac{\partial \log L}{\partial b_j} = \frac{-1}{\sigma^2} \sum_{i=1}^n x_{ij} (y_i-x_i^T b). \] 第一式より、\( \sigma^2 = \frac{1}{n}\sum_{i} \varepsilon_i^2. \) 第二式は(p+1)変数 \( b \)に関する (p+1)個の１次方程式となり、解ける。 しかも、その連立方程式は 最小二乗法と同じ である。 従って、\( X^T=(x_1,...,x_n) \) とおけば \( b \)の最尤推定値 \( \hat{b}=(X^T X)^{-1} X^T y \)となる。さらにハット行列を \( H=X (X^T X)^{-1} X^T \)と定義すれば、\( \varepsilon \)の最尤推定値 \[ \hat{\varepsilon} = y-\hat{y}=y-X\hat{b}=(I - H)y = (I - H)(y- \bar{y}). \] それを第一式に代入して、誤差等分散の前提と前式より \[ \hat{\sigma^2}=\frac{1}{n}\sum_{i} \hat{\varepsilon_i}^2= (1-\frac{tr(H)}{n})\sigma^2=(1-\frac{p+1}{n})\sigma^2 \] と求まる。最終的に、対数尤度は次の通り。 \[ \log L(\hat b, \hat{\sigma}^2|x,y) = - \frac{n}{2} \log(2 \pi \hat{\sigma}^2) - \frac{1}{2 \hat{\sigma}^2} \sum_{i=1}^n \hat{\varepsilon_i^2} = - \frac{n}{2} \log(2 \pi \hat{\sigma}^2) -\frac{n}{2}. \]
これより devience と AIC は次の通りである。 \[ D = n \log(2 \pi \hat{\sigma}^2) + n , \] \[ AIC = D + 2 (p+2) = n \log(2 \pi \hat{\sigma}^2) + n + 2(p+2). \]
Residual Devienceを比較する検定には、anova ではなく car::Anovaを 用い、type=2 によって Type II により検定すること。また、unbalanced design の場合は Type III も実行して比較すること。
glm object については、method: coefficients, deviance, df.residual が用意されているので、objectの構造を知る必要はない。
method: predict をそのまま用いると、予測logit 関数値が返る。
予測確率値を求めるには、option: type=“response” を付加するか、
method fitted を使用する。
例えば、glmの結果obj から、dispersion の程度を見るには次のようにする。
(dp <- sum(residuals(obj, type = "peason")^2)/df.residual(obj))
狭義の Logistic Regression は、応答変数が01の２値のものをよぶ。
広義の場合は、応答変数が比率であり、それにlogistic curve をfit
させる回帰を呼ぶ。GLMであつかう限り、どちらも同じである。
２値とはいえ、0,1は数字としての意味がある。
さらに、GLMにおいては、説明変数に数値変数のみならず、名義尺度
（01などの数字表現を含む）が
まざっていても同じである（解釈に気をつけるべきところがある）。
plot.lm は、class glmにも適用され、診断プロットを行う。診断は ６通りあり、which= に指定する。Rcmdrで古典的４通りを出すこともできるが、なるべく６通りすべて出すのが好ましい。このためには
opar <- par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
for (cnt in 1:6) plot(obj, which = cnt, id.n = 5, cex.id = 1)
par(opar)
とでもすることに成る。ここで、RStudioの場合にも合うように、2行3列の横長の形式とした。２段組の論文ならmfrow=c(3,2)がよいかもしれない。 またid.n=5 (default は3）とは、ラベルを 出す外れ値の件数である。そのサイズを cex.idで指定する（defaultは 0.85）。単に plot(obj)とすると、入力待ちになる。
dataframe dat に、２値あるいはm値のデータが数字で表現されている 列名 cola, colb があるとする。属性は通常 INT であろう。この２列の クロス集計を取るときに、個々を因子にした列を追加する方式で 行うと、Rcmdr では再度 dat を読みこまねばならず、一つしか データセットがないときには、次の２つの方法のどちらかとなる。
いずれも至って面倒である。直接つぎのようにするのがよいであろう。
(ab.t <- with(dat, table(namea = as.factor(cola), nameb = as.factor(colb))))
すると、dat$cola と dat$colb のクロス集計が取られ、注釈として namea, nameb が表示される。namea は cola とは別の文字列でも良い（引用符はつけない）。この程度は Rcmdr に依存せずにコマンドでしたほうが良い。