シグモイド曲線とは

たとえば、塩水の濃度を舌で感じる強さは、極端に薄い濃度の範囲では塩分濃度を上げても塩味はあまり感じないけれど、ある濃度あたりから濃さの応じて比例した塩分の強さを感じるようになります。あまり濃くなってくると今度は濃度を増やしてもただ濃い濃度にしか感じない領域に入ります。

薬などの投薬効果も同じようにシグモイド曲線になることが知られています。我々の周囲にはシグモイド曲線を描く多くの事例があります。増加が頭打ちになる曲線はたいがいこれです。商品の売上数のライフサイクル前半とかも入るかな。

シグモイド曲線は、ジャンルによって使用上の目的が違うので、いろいろの名前で呼ばれます。成長曲線、需要曲線、微生物増殖曲線、死亡率曲線、サイフサイクル曲線などなど。あ！、ラだ。

同じＳ字でも、適用する関数式がいくつかあり、分野によって好まれているやり方があります。近似式の対象が、何らかの方法で本当にその関数になると証明されていないかぎりどの関数を近似式に選ぼうと自由です。しかし違う方法を適用すると、その分野では自分だけなぜその関数を適用したかを説明せざるをえないハメに陥ります。

関数式は次の３種がよく使われています。シグモイド形状を何とか直線回帰の係数計算にもちこむためにそれぞれ固有のデータ変換方式があります。はじめの二つはＳ字形状が回転対称を前提にしていますが、ゴンペルツは自由です。それぞれ制約なりルールがあります。

累積正規分布曲線　：プロビット変換
ロジスティック曲線：ロジット変換
ゴンペルツ曲線　　：対数化した差分

とりあえず　プロビット変換の説明

死亡率 p からプロビット変換するには、エクセルなら関数に準備されているよ。同じマウス数に対して薬品の投与量を2倍、4倍、8倍・・・といった指数比で水準を変えてそのときの死亡率または有効率を p とします。投与量の単位はいろいろです。

fx ボタンを押す>統計>NORMSINV（標準正規累積分布の逆関数）>OK>ｐのセル参照

各水準の供試マウスの数が同じ時には通常の回帰分析が使われる。えられた回帰式の p = 0.5 (半数死亡または半数有効)のときの投与量が LD50(50% Lethal Dose) または EC50(50% Effective Concentration) 値となる。

元の確率 p に戻すには 5 を引いておいて NORMSDIST 関数かな。標準化の系で考えるのがイヤという人は NORMINV 関数を使えばいいよ。5 ×標準偏差値を足すことになるけど。面倒でしょ?

 ゴンペルツ曲線の求め方

 ゴンペルツ 攻略作戦

我々の周囲で完全に回転対称になる曲線は限られているから、ここではその中でも比較的適用範囲が広く、ここで勉強した範囲に近い方法として、ゴンペルツ曲線（Gompertz curve）の当てはめ方を学ぶことにしましょう。
ただし、複雑な形状の曲線なので、これまで勉強した回帰計算から直接最終の式の係数がえられるわけではありません。次の作戦で進みます。

1. 次の（１）式がゴンペルツ式です。最終目標は、この式の係数または定数　K、a、b の値を求めてこのゴンペルツ式を数値の入った式にすることです。

   

2. まず最初は、 K、b の値を回帰計算から求めます。たとえば、去年の売上と今年の売上の差の数値を △売上 で表すやり方で、（１）式は（２）式に変換できます。この式をよく見ると、独立変数 log y と、目的変数 △log y との単回帰式になっています。  証明を開く --- 再クリックで閉じます

   ---

   （証明）

    

   （１）式は両辺の常用対数をとると、次のように変形されます。

   log y₁ = log K ＋ b^x･log a ・・・・・・（４）

   ここで、ｘ の変化量が 1 ずつ変わるデータを考えます（これが基準変更の根拠です）。そうするとｘが微小量 １だけ増えた時の（４）式は （５）式となります。

   log y₂ = log K ＋ b^x＋1･log a ・・・・・（５）

   両辺についてそれぞれ差（５）−（４）をとります。左辺を△log y とおくと

   △log y = log y₂ − log y₁ = b^x･log a ( b − 1 )

   （４）式を代入して x , a を消去し整理すると

   △log y = ( 1 − b )log K − ( 1 − b )log y₁ ・・・・・（６）

   （６）式は（２）式と同じ。証明済！

   b も K も定数だから、（６）式は△log y を目的変数としてlog y を従属変数とする回帰式になっていることがわかります。ｘ の変化量が 1　ずつの場合のゴンペルツ曲線なら回帰計算で係数を求めることができることが証明されました。

   変数が等間隔での変化量ならば、変数変換して変数を１ずつ変化する量に変換しておけばこの方法から近似式を得ることができますね。

    

   （混乱しないように）

   （６）式から b , K を求めるということは、最初の従属変数 x は１ずつ変化することを前提にした近似式であることになります。どのような従属変数 x データを持ってこようが関係ないところで、 x が 1 ずつ変化しているときのシグモイド形状として目的変数だけで計算されています。
   しかしあなたが知りたいのは観測した変数と目的変数 y との関係ですね。最初のシグモイド曲線を引き伸ばしたり位置をずらしたりしてあわせる作業が残っています。

   だからあなたは、１ずつ変化する x と、あなたの観測した等間隔の従属変数との関係式を求めて代入することになります。
   でもこれだけでは、シグモイドの形状が決まるだけです。まだ、あなたが観測したシグモイド曲線の位置 a は決まっていません。シグモイド曲線の位置は x = 0 の位置に対応する目的変数の値から方程式で算出されます。 x = 0 の位置は実績グラフを描き、変曲点を何らかの方法で推定するわけです。タコの方法は従属変数も観測データがばらつくことを予想して、変曲点近傍の3点のデータは直線であることを前提にその3点の目的変数の平均値で推定しています。

   （１）式に K　と b との数値を代入した式からエクセルなら、ゴールシークで求めることができます。

   直線になるはずのあたりのデータがあまりにばらついているのなら、手でシグモイドカーブを書いて変曲点を目で見て推定するという原始的なやり方が実務的な近似には有効です。発表には不向きですが。

   ---

   　

   △log y = ( 1 − b )log K − ( 1 − b )log y 　･････　（２）

    log y = X ; △log y = Y ; −( 1−b ) = A ; ( 1−b )log K = B とおくとこの式はつぎのようになる。

   Y = AX ＋ B

   第５章で説明したように、回帰計算できるだろ？
   回帰計算結果の偏回帰係数 A が数値で得られるから、 - (1 - b) の値が求まりますね。この値から簡単な計算で b の値が決定します。

   （２）式の誘導経過から、差 △log y の数値をいれる欄は、log y₂ - log y₁ のときは y₁ 側（上の欄）に入ります。

3. 同時に回帰計算出力結果の定数項から B の値が得られます。 B = ( 1 − b )log K です。 b はすでに決まっているから、K　の値も簡単な計算で決定できます。これで（１）式の K、b が決定できました。

4. さて、残る a です。ちょっと工夫が要ります。
   x が 0 の場合、（１）式の b^x は b⁰ = 1 となります。X = 0 の時の対応するデータを Y₀ とすると式が
   　Y₀ = K a となり　Y₀ 、 K が既知のため a を計算できます。 x が０の場合というのは変曲点にあたります。だから変曲点にあたる y を推定できれば y₀ になります。

   （注） 下の計算例で基準変更の欄があります。そこで 2001 年に数値が０となっています。このことを説明しています。

5. 指数計算ですから、少しの誤差が大きく結果に響くので、 x = 0 の位置をグラフの直線的に増加している中央（変曲点の近く）にして、そのときの観測データ y、K の値から （１）式より a 値を求めます。

6. で、ちょっとまった。データにはばらつきがあるから、1点だけで、 y 値を推定するのは危険です。だからゼロポイントに直線部分を選んでいるのです。x が０前後の3点の 対応する y 値の平均値を、ゼロポイントの y 値とします。

7. ここまでは、変数 x （ 1 ずつ変化する）との近似式になっていますから、実際のグラフは元の観測値にもどすことが必要です。たとえば、1995年を x = 0 としたなら、1996年は 1995 + 1 という換算式をとおして 　x = 1 となります。
   今後の y 値を予測したいときにはこの例の場合、（１）式の x → ( x − 1995 )と置き換えた式で算出することになります。

 じゃあ、計算事例にいこか

 ≪あなたの課題≫

1. 今後の最大売上はどれくらいに見積もられるか。
2. 2007年度の売上見込みはどれくらいに見積もればいいかを知りたい。

 ≪計算の手順≫

1. 独立変数になる log y の計算欄を作ります。1998年の例では、エクセルなら関数をもちいて =log10(13) と計算させればいいです。独立変数の黄色の列をつくります。

2. 次に目的変数となる △log y を計算します。1998年にある数値は、log y₁₉₉₉ − log y₁₉₉₈ で計算されています。
   過去側の年度の欄に差の数値が書き込まれていることに注意してください。

3. 黄色の欄で回帰計算をします、神田ソフトの場合目的変数をいちばん左に置くルールなので、コピーする前に 列の順番を置き換える必要があります。最初から神田方式に列をならべて計算するほうが簡単ですよね。また特に神田ソフトである必要は無く、エクセルの回帰計算で行ってもかまいません。

4. 回帰計算が終了すると、偏回帰係数と定数項の数値が出力されてきます。表の右側に log y の係数、定数項とあるのがその数値です。

5. （２）式の log y の係数は、 − ( 1 − b ) ですね。定数項の下にある二番目のブルーの段の式が成立します。式を解けば、 b の値が決まります。

6. 定数項からも同様の方法で、 K の値が決定されます。

   b = 0.429
   K = 848007

   K の値から今後の最大販売高は、 848,007 千円と見積もられます。

   （注）変曲点以降のＳ字形状の上のほうの実測値がない中では、推定精度はよくないですから、報告書にはあまり断定的に書かないように。散布図の中に出来上がった近似曲線を描いてみれば判断しやすいでしょう。

7. では、 a の値を求めてみましょう。（１）式に K、b を代入して（３）式がえられます。

   

   変曲点は、グラフから見ると2001年あたりにありますから、基準年を2001年にしましょう。上の表には基準変更欄 x に計算値が記載されています。2000年〜2002年の y 値を平均して2001年の推定値 y₀ を算出します。

   ( 110074 + 352934 + 582474 ) / 3 = 348494
   y₀ = 348494

   （注） うまくグラフに近似式がフィットしなければ、手書きで描いたシグモイド曲線の変曲点 x = 0 のときの y 値をグラフから読み取り、 y₀ としてみてください。
   

8. b⁰ = 1 から（１）式は（４）式になるのでこれに代入します。

   a = ₁₀ ( log y₀ - log K ) ･･････　（４）

    a = 10 ( log 348494 - log 848007 )
      = 0.410956513
   

   a = y₀ ／　K　で計算してもよいです。

9. これで、求める数値 K、a、b を全て得ることができました♪　その式は次のとおりですね。

   

   観測値の曲線（青線）と、ゴンペルツ近似曲線（ピンク線）はほとんど重なるほどに一致した。

10. 最後の課題、2007年度の需要予測値は、上式の x 値に 2007 を代入して計算します。

    y₂₀₀₇ = 843334 千円

11. はい,お疲れ様でした。君の解析力の巾がすこしひろがったかな。
    ゴンペルツ曲線の式は、色々あるから私の知っている式ではないと力説しないように。一番大事なことは観測値の曲線をいちばんよく近似させることだからね。

  y がゼロから始まらないときは？　LD50 はどう計算するの？

y 値が、 0 から始まらない場合には、次のようにします。 x 値が 0 の場合の y の観測値（通常はブランク値といったり、コントロール値といったりします）を y_min とします。
技術的に y_min 値が決定できない場合には一般にシグモイド曲線は近似できません。
観測値がシグモイド曲線として単調増加するならば y > y_min です。

y' = y - y_min　･･････（５）

  と変数変換すれば、 y' 値は上表の売上 y の欄に該当するのでこの欄に数値を入れます。

観測データのなかの最小値を、機械的に y_min 値として扱うと、それが比例的に上昇している部分である場合が生じます。 y_min 値は、シグモイド曲線の下側の水平部分の値である必要があるので、こういったときにはシグモイド近似はひどくあわないものになってしまいます。

（５）式で変換したときには、算出されてくる K 値は次の(6)式の値になっています。報告書のなかで K 値で議論しているのか、 y_max 値で議論しているのかはしっかりと意識しておく必要があります。観測座標の中に推定データのグラフを並べて表示すると間違いを発見しやすいです。

K = y_max - y_min ･･････　（６）

  LD50 などの計算をしたいときには、 y_min 値がゼロから始まるケースに換算した系で、 y = K/2 と置いたときの（１）式を移項して変換した（７）式から算出します。算出される x 値は、基準変更した数値になっていますので、あなたが設定した換算式で元に戻す必要があります。

  どうも、ゴンペルツは皆が知らないから流行じゃないし、ロジスティック回帰分析をやってみた いというひとが多そう。主成分分析の流れで、関数の当てはめをおまけにつけたんだけど、どうしてもロジスティックがいいという人もいる。

ん〜っ。私としては、なぜか二流扱いされている主成分分析のすごさをみんなに知ってもらいたかったのだけど、お客様には負けるなぁ。残念ながら最尤法を扱えるソフトをもってないのよ。お金ないしね。
仕方ない。計算事例は諦めてね。とりあえず、ロジスティック回帰分析の原理と手順を示しておけば、他のもっと数学的に解説してある文献も読みこなせるようになるかな。

でも、ほんとにすごい主成分分析もやってみてね。主成分分析は情報を集約する方法などという昔の使い方じゃない最新の方法をアップしているんだからね。お約束だよ。

ロジスティック回帰分析をなんに使うか。ロジスティック回帰はお医者さんが好きだから、病気診断の判断ツールとして使うことが多い。病気の発生確率を関係する他の因子（コレステロール値等）から計算する。などよく使われているな。

ロジスティック回帰分析は、発生する・しないなどという２分した事象の確率とその他の因子群との関係を回帰式にする場合に使われます。ロジスティック回帰分析では、分類型のデータを集計して発生確率として扱っていることに注目してね。
確率を計算するには、多くの同一条件での発生数をカウントすることが必要になります。ひとつひとつのデータの組で計算する重回帰分析と多数のケースを集計した値で計算するロジスティック回帰分析とはデータの質がちがいます。

さらに、通常はフィールドの調査結果では、ケースごとにサンプル数が異なることが普通です。だから、ロジスティック回帰式の係数計算には通常の重回帰分析ソフト（最小二乗法）では計算せず、最尤法という方法で求められています。
ケース毎のサンプル数を同じにできるマウスなどの実験なら、通常の重回帰分析で対応できます。

 基本知識１：　オッズとロジット変換

オッズがよく使われてきたのはきったはったの勝負の世界。勝ち確率／負け確率　の比がオッズです。勝つ確率は負ける確率の何倍かという数値ですな。オッズが大きいほうが勝ちを呼ぶわけだね。

地域によっては、分子と分母が逆の場合のオッズもあるから、そういう時はオッズが小さい方が勝ちを呼ぶ指標です。

そしてオッズの自然対数値を ロジット と呼びます。これも対数値に名前を付けただけなのだから深く考えないこと。ロジットの Ｌ のギリシャ語は λ（ラムダ）なのかな。ロジットを λ と表現します。
（１）式の計算でオッズをロジット値にすることをロジット変換といいます。

 基本知識２：　ロジットを目的変数にした重回帰分析が ロジスティック回帰

マウスで実験する場合なら、マウスの検査数を同じに出来るからこの計算は簡単です。いろいろな実験条件（= 多変数）を複合させて試験して、マウスの死亡率（または治癒率など）を求めると、通常の重回帰分析で計算できます。

目的変数をロジット値に変換しておいて重回帰分析すれば（２）式が完成します。ロジット値の９５％信頼区間は普通に重回帰分析から計算できるから、確率 p に逆算すれば確率値での信頼区間が計算できるね。

オッズは多くの因子によって、変わってゆくわけだけど、見方を変えるとオッズの対数値は多くの因子の合成変数なのです。
カンがいい人は競馬のいろいろな情報（馬体重、つやなど）をいれてロジスティック解析すれば勝つ確率が計算できるじゃんと思うよね。入れる変数が勝ちの変数じゃないと馬く（笑って♪）いかないから、筆者は責任をとらないぞ。

ロジスティック回帰で変数選択したい人たちがいます。

お医者さんは、下で説明するオッズ比を使いたいから、ロジスティック回帰式の係数の推定精度が正しいことを願うはずです。この場合は、重回帰分析の変数選択法を用いるやり方が推奨されます。

ただし、関係の強い変数を犠牲にしての結果ですから、選ばれてきた変数だけが病気の発生確率に関与する因子であると報告してはいけません。ステップワイズ法などの変数選択法では、係数の推定精度を上げるためにとりあえず関係の深い二つの変数の一つを捨てただけなのですから。

λと X1、X2、･････Xnとで主成分分析をおこない変数選択（手法詳細は第3章へ）をすれば、肝心の変数を排除する重回帰的変数選択よりわかりやすい変数群が得られます。（詳細理由は第5章の多重共線性へ）

 基本知識３：　ロジスティック曲線とLD50

（３）式 λ のところに（２）の重回帰式の右辺が入ると、ロジスティック回帰式というわけです。

回帰が投薬量と死亡率の2変数の場合には LD50 が簡単に求められます。 ５０％死亡確率（ p= 0.5 ）のときのλ値が（３）式より求められるから、（２）式に代入して投薬量を求めると LD50 になります。

（３）式を見れば、なぜロジット値を対数値ではなく e を底にした自然対数値にしたかわかりますね。簡単で美しい式にできるからです。

 基本知識４：　オッズ比

ロジスティック回帰式のどれかの変数 Xi についてだけ,その値を a ← b にしたときのロジット値 λ を λa 、λb とします。そのとき対応する p 値を p_a 、 p_b とします。その他の変数は同一条件（同じ値）だと仮定するとロジスティック回帰式はとても便利な性質があります。 （２）式を使ってロジット値の差をとると、その他の変数は同じ条件だから消えて（４）式ができます。

  （４）式はオッズ a とオッズ b の比の自然対数値を示していいます。オッズ a とオッズ b の比はオッズ比と呼ばれます。

ロジスティック回帰式はすでに計算されていて、係数 βi 値は既知だし、 a 、 b 値は私達が与える任意の既知の値だからオッズ比は（４）式を変形した（５）式で計算できます。

次に説明するオッズ比の利用の説明の都合で、（５）式の中央の式は変形されているけれど、これはよく見るとオッズ比になっているだろ？
オッズも確率の比だけど、オッズ比は「二つのオッズの比」だから混乱しないように。オッズは比で、オッズ比は「比の比」だからね。

  もし p_a 、 p_b が、非常に小さい確率の時には中央の式の右の項は近似的に１になります。だからこの場合オッズ比は近似的に p_a と p_b の比になります。

たとえば通常の健康な人の Ｘi の測定値が b のとき p_b が健康な人の病気発生確率ですが、測定値が a になった場合には病気発生率は何倍になるかという数値としてオッズ比が使えます。ある処置をとって相対危険度（RR）はいくらぐらいかという判断に使うようです。多変量のロジスティック回帰式を求めておくと、一つの式で色々な変数に対して、オッズ比はすぐに計算できて便利ですね。