オッズを使ってロジスティック回帰してみる

ロジット・モデルと言うと個票に応用する二項選択モデルを最尤法などで推定することを思い浮かべる人が多いと思いますが、都道府県ごとの率など集計データに一般線形回帰をかける事で推定する事もできます。観測数不明の集計データしか取れない場合も多いので、こちらの方が重宝するかも知れません。試しに使ってみましょう。

1. ロジット・モデルの線形化

ロジットモデルは以下のような累積分布関数から事象の発生確率 $p$ を定めます。

$p = \frac{ e^{x \beta} }{ 1 + e^{x \beta} }$

$x$ は条件を表す説明変数のベクトル、 $\beta$ はその係数です。 $x \beta$ の値に関わらず $p$ は0から1の間に収まるので、1超や負のあり得ない確率は出てこなくなるのが利点です。
このときオッズ $p/(1-p)$ を整理すると、以下のようになります。

$\frac{p}{1-p} = e^{x \beta}$

対数を取ると以下のように、線形モデルになります。

$\log \frac{p}{1-p} = x \beta = \alpha + \beta_1 x_1 + \beta_2 x_2 + \cdots$

これに誤差項がつくと仮定して、重回帰をかけます。

2. 試しに推定してみる

手ごろなデータセットが『ブール代数分析による社会的カテゴリーの研究:「日本人」カテゴリー認識の分析』*1の中にあったので拝借してきて、推定してみます。日本国籍を持ち社会学関係科目を履修している関西の大学生205名が、日本国籍・日本在住・日本人血統・日本語話者の4属性から日本人か否かをどう判別するかアンケートを取ったもので、サンプリング・バイアスと血統属性のトートロジー感（親が日本国籍？）が気になるデータですが、推定自体は問題なくできます。

### データセット ###
# N:日本国籍，R:日本在住，B:血統，L:日本語話者，ANS:日本人とする率
df1 <- data.frame(
  N = c(1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0),
  R = c(1,1,1,0,0,0,0,0,1,1,1,1,0,0,0,0),
  B = c(1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0),
  L = c(1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0),
  ANS = c(0.995,0.922,0.922,0.659,0.995,0.756,0.82,0.502,0.639,0.4,0.278,0.059,0.493,0.18,0.049,0.015))
# オッズを計算する
attach(df1)
df1$ODDS <- df1$ANS/(1-df1$ANS)
detach(df1)
# 回帰分析をかける
r_ols <- lm(log(ODDS) ~ N + R + B + L,data=df1)
summary(r_ols) # 結果確認

指数をとらないと確率にならないので、推定結果だけでは分かりづらいですね。
親は日本人では無いものの、日本国籍がある日本語話者が日本人認定される確率を計算すると以下になります。

e <- exp(sum(coef(r_ols) * c(1, 1, 0, 0, 1)))
e/(1+e)

全ての説明変数が0を基準に限界効果を計算してみましょう。

e <- exp(sum(coef(r_ols) * c(1, 0, 0, 0, 0))) # c(...)の中の最初の1は切片項のため
coef(r_ols)[2:length(coef(r_ols))] * e/(1+e)^2 # e/(1+e)^2が密度関数になる

説明変数が不連続なので非直観的な数字になります。
これでは解釈しづらいので、0から1に変わったときの確率の計算もします。

l <- length(coef(r_ols))
m <- matrix(1, l-1, l)
m[,2:l] <- diag(l-1)
e1 <- exp(m %*% matrix(coef(r_ols), l, 1))
e2 <- exp(sum(coef(r_ols) * c(1, 0, 0, 0, 0)))
PCT_0to1 <- e1/(1+e1) - e2/(1+e2)
rownames(PCT_0to1) <- c("N", "R", "B", "L")
PCT_0to1

3. 多値選択モデル（Multinomial Logit Model）への対応

3値以上の多値選択モデル（Multinomial Logit Model）でも、複数の重回帰を行なう事で推定値が出ます。
選択肢の数を $J$ 、選ばれた選択肢を $m$ として累積密度関数は、

$p_m = \frac{e^{x \beta_m}}{1 + \sum^{J-1}_{j=1} e^{x \beta_j} }$

$p_0 = 1 - \sum^{J-1}_{j=1} p_j = \frac{1}{1 + \sum^{J-1}_{j=1} e^{x \beta_j} }$

となります。ここで $p_0$ は、直接推定されない参照カテゴリーが選ばれる確率です。 $\beta$ が選択肢ごとに変わることに注意してください。
オッズは、

$\frac{p_m}{p_0} = e^{x \beta_m}$

です。この式の対数をとって、 $m=1,\cdots,J$ の分だけ重回帰をかけると推定が終わりです。

となります。

4. 欠点

集計カテゴリーで差が出ないデータだと、集計することで説明変数と被説明変数の差異が小さくなります。この場合、検出力が露骨に落ちます。個票が使える場合は、なるべく個票を使いましょう。

おまけ：力技的に元データを復元

実はこのデータは、全ての説明変数の組み合わせの回答者数がそれぞれ205だと分かっているので、観測数16の集計データから、観測数16×206=3280の質的データに変換した上でglm関数で推定する事も可能です。

n <- 205
m <- 2^4 # 説明変数の数が4
df2 <- data.frame(
  N = numeric(n*m),
  R = numeric(n*m),
  B = numeric(n*m),
  L = numeric(n*m),
  ANS = numeric(n*m)
)
attach(df1)
for(i in 1:m){
  n1 <- round(n*ANS[i], 0)
  n2 <- n-n1
  b <- (i-1)*n + 1
  e <- i*n
  df2$N[b:e] <- rep(N[i], n)
  df2$R[b:e] <- rep(R[i], n)
  df2$B[b:e] <- rep(B[i], n)
  df2$L[b:e] <- rep(L[i], n)
  df2$ANS[b:e] <- c(rep(1, n1), rep(0, n2))
}
detach(df1)
r_glm <- glm(ANS ~ N + R + B + L, family = binomial(link = "logit"), data=df2)
summary(r_glm)