ロジット変換に関して

ロジット変換に関して。

割合を目的変数とする厳密な解析では、ロジスティックモデルをあてはめたり、割合の発生率が小さな場合はポアソンモデルをあてはめたりする。目的変数が二項分布やポアソン分布に従うと仮定するからである。

だが、そういった目的変数に対しても、回帰モデルや分散分析モデルで解析を行う場合は多々あると思う。こういったモデルの方が分かりやすいだとか、ソフトウェアの都合上、こういったモデルの出力結果の方が充実している場合等考えられる。

目的変数は二項分布に従うと考えられるが、それにロジット変換を施して回帰/分散分析モデルで解析を行う場合に関してのメモを記す。

本業で工業系のデータを扱うことも多いので、不良品発生率を目的変数として書く。

不良品発生率の他に、化学品製造だと収率を目的変数とするような場合もある。

不良品発生率が真ん中あたりでばらついている場合（例えば50%）のデータだと変換を行わなくても問題がない場合も十分あると思う。
だが、不良品発生が0%付近に分布しているだとか、収率が100%付近に分布しているだとかの偏った分布の場合、そのまま回帰/分散分析モデルをあてはめると、残差の等分散性が全く満たされていなかったり、非線形な関係だったりする。こういった場合、変換を施すと改善される場合がある。

n個のサンプルを調査し、f個の不良品を発見したとする。

fは二項分布に従うと仮定する。
f ～ Bin( n , π)
π：母不良率

fの期待値は、
E( f ) = n π
fの分散は、
V( f ) = n π(1 - π)
となる。

不良品発生率は、p = f / n と定義できる。

よって、pの期待値は
E( p ) = π
Pの分散は、
V( p ) = π(1 - π) / n
と計算できる。

nがある程度大きいときは、中心極限定理によりpは正規分布に従うとする。
p ～ N( π, π(1 - π) / n )

---------- xのロジット変換 ----------
Logit(x) = ln{ x / (1 - x) }
オッズ：x / (1 - x)
対数オッズ：ln{ x / (1 - x) }
-----------------------------------------------

pにロジット変換を施す。
Logit(p) = ln{ p / (1 - p) }

デルタ法により、Logit(p)の期待値、分散は近似的に以下に計算できる。
E( Logit(p) ) = ln{ π/ (1 - π) }
V( Logit(p) ) = 1 / nπ(1 - π)

正規分布に従うとする。
Logit(p) ～ N( ln{ π/ (1 - π) }, 1 / nπ(1 - π) )

---------- デルタ法 ----------
確率変数xがE(x) =μ、V(x) =σ2のとき、微分可能な関数f()に対して以下が成り立つ。
y = f(x)
E(y) = f(μ)、V(y) = f ’(μ)^2 σ2
------------------------------------
デルタ法に関しては、こちらの資料など参考になると思う。
http://www012.upp.so-net.ne.jp/doi/biostat/CT39/delta.pdf

不良品発生率p をロジスティックモデルでモデル化する。
p = 1 / 1 + exp{ -(a+bx) }

簡単にするため、説明変数1つ（x）の単ロジスティックとする。

pにロジット変換を行うと直線になる。
Logit(p) = a + bx

それぞれ、
p ～ N( π, π(1 - π) / n )
Logit(p) ～ N( ln{ π/ (1 - π) }, 1 / nπ(1 - π) )
の分布に従い、

a、bを任意の値とし、不良品発生率pとそのロジット変換を縦に、横をxとして線を描き、その周りに「1.96×標準偏差」をプロットしてみる。

πの推定値はp。

n数は30とする。
n <- 30

切片、傾きをa = 0、b = 0.5 とする。
a <- 0
b <- 0.5

横軸ベクトル。
x <- -100:100 / 10

ロジスティックモデルでpをモデル化。
p <- 1 / ( 1 + exp(-(a+b*x)) )

pの「1.96×標準偏差」の上側下側。
pU <- p + 1.96 * sqrt( p*(1-p)/n )
pL <- p - 1.96 * sqrt( p*(1-p)/n )

plot( x, p, type="l", xlim=c(-10,10), ylim=c(0,1), main="n=30、両側に標準偏差の1.96倍", xlab="x", ylab="p" )
pの「1.96×標準偏差」の上側下側を重ね描き。par( new=TRUE )でplot()関数の重ね描きができるよう。
par( new=TRUE )
plot( x, pU, type="l", col="red", xlim=c(-10,10), ylim=c(0,1), xlab="", ylab="" )
par( new=TRUE )
plot( x, pL, type="l", col="red", xlim=c(-10,10), ylim=c(0,1), xlab="", ylab="" )
グリッド線を追加。
abline( h=(1:10)/10, v=(-5:5)*2, lty=3 )

pが0.5のときのpの分散が最も大きく、pが0や1に近いほどその分散は小さくなる。

pにロジット変換を行ったものにも同様のプロットを行う。
Logit(p)、直線になる。
logit_p <- a+b*x

Logit(p)の「1.96×標準偏差」の上側下側。
logit_pU <- logit_p + 1.96 * sqrt( 1/(n*p*(1-p)) )
logit_pL <- logit_p - 1.96 * sqrt( 1/(n*p*(1-p)) )

pの場合と同様にプロットを行う。
plot( x, logit_p, type="l", xlim=c(-10,10), ylim=c(-5,5), main="n=30、両側に標準偏差の1.96倍", xlab="x", ylab="Logit(p)" )
par( new=TRUE )
plot( x, logit_pU, type="l", col="red", xlim=c(-10,10), ylim=c(-5,5), xlab="", ylab="" )
par( new=TRUE )
plot( x, logit_pL, type="l", col="red", xlim=c(-10,10), ylim=c(-5,5), xlab="", ylab="" )
abline( h=-5:5, v=(-5:5)*2, lty=3 )

 

このように、中心付近でのLogit(p)の分散は小さく、端に行くほどLogit(p)の分散は大きくなる。

ちなみに、p = 1 / ( 1 + exp{ -Logit(p) } )なので、Logit(p)とpの対応は以下。
LogitP <- 5:-5
P <- 1 / ( 1 + exp( -LogitP ) )
data.frame(LogitP, P)

 　
   LogitP           P
1       5 0.993307149
2       4 0.982013790
3       3 0.952574127
4       2 0.880797078
5       1 0.731058579
6       0 0.500000000
7      -1 0.268941421
8      -2 0.119202922
9      -3 0.047425873
10     -4 0.017986210
11     -5 0.006692851

最小2乗法でLogit(p) = a + bxのa、bを推定することは理論的には正しくない。最小2乗法による推定は直線の周りに等分散を仮定している。

やったことはないのだが、重み付最小2乗法といったものが用いられると思う。

でも、実務では最小2乗法でやってしまうこともあると思う。データにもよるが、十分傾向は把握できる場合もある。でも、注意が必要。

最小2乗法は直線の周りでデータのばらつきを平等に扱うので、目的変数pの割合が0や1近く偏っているデータだと、そのばらつきを過小評価した推定になってしまう、ということだと思う。要するに、0や1近くのデータがいっぱい固まっている箇所をうまくモデルで捉えることができていないのでは、ということだと思う。

ちなみに、p = f / nのロジット変換（Logit(p) = ln{ p / (1 - p) }）だとpが0や1だと変換できないので、修正を施し、修正p（p*）をp* = ( f+0.5 ) / ( n+1 ) とし、そのロジット変換を経験ロジットと呼ぶそう。
Logit(p*) = ln{ p* / (1 – p*) } = ln{ ( f+0.5 ) / ( n-f+0.5 ) }

ロジット、経験ロジット変換に関して取り上げられている文献は以下。

入門 統計解析法

統計的方法のしくみ―正しく理解するための30の急所