Haggstrom, G.W. (1983) Logistic Regression and Discriminant Analysis by Ordinary Least Squares. Journal of Business & Economic Statistics, 1(3), 229-238.
勤務先の仕事の都合で、多数の多項ロジスティック回帰モデルのパラメータ推定値を、ダミー変数に対する線形回帰モデルのパラメータ推定値へと大急ぎで変換しないといけないという謎の用事ができてしまい(自分でもこいつ何言ってんだと思う)、慌てて読んだ奴。事情はちょっと書けないけれど、なぜそんなシュールな事態に陥ったものかと、正直、途方に暮れた。ところが読んでいる途中で、さらなる別の事情によって必要性が消滅したもので、この論文のほうは続きを読む気が失せてしまった。整理の都合上、読了としておく。
イントロダクション
母集団は\(m\)個の排反な下位母集団\(\pi_1, \ldots, \pi_m\)からなるとする。個体について\(q\)次元ベクトル\(x\)が測定されており、これに基づいて個体を下位母集団に分類したい。
所属先番号を\(y\)とする。事前確率を\(p_j = P(y=j)\)とする。\(\pi_j\)における\(x\)の条件付き密度を\(f_j (x)\)とすると、事後確率は$$ p(j | x) = P(y = j | x) = \frac{p_j f_j(x)}{\sum_k^m p_k f_k(x)}$$ である。事後確率が最大である所属先に所属させれば、正分類確率が最大になる。
というわけで、どうやるかはともかくとして、条件付き確率を単純な関数として表現し、そのパラメータを\(n\)個の観察事例でもって推定したい。
ロジスティック回帰モデルは$$ p(j | x) = \frac{\exp(\gamma_j + \delta^\top_j x)}{\sum_k^m \exp(\gamma_k + \delta^\top_k x)} $$である。\(m=2\)ならもっと簡単になって$$ p(1 | x) = \frac{1}{1+ \exp(-\alpha + \beta^\top x)} $$ ただし\(\alpha = \gamma_1 – \gamma_2, \beta = \delta_1 – \delta_2\)である。
識別のために、\(\sum_k \gamma_k = 0, \sum_k \delta_k = 0\)とか\(\gamma_m = 0, \delta_m = 0\)といった制約をかけたり、\((x, y)\)の同時分布を指定するパラメータを考えて、\(\gamma_j, \delta_j\)はその関数だという制約をかけたりする。
\((x, y)\)の同時分布がロジスティック回帰モデルを満たすのはどんなときか。
第\(m\)群を除く全ての群で\(\log\left( \frac{f_j(x)}{f_m(x)} \right)\)が\(x\)の線形関数であれば、同時分布はロジスティック回帰モデルを満たす。[ううむ… 残念ながら学力不足で、なぜこういえるのかわからない。すべての群で \(\log f_j(x)\) が \(x\)の線形関数であれば同時分布はロジスティック回帰モデルを満たす、というのならわかるんだけど、これはそれより弱い主張ですよね]
これを常に満たすような\(f_j(x)\)の定式化は複数あるが、そのひとつは$$ f_j(x) = c_j \exp \left( \frac{-(x-\mu_j)^\top \Sigma^{-1} (x-\mu_j)}{2} \right) \phi(x) $$ という定式化である。ただし、\(\Sigma\)は\(x\)の共分散行列で正則とする。[ううむ、ここも学力不足でよくわからない。これは同時分布がロジスティック回帰モデルを満たすことの十分条件であって必要条件ではない、という理解であっているだろうか]
各群で\(x\)が多変量正規分布に従い、その共分散行列が群間で共通に\(\Sigma\)だとしよう。このとき$$ f_j(x) = (2\pi)^{-q/2} |\Sigma|^{-1/2} \exp \left( \frac{-(x-\mu_j)^\top \Sigma^{-1} (x-\mu_j)}{2} \right)$$であるから、ロジスティック回帰モデルが満たされる。 ロジスティック回帰モデルと見比べると、$$\delta_j = \Sigma^{-1} \mu_j $$ $$ \gamma_j = \log p_j – \frac{\mu^\top_j \Sigma^{-1} \mu_j}{2} $$となっている。二値だったらもっと簡単で、$$\beta = \Sigma^{-1} (\mu_1 – \mu_2) $$ $$\alpha = \log(p_1/p_2) – \frac{\beta^\top (\mu_1 + \mu_2 )}{2} $$である。[ここは信じて写経した]
[細かいところわかんないけど、まあいいや、とにかく、各群の\(x\)が同じ共分散行列のMVNに従うという仮定があればロジスティック回帰モデルが満たされますってことね。なんだか不必要に強い仮定なような気がするけど、気にせずに先に進もう]
では、パラメータをどうやって推定するか。ふたつのケースを考えよう。
ケース1、\(y_i\) が確率変数である場合。尤度関数は、\(y_i = j\)であることを示すインデクスを\(v_{ji}\)として、$$ L = \prod_i^n \prod_j^m (p_j f_j(x_i))^{v_{ji}} = \prod_j^m p_j^{n_j} \prod_i^n (f_j(x_i))^{v_{ji}} $$ ケース2, \(y_i\)は定数で、\(m\)群の下位母集団から固定サイズ\(n_1, \ldots, n_m\)の標本をとってきた場合[ケース・コントロール研究みたいなものかな]。この場合は尤度関数から\(p_j\)が消える。
ケース2の場合、\(\mu_j\)のMLEは標本平均ベクトル$$ \hat{\mu}_j = \bar{x}_j = \frac{1}{n_j} \sum_i v_{ji} x_i$$であり、\(\Sigma\)のMLEは$$ \hat{\Sigma} = \frac{1}{n} \sum_i^n \sum_j^m v_{ji} (x_i-\bar{x}_j)(x_i – \bar{x}_j)^\top$$となることが知られている。ケース1で\(p_j\)が未知ならばそのMLEは\(\hat{p}_j = n_j / n\)である。というわけで、無事にパラメータのMLEが手に入りました。
ここからは、これらのパラメータMLEを、どうにかして線形モデルのOLS推定量として手に入れる方法を考えます。
二値の場合
ここでは\(\pi_i\)に属するかどうかをダミー変数ベクトル\(y\)で表す。データがあたかも線形モデル$$ y_i = \alpha + \beta^\top x_i + e_i$$を満たすかのように考えたときの\(\alpha, \beta\)のOLS推定値を\(a, b\)とする。なおこういうのを”intermediate least squares” (ILS)推定値という。残差二乗和を$$ SS_e = \sum_i^n (y_i – a – b^\top x_i)^2$$とする。
定理1。ロジスティック回帰係数のMLEとILSのあいだには以下の関係がある。\(K = n/SS_e\)として$$ \hat{\beta} = Kb$$ $$ \hat{\alpha} = \log(\hat{p}_1/\hat{p}_2) + K(a – 1/2) + n(1/n_1 – 1/n_2)/2$$ [証明。めんどくさいからスキップ]
[ILSから\(H:\beta = 0\)のF統計量を求める方法とか、\(\beta\)のSEやt統計量の求め方とか。だんだん「いいからロジスティック回帰のソフトを使えよ」という気分になってきちゃったので読んでないけど、とにかく求められるらしい]
多項の場合
先ほど示したように、ロジスティック回帰では $$ p(j | x) = \frac{\exp(\gamma_j + \delta^\top_j x)}{\sum_k^m \exp(\gamma_k + \delta^\top_k x)} $$ だが、分母と分子を\(\exp(\gamma_m + \delta^\top_m x)\)で割って、\(\alpha_j = \gamma_j – \gamma_m, \beta_j = \delta_j – \delta_m\)とすると、\(m\)を除く\(j\)について $$ p(j|x) = \exp(\alpha_j – \beta^\top_j x ) / \left( 1+ \sum_k^m \exp(\alpha_k + \beta^\top_k x) \right) $$ \(m\)については $$ p(m|x) = 1 – \sum_{j}^{m-1} p(j|x) $$ となる。でもって$$ \beta_j = \Sigma^{-1} (\mu_j – \mu_m)$$ $$ \alpha_j = \log (p_j / p_m) – \frac{\beta^\top_j (\mu_j + \mu_m)}{2}$$である。ここまではいいっすね？
[ここから約1ページにわたって説明があるんだけど、疲れちゃったのですっとばして…]
というわけで、
定理3。多項ロジスティック回帰のパラメータMLE\(\hat{\alpha}_j, \hat{\beta}_j\)は、ダミー変数についての線形モデルのOLS推定値を\(a_j, b_j\)、残差平方和を\(SS_j\)、\(K_j = n/SS_j\)として$$ \hat{\beta}_j = K_j b_j $$ $$ \hat{\alpha}_j = \log (\hat{p}_j / \hat{p}_m) + K_j (a_j – 1/2) + \frac{n(1/n_j – 1/n_m)}{2} $$ となる。
[t統計量はどうなるか、とか…]
ここまでで論文の2/3くらい。あとはめんどくさくなって目を通していないので、なにが書いてあるのか把握してない。分布が非正規だったらなにがどうなるかというような話かな？
正直、素直にロジスティック回帰のソフトを使えばいいじゃない… なんなのこの無駄な知識… いったいいつ活かすの…という気がしてきてしまったのだが、全然別の文脈で最近見かけた話で、どうしても線形確率モデル(二値変数を目的変数にとった線形モデル)を使いたいことがあるよね、しかしモデルによる予測値が0-1の範囲の外に出るので困るよね、これを防ぐ工夫を考えました、という説明の際にこの論文が引用されていた。あれって、パラメータ解釈は線形確率モデルでやるけど、予測の時だけロジスティック回帰モデルに変換するってことだったのかな。よくわからんが、まあとにかく、なんであれまるきり無駄なことというのはないのだろう。