ロジスティック回帰分析

作者: 太郎丸博
出版社/メーカー: ナカニシヤ出版
発売日: 2005/07/01
メディア: 単行本
購入: 1人クリック: 26回
この商品を含むブログ (12件) を見る

　統計学の教科書に載っている女性の年齢と未婚率の関係の分析（太郎丸 2005: 173）を再現してみる。まずは個票データに対してロジスティック回帰分析を適用するやり方から。

Y	<- c(rep(0,  1), rep(1, 20),		# 0＝既婚，1＝未婚
             rep(0, 29), rep(1, 25),
             rep(0, 56), rep(1, 10),
             rep(0, 49), rep(1,  3),
             rep(0, 59), rep(1,  2),
             rep(0, 70), rep(1,  3),
             rep(0, 76), rep(1,  3))
X	<- c(rep(22, 21), 			# 年齢カテゴリの階級値
	     rep(27, 54), 
	     rep(32, 66), 
	     rep(37, 52),
             rep(42, 61), 
	     rep(47, 73), 
	     rep(52, 79))
d1	<- data.frame(Y, X)

　未婚か既婚かをあらわす変数Yを年齢Xに回帰させるロジスティック回帰式を立てて、最尤推定法をつかって推定する。

l1	<- glm(Y~ X, data = d1, family = binomial)
summary(l1)

Call:
glm(formula = Y ~ X, family = binomial, data = d1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5388  -0.4826  -0.1876  -0.1158   3.1648  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.08655    0.78789   6.456 1.08e-10 ***
X           -0.19400    0.02472  -7.848 4.23e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 360.44  on 405  degrees of freedom
Residual deviance: 251.52  on 404  degrees of freedom
AIC: 255.52

Number of Fisher Scoring iterations: 6

　推定結果を見ると $\mathrm{logit}(\hat{p}) = 5.08655 - 0.19400 X$ となっている。切片と年齢の回帰係数の推定値はいずれも1％水準で統計的に有意である。

　次にクロス集計された結果に対してロジスティック回帰分析を適用する方法を実行してみる。この場合、あらかじめクロス表形式の集計データを用意しておく必要がある。

X1	<- c(22, 27, 32, 37, 42, 47, 52)
Y1	<- c(1, 29, 56, 49, 59, 70, 76)
Y2	<- c(20, 25, 10, 3, 2, 3, 3)

d2	<- data.frame(X1, Y1, Y2)
d2						# データを確認

  X1 Y1 Y2
1 22  1 20
2 27 29 25
3 32 56 10
4 37 49  3
5 42 59  2
6 47 70  3
7 52 76  3

　結婚状態のカテゴリのうち、未婚（Y2）を1、既婚（Y1）を0として年齢×結婚状態のクロス表にロジスティック回帰モデルを当てはめる。

l2	<- glm(cbind(Y2,Y1)~ X1, data = d2, family = binomial)
summary(l2)

Call:
glm(formula = cbind(Y2, Y1) ~ X1, family = binomial, data = d2)

Deviance Residuals: 
       1         2         3         4         5         6         7  
 2.99032   0.01082  -1.87446  -1.31202  -0.47296   1.31558   2.35751  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.08655    0.78789   6.456 1.08e-10 ***
X1          -0.19400    0.02472  -7.848 4.23e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 130.614  on 6  degrees of freedom
Residual deviance:  21.689  on 5  degrees of freedom
AIC: 47.495

Number of Fisher Scoring iterations: 5

　ロジスティック回帰式の推定結果は個票データを分析したときと同じだが、Null devianceとResidual devianceの数値が先ほどの結果と異なっている。このResidual devianceの値が教科書の $G^2$ の値と一致するから、ロジスティック回帰式から推定される未婚者の期待度数と実際のセル度数との乖離の程度をあらわしているのだと思う。期待度数 $\hat{\mu}_{ij}$ と観測度数 $n_{ij}$ から尤度比統計量 $G^2$ を、次の式

$G^2 = 2 \sum_{\mathrm{cells}} n_{ij} \log_e \frac{n_{ij}}{\hat{\mu}}_{ij}$

によって計算することができるので、一応、計算結果を確認しておく。

P1	<- 1 - fitted(l2)
P2	<- fitted(l2)
M	<- Y1 + Y2
E1	<- P1 * M
E2	<- P2 * M
LR1	<- sum(Y1 * log(Y1 / E1))
LR2	<- sum(Y2 * log(Y2 / E2))
G2	<- 2 * (LR1 + LR2)

G2
[1] 21.6894

　後者のやり方だとモデルの適合度がただちにもとまるから便利。ロジスティック回帰分析について、もっときちんと勉強しないと。