【統計一口メモ 第32話】ロジスティック回帰って?

名古屋市立大学大学院医学研究科 非常勤講師 薬学博士 松本一彦

回帰分析というと、直線回帰、非直線回帰、重回帰などよく目にします。ロジスティック回帰も基礎から臨床まで広く使われていますが、どのような時に使うのでしょうか?

 ...

表1

 ...

投与量(x)は等比級数で6段階にして、それぞれn=10匹のラットに薬剤を投与し、毒性の発現がみられた例数(f)を観察した例です。このように、ロジスティック回帰は横軸が量的変数=投与量、縦軸が質的変数=%の場合に用いられる回帰分析です1)。

§1.直線回帰とロジスティック回帰の違い

表2 直線回帰

 ...

図1 直線回帰

 ...

図2 ロジスティック回帰

 ...

まず、図で確認すると一目瞭然でロジスティック回帰はS字状になっています。

この図を描くモデル式が直線回帰では y=a + bxで、ロジスティック回帰では後に述べる複雑な式になります。

§2.尤度と確率そして最尤法

確率(π)」はnが10でπ(%)が分かっているときの例数fを求めます。π=20%のときはf=2です。一方、「尤度」はnが10で例数(f)が分かっているときのπ(%)を求めます。

f=2のときはπ=20%です。例えばn=100でf=28が得られたとき、確率πを変化させて、πがいくつのときf=28という結果が出やすいかを考えます。その尤度が最大となるπを推定する方法を「最尤法」と呼びます。また、尤度の対数は「対数尤度」といいます。

fiが独立な分布にしたがうときに、全体の尤度は各fiの尤度の「積」で表されます。この「積」を最大化するということは、対数尤度の和を最大化することと同じです。そこで対数の符号を逆にする、すなわち、「-1をかける」と最大化が最小化という言葉に置き換わり、最小二乗法と同じ考え方に立つことができます。ここで注意しなければいけないことは、対数をとるときには常用対数(log)ではなく自然対数(ln)を用いることです。さらに、2倍しておくと、検定などに便利であるということで「-2ln(尤度)」の和が最小になるようにパラメータを決めることになります。-2ln(尤度)は尤度=likelihoodのLを使って-2Lと表します。

§3.ロジスティック曲線のあてはめ―ソルバーを使いこなす―

S字曲線(ロジステック曲線)のモデルは次の式で表されます。pの予測値(phat)を求める

ときに使います。

 ...

ロジスティック回帰は最尤法を用います。すなわち、-2Lを使うということです。エクセルでの計算にはソルバー関数を使って計算します。それでは、-2Lの求める手順を見ていきましょう。

表3 -2L値計算

 ...

手順1.パソコンにソルバーを前もってアドインしておく必要があります。ファイル →  オプション → アドイン → ソルバー → OKです。上のリボンの「データ」にある「分析」に入っていれば成功です。

手順2.まず、a、bに「適当な初期値を入れる」と教科書には書かれていますが、何を入れていいのかわかりません。そこでaは最終的にX50を求めることになるので用量xの中間あたりの数値2を入れておきます。bは全く想像がつかないので、fの最高値10としておきます。表3はソルバーを実施した後の結果を示しています。

手順3.pの予測値phatを求めます。トップのセル(ここでは[0.048]のセル)にモデル式('=1/(1+EXP(-(a+b*LOG(x))))にa、b、xを使って計算します。

実際には「'=1/(1+EXP(-(N$23+N$24*LOG($J16))))」のようにa、b、xのセルを選択します。$の位置に要注意です。ここでは、a=2、b=10、x=1.0のセルが選択され、0.048と計算されます。これを6行全部にコピーします。

手順4.-2L値を求めます。BINOMDIST関数を使って2項分布の尤度を求めて、それを-2倍します。式は =-2*LN(BINOMDIST( f, n, phat, FALSE ))を使って計算します。これを6行全部にコピーします。

手順5.Lセル(目的セル)にSUM関数を使って、-2L値の合計を求めます。

手順6.ソルバーでLを最小とするaとbを求めます。ソルバーを選択すると次のような画面が表示されます。

図3 ソルバーのパラメータ

 ...

「目的セルの設定」はLセルを選択し、「最小値」に〇、「変数セルの変更」は a、bを選択しますが、実際には、自動的に$を付けた値が表示されています。

「□制約のない変数を非負数にする」のチェックを外します。これを忘れないで。「解決」をクリックします。次の図のような「ソルバーの結果」が表示され「OK」で表3の結果が表示されます。

図4 ソルバーの結果

 ...

手順7.ここで、目的のX50値の計算をしましょう。今までのロジスティック曲線のphatのモデル式をX50値のモデル式に変更するだけでOKです。モデル式は次のようになります(エクセルではlogは自然対数LN、log10は常用対数)。

 ...

先のロジスティック曲線のaに代わってX50とします。初期値はaに2、bに10を入れて-2L値はロジスティック曲線と同じ式を使います。結果は表4のようになり、X50値=1.837が求まりました。

表4 W50値計算

 ...
§4.JMPでの解1)

JMPでは「非線形回帰」を用いて解析します。結果はエクセル解と同じです。

図5 非線形回帰のあてはめ

 ...
§5.Pharmaco Pharma(β版)での解2)

図6

 ...

表5

 ...

表6

 ...

いずれのソフトでもX50は1.837となりました。Pharmacoではさらに逆推定としてX値から投与量を推定することもできます。

図7

 ...

図8

 ...

以上のように、ロジスティック回帰はエクセルでもソルバーを使えば解析できることを説明しましたが、実務者は市販のソフトを駆使されると思います。ただ、JMPソフトでもかなり精通していないと取り組むことは難しいと思います。逆推定も同じです。

なお、ロジスティック回帰を使う目的の一つに、複数の薬剤の効力比を求めることがあると思います。それは、次回以降の一口メモで取り上げてみたいと思います。

以上

  • 1)芳賀敏郎「医薬品開発のための統計解析 第3部」サイエンティス社 2016年
  • 2)Pharmaco工房:Pharmaco Pharma(β版); https://pharmaco.club.com/