統計分析の手法は教科書で一度勉強したらそれでOK、もう二度と勉強しないでいい、ということはなくて、しばらくしたらまた忘れてしまう。だから、再び同じ手法を使う機会が出てきたらもう一度同じ本を勉強し直さないとならない。もっと要領よくやってる人もいるのかもしれないけれど、私は要領が悪いからそうやっている。統計学は、忘れては覚え忘れては覚えの「賽の河原」方式で覚えていくものだと思ってる。

今回の本は、前にロジスティック回帰を使うことがあって一応全部読んでいるのだけど、しばらく使う機会がなかったのでだいぶ忘れている。それで復習して、必要最低限のことは思い出せた。ただ、理屈ばかり思い出して、具体的に分析をどう進めればいいかというのはちょっと曖昧ではある。そこで、細かい理屈はすっ飛ばして、とにかくどうすればいいのか、どの数字を見ればいいのか、ということに着目して整理していきたい（もちろん整理の仕方は自己流で、自分の考えもちょこちょこ入れてる。文句があるなら本を買ってください。古本しかなくてバカ高いけど）。¹

ロジスティック回帰のモデルはややこしい

何らかの変数（説明変数）で別の変数（目的変数）のばらつきを説明したり予測したりする手法が回帰分析だ。だけど、目的変数がYes or Noのように2つの値しかとらないダミー変数である場合、普通の回帰分析を使ってもモデルがデータにうまく適合しない。というのは、「普通の回帰分析」とは正確には「線形回帰分析」というやつで、説明変数と目的変数が線形の関係じゃないと扱えないからだ²。

目的変数が連続変数ではなく、2つの値しかとらない場合、ロジスティック回帰分析という手法を使う。目的変数が連続変数か2値変数かなんて些末な問題のように見えるかもしれない。でも、目的変数が2つの値しかとらない、という制限を加えるだけで、モデルはかなり複雑なものになる。だからロジスティック回帰分析だけで本書のように600ページ近い教科書が必要になるのだ。³

ロジスティック回帰のモデルは以下のようなものだ（右端の方の括弧が上手く表示されなくてへんになってるけど、雰囲気だけわかればいい）。

$P(D =1 | X_{1}, ... , X_{k}) = 1 / (1 + \exp [-(\alpha + \Sigma \beta_{i}X_{i})$ ) ]

この時点で、なんだかやっかいなことになってきたなあ、と思う人もいるかもしれない。

ここで、左辺のPは確率だ。つまり、説明変数がある値をとったときに目的変数（D)が1になる確率を意味する。問題は右辺だ。なんだこのややこしい式は。

線形回帰のモデルだったらもっとシンプルにこうなる。

$Y = \alpha + \Sigma \beta_{i}X_{i}$

これは、Xが大きくなればYが大きくなる（あるいは小さくなる）という関係性を表している。係数のβが正なら正の関係だし、負なら負の関係だ。そして、係数の絶対値が大きければ、その説明変数が目的変数に与える影響は大きくなる。わかりやすい。

一方、ロジスティック回帰のモデルは、線形回帰のこの式を右辺の分母のところに入れている。その上、マイナスをかけて指数変換までしている。だから、ぱっと見、何がなんだかわからない。Xが大きくなったら何がどうなるのかとか、係数の意味はいったい何なのかとか、ぜんぜんわからん。

ロジスティック回帰では係数をみてもしょうがない。オッズ比を求めよ

そしてさらにやっかいなことに、実はロジスティック回帰では、P(の推定値)そのものは求められない。そのためには、「追跡研究」という研究デザインにしなくてはならない。それはたとえば、100人の被験者を対象にして、現在の説明変数が、数年後の目的変数にどう影響するかを見る、という研究デザインだ。お金をかけた研究ならそんなこともできるかもしれないけれど、普通はそんなことしないし、社会科学の場合ならなおさらやらない。だから、Pは求められないのだ。

そこで、ロジスティック回帰では確率の代わりに「オッズ」という考え方を使う。オッズとは次のようなものだ。

$オッズ = \frac{P}{1-P}$

確率とは表現が違うけれど、まあ、確率に似たものだと考えていいだろう。確率（P）が高ければオッズは大きくなるし、確率（P）が低ければオッズは小さくなる。ただし、確率とちがって、最大値は1ではない。理論的には無限大まで大きくなることができる。

そして、ロジスティック回帰で注目するのはこのオッズによって作られた「オッズ比」だ。

$オッズ比（OR）　= \frac{注目している説明変数が1のときに目的変数が1となるオッズ}{注目している説明変数が0のときに目的変数が1となるオッズ}$

「説明変数が1のときに」のところは、書き方が正確ではない。ここの値が変わる場合もある。とりあえず雰囲気を示すためにザックリ書いているだけだ。

このオッズ比を求めるのが、ロジスティック回帰のひとつの目標だ。オッズ比は確率（P）そのものではない。ある説明変数が大きな値を取ったときに、目的変数が「1」になる「なりやすさ」を示している。説明変数が大きな値を取ったときに、目的変数が「1」になるオッズが大きくなるとしたら、オッズ比は大きくなるだろう。逆に、説明変数が大きな値を取ったときに、目的変数が「0」になるオッズが大きくなるとしたら、オッズ比は小さくなるだろう。だから、目的変数が「1」になりやすいとき、オッズ比は大きくなるといえる。確率（P）そのものは求められないので、オッズ比を使うことで、その説明変数が目的変数にどんな影響を及ぼすかを見ようというわけだ。

では、オッズ比はどうやれば求められるのか？　交互作用項があるときはちょっとやっかいだけど、基本的には次の式で求められる（ただし、コードの仕方で変わってきます。ここでは、説明変数が0, 1でコードされていると仮定している）。

$オッズ比(OR) = \exp(\beta_{i})$

つまり、説明変数Xiが目的変数にどんな影響を及ぼすか、というオッズ比をみたいのなら、上の式にその説明変数の係数を代入してやればいいのだ。

逆に、ロジスティック回帰の場合、係数であるβそのものには意味がないことに注意が必要だ。線形回帰の場合だと係数を見れば「なるほど、体重が1キログラム増えるごとに、かけっこのタイムは0.1秒伸びるんだな」みたいな判断ができる。だけど、ロジスティック回帰の場合は、係数をいったん指数変換してオッズ比にしてやらないと、意味のある解釈ができないのだ。

検定の前に…最尤法って？

検定の前に、まず、尤度という考えを知っておく必要がある。

線形回帰では最小自乗法という方法を使って、モデルをデータに適合させる。これは、モデルによる予測値と観測データのズレの二乗和がもっとも小さくなるようにモデルの切片と係数を決めるという方法だ。

しかしロジスティック回帰では最小自乗法は使えない。なぜ使えないかという理由は省く。代わりに使われるのが最尤法だ。これは、尤度関数が最大になるようにモデルをデータに適合させる方法だ。

尤度関数とは、実際に観測されたデータについて、そうしたデータが観測されるもっともらしさ（尤もらしさ）を表す関数だ。なんのこっちゃって感じなので、シンプルな例でイメージを見てみよう。これは、コイン投げを100回して、表が75回というデータが観測された状況だ。Pはわからない。この場合の尤度関数は次のようになる。

$Pr(表=75 | n =100, P) = 定数 \times P^{75} \times (1-P)^{100 - 75}$

つまり、Pの同時確率で表せるということだ。最尤法では、この尤度関数が最大になるようにPを決める。この場合、Pは0.75になる。つまり、このコインは歪んでいて、表が出る確率が75％もある、というわけだ。「こうした結果が出るのが尤もらしいと思える度合い（尤度）を最大にするように確率（P）を推定する」のが最尤法だ。

ロジスティックモデルで同じことをやろう。この場合、Pではなく、ロジスティックモデルの切片や係数をいろいろいじって、尤度関数を最大化することになる。もちろんそんな計算はコンピュータがやってくれる。また、ロジスティックモデルを使った尤度関数はやたら複雑だけど、そんなのいちいち自分で作らなくてもいい。最尤法の基本的なイメージがわかっていればそれで十分だ。

検定の方法1：尤度比検定

最尤法を使ってロジスティックモデルの切片と係数を推定するということは、統計ソフトでロジスティック回帰を実行すると最尤値も出てくるということだ。最尤値そのものには意味はない。でも、モデル同士の最尤値を比較することで、係数の検定をすることができる。

たとえばロジスティック回帰で、目的変数は同じだけど、説明変数が異なる次の2つのモデルをつくり、最尤法でデータに適合させてみたとしよう。

モデル1
- 説明変数はa, b, c, d, eの5つ
モデル2
- 説明変数はa, b, cの3つ

ここで、モデル1とモデル2の最尤値を比較してみよう。もし、両者に有意な差があるとしたら、モデル1から説明変数d, eを取り除くべきではない。取り除くと、最尤値が有意に減ってしまうからだ（モデル1はモデル2を包含する関係にある。小さいモデルの方を包含する大きいモデルの方が最尤値は大きくなる）。逆に、有意な差がないのなら、d, eは取り除いてもいい。これを「尤度比検定」という。

なお、細かい話だけど、尤度比検定は、最尤値そのものを比較するのではない。尤度比を使ってLR統計量というものを求めるのだ。

$LR = -2 \ln(小さい方のモデルの最尤値) - (-2 \ln(大きい方のモデルの最尤値))$

LRは、カイ二乗分布で近似できる。自由度は0と設定したパラメーターの数、つまり、小さい方のモデルで取り除いた説明変数の数だ。

なお、さっき書いたように「大きい方のモデルの最尤値」の方が「小さい方のモデルの最尤値」よりも大きくなるけど、LRの式ではそれぞれの最尤値の対数にマイナスをかけているので、大小関係は反対になる。だから、LRは必ず0以上となる。

検定の方法2：Wald検定

説明変数1つずつについて検定したいときはWald検定をしよう。Wald検定では、次のようにZという統計量を求める。

$Z = \frac{\hat{\beta}}{S_{{\hat{\beta}}}}$

ここで、分母は推測された係数βの標準誤差を示す。このZを二乗した統計量は自由度1のカイ二乗分布に従うので検定ができる。

だったらWald検定だけやればいいじゃないの？　なんで尤度比検定なんて面倒なことをするの？　そう思われるかもしれない。

実は、サンプル数が多いときは問題ないのだけど、サンプル数が中規模以下だと、次のようになってしまうのだ。

$LR \not\approx Z^{2}$

つまり、尤度比検定とWald検定で、統計量が変わってしまうということだ。こういう場合、LRの方が優れていることが明らかにされている。つまり、サンプル数が少ないときは、Wald検定よりも尤度比検定の方をやるべきなのだ。

区間推定

オッズ比を求め、有意かどうかを判定したら、95％信頼区間もちゃんと求めよう。「オッズ比が有意だからいいじゃん」と思考停止になるのではなく、その値の精度を理解するためには信頼区間が必要だ。有意であっても信頼区間がやたらと広かったら困る。たとえば点推定されたオッズ比が3.0であっても、その信頼区間が1.2～5.9だったりしたらオッズ比がどれくらいの大きさだと考えればいいのかよくわからなくなるだろう。社会科学なら別に良いかもしれないけれど、疫学みたいに人の命が関わる分野でこんなに推定がガバガバだったら結果を現場にどう生かせばいいか判断できない。だから信頼区間を求める事は重要なのだ。

95％信頼区間は次のように求める。

$オッズ比の95％信頼区間 = \exp [\beta_{i} \pm (Zのパーセント点 \times 標準誤差の推定値$ ]

で、「Zのパーセント点」は、95％信頼区間の場合は1.96となる。

ただし、これは説明変数をどうコードするかで変わってくるので注意。たとえば、(0, 1)ではなく(-1, 1)とコードされていたら、1と-1の差は2なので、上の式のexpの中身は2倍になる。

あと、交互作用があるケースは、標準誤差の部分に他の係数の分散・共分散が影響してくるので、非常に複雑になります。

モデルの適合度をどう評価するか：逸脱度とHL統計量

さて、モデルがデータにどの程度適合しているうのか、どう評価すればいいのだろう？　つまり、回帰分析でいうところの「決定係数」に当たる指標はどんなものなのか、という話だ。

「逸脱度」という指標がある。これは、評価対象のモデルの最尤値と、飽和モデルの最尤値を比較した指標だ。こいつはカイ二乗統計量に近似できるので、検定ができる。

$逸脱度 = -2 \ln (評価対象モデルの最尤値 / 飽和モデルの最尤値)$

でもその前に、「飽和モデル」ってなに？　これは、結果を100％完璧に予想できるモデルのことだ。たとえば、各被験者ごとにその被験者であるかどうかを示すダミー変数をつくり、それを全部モデルに説明変数としてぶちこむ。たとえばAさんを表すXaというダミー変数は、Aさんに関しては1で、他のBさんやCさんに関しては0という値を取る。目的変数が「ガンである」であって、Aさんがガンだとしたら（「ガンである」=1）、Xaの係数を1にすればいい。そうすると、Xaというダミー変数によって、Aさんの目的変数の値を完全に予測できる。これを、すべての被験者についてつくるのだ。そうすれば、結果を100％完璧に予想できることになる。

そして、もちろんこんなモデルは無意味だ。普通、モデルをつくるのは、より少ない変数で結果を予測したり説明したりするためだ。被験者の数だけ説明変数を投入しても、興味のある現象に関して何も知見を得ることができない。だけど、「完全なモデル」という意味では一種のベースラインとして使える。この飽和モデルをベースラインとして、そこから大きく逸脱していないモデルなら、そのモデルはデータによく適合している、といえるのだ。

だから上記の逸脱度で検定をするときは、有意になってはいけない。有意になったら、それは飽和モデルから大きく逸脱しているということであり、適合度が低いということなのだ。

ところが、逸脱度には問題がある。共変量パターン別のサンプル数が1に近い状況になると、カイ二乗統計量に近似できなくなるのだ。「共変量パターン」というのは、たとえば「年齢が50歳で、性別が男で、タバコを吸っていて…」といった、各変数の値のパターンのことだ。もし説明変数に連続変数が含まれていると、共変量パターン別のサンプル数が1になりやすくなる。なぜなら、たとえば年齢のような連続変数は「51歳、52歳、53歳…」と非常に多くの値を取り得るからだ。そのため、共変量パターンが同じになるサンプルが少なくなってしまうのだ。

そこで、問題のある逸脱度よりもHL統計量を使った方がいい。HL統計量の求め方は面倒だが、全部ソフトがやってくれるので覚えなくていい。基本的な考え方は逸脱度と同じで、有意でない方が適合しているといえる。ちなみに、HL統計量もカイ二乗統計量に近似できる（ただし自由度の求め方はちょっとややこしい）。

目的変数が3つ以上の値を取るときは？

目的変数が2値ではなく、3つとか4つとかの値をとることがある。こういうときは多項ロジスティック回帰を使おう。

たとえば、目的変数が(0, 1, 2)の3つの値を取るとする。この場合、これらの値のうちのどれか（普通は0）をベースラインに設定して、それと1または2という別のカテゴリーとを比較する、というやり方をする。

「比較」というのは、具体的には次にように、オッズ（のようなもの）を求めるのだ。

$0と1を比較したオッズ（のようなもの） = \frac{P(D=1)}{P(D=0)}$

$0と2を比較したオッズ（のようなもの） = \frac{P(D=2)}{P(D=0)}$

これらはオッズそのものではない。下がオッズの正しい式だ。

$オッズ = \frac{P}{1-P}$

だから、オッズを求めるのなら、本当は次の式でないといけない。

$1以外と1を比較したオッズ = \frac{P(D=1)}{P(D=0) + P(D=2)}$

だけど、これだと「0と1」の比較ではなく、「1以外と1の比較」になってしまっている。だから、オッズを使ってしまうと主旨が変わってくるのだ。

多項ロジスティック回帰の場合、オッズ比を求めるなどの計算に使うのは「オッズ」ではなく「オッズ（のようなもの）」の方だ。具体的には次の式を使う（例によって、「説明変数が～のときに」の部分は、ケースによって変わってくる）。

$オッズ比 = \frac{注目している説明変数が1のときに目的変数がgとなるオッズのようなもの}{注目している説明変数が0のときに目的変数がgとなるオッズのようなもの}$

ここでgは、目的変数のカテゴリーの値だ。この場合、1か2の値を取る。

オッズ比がこのように求められるということは、gの値がなんであるかによって、オッズ比が変わってくるということになる。つまり0と1を比較したときのオッズ比と、0と2を比較したときのオッズ比は異なるのだ。多項ロジスティック回帰で気をつけるべきことはここだけだ。

別の例でもう少し見てみよう。たとえば、目的変数のカテゴリーが(0, 1, 2, 3)の4つだとする。そして、0をベースラインの設定する。このとき、オッズ比は目的変数のカテゴリーの比較の組み合わせごとに求めなければならない。つまり、次の3パターンのオッズ比を求めなければならないということだ。

目的変数＝0と目的変数＝1を比較したときのオッズ比
目的変数＝0と目的変数＝2を比較したときのオッズ比
目的変数＝0と目的変数＝3を比較したときのオッズ比

検定したり信頼区間を求めたり、といったことは2値ロジスティック回帰のときと同様にやることができる。でも、「目的変数のカテゴリーのどれとどれを比較するか」ごとにオッズ比を求めなければならない、というのが面倒なところだ。ソフトが勝手に計算して出してくれるけど、それを論文掲載用の表に直すのはなかなかめんどくさいかもしれない。

ところで、こんな風に比較の組み合わせごとにオッズ比を求めないとならないのなら、2値ロジスティック回帰を繰り返し実行しても同ようなものじゃないか？　と思われるかもしれない。つまり、こんな風にやるのだ。

目的変数を「1」と「1以外」にコードし直して2値ロジスティック回帰を実行してオッズ比を求める
目的変数を「2」と「2以外」にコードし直して2値ロジスティック回帰を実行してオッズ比を求める
目的変数を「3」と「3以外」にコードし直して2値ロジスティック回帰を実行してオッズ比を求める

これでも確かに同じようにオッズ比は求められる。でも、結果が多項ロジスティック回帰と変わってくることがある。というのは、多項ロジスティック回帰でモデルをデータに適合させるときの尤度関数は、目的変数のカテゴリーが3つ以上であることを踏まえたものになっているからだ（尤度関数の式の導出は簡単だけど、書くのがめんどくさいので書かない）。一方、2値ロジスティック回帰を繰り返し実行するやり方だと、目的変数のカテゴリーが3つ以上あるという情報が無視されてしまう。だから、目的変数のカテゴリーが3つ以上なのなら、2値ロジスティック回帰を繰り返すのではなく、素直に多項ロジスティック回帰を行うべきなのだ。⁴

なお、本書ではモデル構築のプロセスについてかなりページ数を割いているけど、全部カットしています。というのは、内容が交互作用や交絡をどう扱うか、というものだから。この本は疫学の研究者を対象にしたもので、疫学だと年齢とか性別とか喫煙状況とかが交絡になってくるのだと思うけど、自分の研究の場合、交絡って発想がそもそも必要なのかよくわからない。年齢や性別は単なる交絡じゃなくて、それら自体が興味のある変数だ。で、交絡について考えない場合、本書で説明されているモデル構築プロセスはほとんど不要になってしまう。いや、もしかしたら交絡を考えないようなモデルをつくること自体が間違っているのかもしれないけど…。でも、そもそも関連する既存研究がほとんどないことやってるから、何を交絡と考えればいいかもよくわからない。もうちょっと勉強していったら「やっぱりあのモデルづくりは不適切だった」となるかもしれないけど、とりあえず今回はモデル構築について考えないことにする。↩
「線形」ってのは、X軸に説明変数、Y軸に目的変数を設定してプロットしたら、右上がりや右下がりの直線的な関係になるということ。ちなみに、Y=X²みたいな関係でも線形と判断できる。なぜなら、X²＝Zという風に置き換えてしまえば二乗の部分が消えるからだ。こういう風な置き換えをしても直線的な関係にならない関係のことを「非線形」という。そして、今回扱うロジスティック回帰はまさにこの非線形な関係を扱うものだ。↩
それはこの教科書が馬鹿丁寧だから、というのもあるのだけど、類書の『データ解析のためのロジスティック回帰分析』でも500ページ近くある。洋書の教科書も1冊持ってるけど、これも350ページ近くあるよ。『Rによる多変量解析入門』でロジスティック回帰を説明する章は、理屈の説明をほとんど省いているので20ページくらいしかないけれど、知識がない人が読んでも何が何だかわからないと思う。↩
なお、目的変数が順序変数の場合は順序ロジスティック回帰というのをやるといい。ただ、今回は分析に使う予定がないので端折ります。↩