ロジスティック回帰３

前回は、残差分析を行うといいつつ、予測値の信頼区間をもとめてしまいました。さて、今回は、前回の課題であった、残差分析を行います

まず、モデルの全体的な当てはまりは、次の式を用いて確認することが出来ます

(1) X^2 = Σ（観測度数 - 当てはめ値）^2 / 当てはめ値

(2) G^2 = 2Σ（観測度数)ln(観測度数 / 当てはめ値)

(1)は、普通のカイ二乗分析と同じです。式を見れば分かるとおり、(1)であらわされるカイ二乗統計量が小さければ小さい程、モデルがデータに良く当てはまってることが分かります。このときの自由度は、

カテゴリー数 - モデルのパラメータ数

で与えられます

(2)は、パラメータを計算する尤度比に基づく方法です。こちらも、カイ二乗分布に従います
ここで重要なことは、(1)も(2)も大標本下で、カイ二乗分布に従うということです。よって、前回計算したように、独立変数が、ほぼ連続変数の場合、この適合度検定は、不適切です。

と言うことで、今回は、独立変数を、カテゴリカルな変数に変換した上で、 モデルを再度当てはめます。その上で、個々のカテゴリーに対して、ピアソン残差を計算し、その結果が、(1)と一致することを確認します。

============================================================
まず最初に、忘却の彼方にあるデータ形式の確認

head(crab.data)
color spine width satell weight satell.flag width.cate
1 3 3 28.3 8 3050 1 28.25-29.25
2 4 3 22.5 0 1550 0 23.24
3 2 1 26.0 9 2300 1 25.25-26.25

そう、こんな感じでしたね。
さて、一番右端の変数ごとに、widthの平均値を取る

> aggregate(crab.data$width,list(y=crab.data$width.cate),mean)
y x
1 23.24 22.69286
2 23.25-24.25 23.84286
3 24.25-25.25 24.77500
4 25.25-26.25 25.83846
5 26.25-27.25 26.79091
6 27.25-28.25 27.73750
7 28.25-29.25 28.66667
8 29.25 30.40714

xが、カテゴリーの平均値です。これを、独立変数にします

NofRec NofSatell SampleRatio width.cate width.mean
23.24 14 5 0.3571429 23.24 22.69286
23.25-24.25 14 4 0.2857143 24.00 23.84286
24.25-25.25 28 17 0.6071429 25.00 24.77500
25.25-26.25 39 21 0.5384615 26.00 25.83846
26.25-27.25 22 15 0.6818182 27.00 26.79091
27.25-28.25 24 20 0.8333333 28.00 27.73750
28.25-29.25 18 15 0.8333333 29.00 28.66667
29.25 14 14 1.0000000 29.25 30.40714

最初のデータから、こんな感じのデータを再作成しておいて、
widht.meanを独立変数、SamplRatioを目的変数とするモデルを当てはめます。ただし、今回は、目的変数に確率を与えるので、ウ エイトにレコード数（＝サンプル数)を、当てておきます。でないと、計算してくれません。変数をアタッチしておけば、後が楽ですね

> attach(crab.category)
> glm(SampleRatio~width.mean,weights=NofRec,binomial)

結果は、こんな感じ

Call: glm(formula = SampleRatio ~ width.mean, family = binomial, weights = NofRec)

Coefficients:
(Intercept) width.mean
-11.5330 0.4654

Degrees of Freedom: 7 Total (i.e. Null); 6 Residual
Null Deviance: 34.03
Residual Deviance: 5.956 AIC: 33.14

パラメーターも、前回と殆ど変わりません
良かったですね

さて、ここで、尤度比検定を確認しましょう
ここでは、β=0としたとき尤度比と、そうでない時の尤度比との差を計算することで、当てはまりのよさ、つまり逸脱度(deviance)が検討できます
まず、β=0としたとき尤度比が

Null Deviance: 34.03

です。ちゃんとNullと書いてあります
またβ≠0の時の逸脱度が

Residual Deviance: 5.956

です
つまり2つの差

Null Deviance - Residual Devian

こそが、βの効果と言うわけです
で、この差は、

34.03 - 5.956 = 28.074

は、2つのモデルの自由度の差(7-6=1)のカイ二乗分布に従います
とうぜん、婚だけ大きな数であれば、1％未満で、βは有意です
glmの計算結果を変数に入れておいて、anovaを実行すれば、
引き算するまでもなく、同じ結果を返してくれます

> anova(crab.category.result, test="Chisq")
Analysis of Deviance Table

Model: binomial, link: logit

Response: SampleRatio

Terms added sequentially (first to last)


Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 7 34.034
width.mean 1 28.078 6 5.956 1.165e-07


devianceが28.078で、ちゃんとあってますね
p値もとっても小さいのでOKです
さて、尤度比による全体の当てはまりが確認できました。
当てはまりの確認は、ピアソン残差でも計算できます
ピアソン残差とは、個別のデータに対して計算します

ei = (yi - niπi) / sqrt(ni*πi(1-πi))

そして、

x^2 = Σei^2

となります。
なんだか、わけの分からない式ですね
でも、しょうがないです。あきらめましょう
glmを当てはめたデータの右端に、予測確率のfittedを付けておきました

NofRec NofSatell fitted
23.24 14 5 0.2745370
23.25-24.25 14 4 0.3925719
24.25-25.25 28 17 0.4993264
25.25-26.25 39 21 0.6206342
26.25-27.25 22 15 0.7181918
27.25-28.25 24 20 0.7983565
28.25-29.25 18 15 0.8591791
29.25 14 14 0.9320432


このデータから、ピアソン残差を計算してみます
それが、こんな感じです

peason.resid <- (NofSatell-NofRec*fitted)/sqrt(NofRec*fitted*(1-fitted))

[1] 0.6925753 -0.8187713 1.1410233 -1.0575784 -0.3792288 0.4270668 -0.3152460 1.0103285

つづいて、二乗和を計算します

sum(peason.resid^2)
[1] 5.016797

この値が、式(1)と同じになるはずです
ということで、この値は、カイ二乗分布に従います
今回は、P=0.4位になって、当てはまりはオッケーです

長くなりましたが、こんな感じです。
なお、ピアソン残差に、てこ比を当てた値を標準化残差といいますが、テキストに式が載ってませんでした。なお、標準化残差は、正規分布に従うそうです。でも、今回は、無視ですわ。
AICは、そのうちって感じですかね。

次回は、多重ロジスティック回帰モデルにチャレンジします

=========================================================
カテゴリカルデータ解析入門
Alan Agresti
サイエンティスト社

のp.153～p.163あたりを挑戦

（感想)
(1) X^2 = Σ（観測度数 - 当てはめ値）^2 / 当てはめ値

式の値が合わなくて苦労した
結局、あわずじまい
でも、テキストが間違ってる気がする

というのも、(1)のX^2は、ピアソン残差の2乗和と一致する
実際、ピアソン残差は、x^2と一致した。
で、ピアソン残差の分母は、標本数に、２項分布の標準偏差をかけたものなのに、(1)式の分母は、単なる当てはめ値としか書いてない
だから、計算が合わなかった。

たぶん、(1)の分母は、ピアソン残差の平方根のことだと思う