Microsoft PowerPoint - JMPの非線形回帰handout.ppt

1

JMPの非線形回帰

の使い方と適用事例

－ JMPer‘s Meeting 2007-2 － 2007年7月19日（木） 日本たばこ産業株式会社 杉山 公仁 前東京理科大学教授 芳賀敏郎 2

本日の発表内容

„ 背景 „ 線形と非線形の違い „ 非線形モデルの解き方 … 基礎編 „ 線形化 „ 最小二乗法について „ JMPの非線形回帰 線形モデル … 本編 „ JMPの非線形回帰 非線形モデル（Emaxモデルを例として） „ 効力比 „ ロジスティック回帰（最尤法） „ JMPのユーザーサポート … モデルライブラリー … サンプルデータ

3

背景（その１）

„ 第3回 医薬品開発のための統計解析講座 2007年10月24日（水）開講 エクセル、JMP（ジャンプ）による基礎から応用統計解析実務者コース （創薬セミナー）

創薬アカデミー・SAS Institute Japan株式会社 JMP事業部 共催

„ セミナーの目的 製薬会社で創薬に携わる研究者、実務者ならびに臨床現場でデータを扱う医療関係者は しばしば統計解析で悩まされています。 社内の統計専門家は恒常的に忙しく、簡単な実験データ解析には相談にのってもらい難い のが実態です。「自分の部署である程度までの統計解析を行いたいがSASやJMPを扱える ほどの要員はいない。エクセルならば部員全員が一応は扱える。エクセルで統計解析はで きないものか？」と悩んでいる研究・実務統括者は多いと思われます。 本講座はそのような悩みを解決すべく、エクセルでゴールシーク、ソルバーを駆使してｔ検定 から多重比較、最尤法、ロジスティック回帰までを修得できるような12ヶ月(月1回）の研修 コースを企画いたしました。3部構成で、基礎、実験計画及び拡張編からなっています。 また、グローバルなSAS社統計ソフトJMPの使い方もマスターできるようなコースです。 本研修終了者は担当部署の業務に関連する統計解析ができるようになり、社内外の統計 家と専門的な会話が可能になるよう、人材育成を目指します。 4

背景（その２）

„

創薬セミナーではナビゲータ方式を採用．

„

「ナビゲータ方式」とは，芳賀が作成したテキ

ストを統計を勉強しているメンバー(セミナー

の一期生を含む）が理解し，理解したことを分

かり易く受講生に伝える方法であり，ナビ

ゲータがつまづくところは受講生もつまづきや

すい点であるため，芳賀が補足説明をすると

いう講義方法のことである．もちろん，ナビ

ゲートするには事前のリハーサルを繰り返し

実施し，十分な準備を行っている．

5

背景（その３）

„ 創薬セミナーの第3部の拡張編では以下の項目を講義している． … 非線形モデル ・ 指数曲線 ・ ロジスティック曲線 ・ 非線形最小２乗法による推定値の標準誤差 ・ 重み付き最小２乗法 … 計数値 ・ 2変量の関係 ・ リスク・オッズ・ロジット ・ オッズとオッズ比 ・ ロジスティック回帰分析・ D50検定 ・ 多項ロジスティック回帰分析 … 生存時間解析 ・ Kaplan-Meier曲線 ・ ログランク検定 ・ 一般化ウィルコクソン検定 ←浜田知久馬先生 6

本発表では，JMPセミナーの一部をご紹介します．

創薬セミナーの第3部からJMPの「非線形回帰」に 関連する部分をJMPのセミナーとして9月に実施す る予定です．

7

線形と非線形

„

次の式1～4は線形？ それとも，非線形？

0 20 40 60 80 100 120 140 0 2 4 6 8 10 12 x y 1 2 3 4 1 0 2 1 2 2 1 0 1 0 : 4 ) ( 50 : 3 : 2 : 1 b x b y b x b y x b x b b y x b b y = − + = + + = + = 式 式 式 式 パラメータの値 b0 b1 b2 式1： 5 10 － 式2： 1 10 -0.2 式3： - 10 5 式4：10 1.1 － 8

解答

„ 式1は直線で，線形です． „ 式2は曲線ですが，線形です． „ 式3は直線ですが，非線形です． „ 式4は曲線で，非線形です． „ 線形と非線形の違いは何でしょうか？ „ グラフの形(直線 or 曲線）ではないことは明らかです． 1 0 2 1 2 2 1 0 1 0 : 4 ) ( 50 : 3 : 2 : 1 b x b y b x b y x b x b b y x b b y = − + = + + = + = 式 式 式 式 0 20 40 60 80 100 120 140 0 2 4 6 8 10 12 x y 1 2 3 4

9

パラメータに注目

„ 式1は でパラメータ b₀と b₁に関して一次 式で表せます． „ 式2は でパラメータ b₀，b₁と b₂に関して 一次式で表せます． „ では，式3と式4はどうでしょうか． „ パラメータに関して一次式では表せません． x b b y = 0 + 1 2 2 1 0 bx b x b y= + + 1 0 b x b y=

つまり，x に関してではなく，パラメータに関して

一次式で表せるモデルが線形モデル，一次式で

は表せないモデルが非線形モデルとなります．

) ( 50 b1 x b2 y = + − 10

問題（その１）：

について，

「任意の反応を示すx の推定値は？」

„

式1は線形モデルでした．

„

では，ある反応(y = 10)となる x を求めるに

はどうしたらよいだろうか？

1.

求めたパラメータを使って，

y = 10

を代入．

2.

式を変形したらどうだろうか？

つまり，求めたい値をモデルにする．

3.

JMPでは逆推定の問題

x

b

b

y

=

₀

+

₁ 0 2 4 6 8 10 12 0 5 10

11

「任意の反応を示す x の推定値は？」

1. 求めたパラメータを使って， y = 10 を代入．

⇒

として解が得られる．

⇒ x

₁₀

の信頼区間は別に求める必要がある．

1 0 10 10 b b x = −

式の変形

)

(

10

_{1 10} 1 0

b

x

b

x

x

b

y

=

+

=

+

−

この式は式3と同じです

2. パラメータとして y = 10 となる x をモデル式に

入れる．

求めるパラメータはどれ？

12

「任意の反応を示す x の推定値は？」

2. パラメータとして x

₁₀

をモデル式に入れる．

„

式変形でパラメータ b

₁

と x

₁₀

については，一次

式で表せなくなり，非線形モデルになった．

⇒JMPでは「非線形回帰」で解析可能 ⇒非線形モデルとして解析すると．原則として，パラメー タの点推定値と信頼区間が同時に求められる． ここで，

x

₁₀ は y = 10 となる x であり，パラメータである．

)

(

10

_{1 10} 1 0

b

x

b

x

x

b

y

=

+

=

+

−

13

「任意の反応を示す x の推定値は？」

3. JMPでは逆推定の問題

⇒「モデルのあてはめ」の逆推定で解析可能．

⇒信頼区間も求められる．

14

逆推定とその区間推定

„ JMP「モデルの当てはめ」と逆推定 „ JMP「非線形回帰」 切片 x 項 2.8571429 1.2261905 推定値 1.300973 0.257631 標準誤差 2.20 4.76 t値 0.0705 0.0031 * p値(Prob>|t|) パラメータ推定値 10.000000 y 5.82524272 予測値 x 4.64490828 下側限界 7.95782316 上側限界 0.9500 1-Alpha 期待応答についての信頼区間 逆推定 16.726190476 SSE 6 DFE 2.7876984 MSE 1.6696402 RMSE b x10 パラメータ 1.2261904762 5.8252427184 推定値 0.25763108 0.55613943 近似標準誤差 0.59578994 4.64491383 下側信頼限界 1.85659101 7.95781372 上側信頼限界 次の値で解く: 解析 NR 解 x y y-hat 1 3 4.08 2 7 5.31 3 5 6.54 4 8 7.76 5 10 8.99 6 12 10.21 7 9 11.44 8 13 12.67 b 1.226 x10 5.825 S 16.726

15

0

2

4

6

8

10

12

14

16

0

5

10

0

2

4

6

8

10

12

14

16

0

5

10

パラメータ x

₁₀

の信頼区間

10

x

y の信頼区間 は点推定値の 上下で等間隔 x の信頼区間 は点推定値の 前後で等間隔 ではない． y = 10 x10 の95%下側 信頼限界 x10 の95%上側 信頼限界 推定値 y の信頼区間^ 非線形モデルの パラメータの信頼 区間は等間隔で はない． 16

信頼限界の求め方１

„ Excel Solver で非線形モデルを当てはめる（前 シート右上の結果が得られる）． „ x10 を推定値にδを加えて固定する．変化させる セルにｂだけを指定してソルバーで解く． „ δとSeの関係をグラフにする． „ 放物線（２次曲線）とはならない． „ ３次式を当てはめる．

17 y = -1.9237x3 + 8.3324x2 - 1.1155x + 16.922 10 15 20 25 30 35 40 45 -2.0 -1.5 -1.0 -0.5 0.0δ 0.5 1.0 1.5 2.0 2.5 3.0 Se

信頼限界の求め方２

18

信頼限界の求め方３

„ 残差平方和Se=16.726 を自由度 6 で割って，平 均平方 Ve=2.788を求める． „ Se+Ve=19.514に横線が引かれている． δ=±0.5 の近くで交わっている． „ または，3次式の2次の項8.3324 から，SeがVe だけ増加するδを√（2.788/8.3324）＝0.578 とし て計算される．これが近似標準誤差 0.556に対 応する． „ Se+F(1,6;0.05)Ve=33.417 の位置にも横線が引 かれている．交点の横座標値を推定値に加える と信頼限界が得られる．

19

「モデルのあてはめ」の逆推定と

「非線形回帰」の使い分け

„

「モデルのあてはめ」の逆推定

…「モデルのあてはめ」は線形モデルのみ …線形モデルにおいて，任意の反応 y を示す x の 値が知りたい場合に有効． „

「非線形回帰」

…線形モデルだけでなく，非線形モデルが解ける． 20

小まとめ

„

線形と非線形の違い

…パラメータに関して一次式でないのが非線形モデル „

線形モデルは｢ニ変量の関係」or「モデルのあて

はめ」で解く．

…逆推定は「モデルのあてはめ」で実行可能 „

非線形モデルは「非線形回帰」で解く．

…非線形モデルのメリットは求めたい値をパラメータと してモデル式にして，直接求められる．

21

非線形モデルの解き方 基礎編

„

線形化による解き方

„

最小二乗法による解き方

„

線形モデルを例に手順の説明

„

JMP「非線形回帰」の使い方

…線形モデル 22

非線形モデルの解き方（線形化）

„

式4

はパワーモデルと呼ばれ，薬

物動態の線形性を評価する非線形モデル

„

薬物動態パラメータAUCとCmaxは対数正規

分布に従うことが知られている．

„

誤差の構造を考慮して，変数変換を行う．

0 20 40 60 y 0 20 40 60 x 0 20 40 60 y 0 20 40 60 x b0=1.0 b0=1.0 b0=0.65 b0=0.65 b0=1.3 b0=1.3 1 0 b

x

b

y

=

0 20 40 60 y 0 20 40 60 x 0 20 40 60 y 0 20 40 60 x b1=0.5 b1=0.5 b1=1.0 b1=1.0 b1=2.0 b1=2.0

23

問題（その２）

„

式4を対数変換してみてください．

) log( ) log( ) log( ) log( 1 0 0 0 1 1 x b b x b y x b y b b + = = =

ここで， log(y)=Y,log(b₀) = B₀,log(x) = X とすると，

X b B Y = ₀ + ₁ となり，先ほどの式1と同じでパラメータ B₀，b₁に関して 線形となります． 0 20 40 60 80 100 0 2 4 6 8 10 x y b₀= 1, b₁= 2 24

問題（その２）

„

式4を対数変換してみてください．

) log( ) log( ) log( ) log( 1 0 0 0 1 1 x b b x b y x b y b b + = = =

ここで， log(y)=Y,log(b₀) = B₀,log(x) = X とすると，

X b B Y = ₀ + ₁ となり，先ほどの式1と同じでパラメータ B₀，b₁に関して 線形となります． パワーモデルでは，このような解析が実施可能． しかし，全ての非線形モデルで有効な解析法ではない． 一般的に，線形化は誤差構造が把握しにくい． 0 20 40 60 80 100 0 2 4 6 8 10 x y 0.01 0.1 1 10 100 0.1 1 10 x y b₀= 1, b₁= 2 log(b₀) = 0, b₁= 2

25

非線形モデルの有用性

„ 求めたいパラメータをモデル式で記述できる． „ 観察された現象やその分野における知見に基づき， 自由にモデル化ができる．

JMP「非線形回帰」の有用性

„ 非線形モデルを線形化せずに，パラメータが直接 求められる． „ 最小二乗法で解く． 非線形最小二乗法を説明する前に，最小二乗法について復習しましょう． 26

回帰分析

„ 次のシートのデータに直線 y = b₀+b₁x を当てはめる． „ 下のスクロールバーを左右に移動すると，b₀ ，b₁が 0.01 きざみに変化し，それに伴い，y-hat が変化する． „ 右のグラフで，観測点 y と直線 y-hat = b₀+b₁x ができ るだけ接近するようにしたい． „ S =Σ(y － y-hat)2_{を不一致度の指標とし，S を最小と} するb₀ ， b₁ を求める． „ これが最小二乗法である．

27

最小二乗法の考え方

x y y-hat 1 3 2 2 7 3.95 3 5 5.9 4 8 7.85 5 10 9.8 6 12 11.75 7 9 13.7 8 13 15.65 b0 0.05 5 b1 1.95 195 S 40.35 0 5 10 15 20 0 5 10 y y-hat 与えた初期値から試行錯誤しながら，最適解を求める． EXCELへ 28

ソルバー

„

試行錯誤を自動的に実行してくれるのが

Excel のソルバー

である．

„

トップメニュー「ツール」から「ソルバー」を選択

する．次シートが現われる．

„

目的セルに S のセルを，目標値に最小値を，

変化させるセルに b

₀

, b

₁

のセルを指定して，

実行する．

„

y = 2.86 + 1.23x, S = 16.73 が得られる．

29

ソルバーの設定画面

EXCELへ 30

線形最小二乗法

「二変量の関係」による解析

„ 回帰分析のモデル y = b₀ + b₁x を例とした，これまでの 解析手順は，非線形最小二乗法の考え方を示すため のものである． „ 線形モデルの場合，実際にはJMP「二変量の関係」(あ るいは「モデルのあてはめ」)で解析する． „ 解析的に，正規方程式を解いて，解を求める． „ 直線を当てはめると次の出力が得られる． 切片 x 項 2.8571429 1.2261905 推定値 1.300973 0.257631 標準誤差 2.20 4.76 t値 0.0705 0.0031 * p値(Prob>|t|) -0.326224 0.5957899 下側95% 6.0405095 1.856591 上側95% パラメータ推定値

31

JMPの「非線形回帰」による解析

手順

1.

JMPデータシートにモデル式を含む列を追加

2.

初期値の設定

3.

数値計算（試行錯誤）の実行

4.

信頼限界の算出

JMPの「非線形回帰」は

Excel のソルバーで最適解

を求めるのと同じ手順

で解析を進める．

信頼限界の算出

JMP へ 32

JMPの「非線形回帰」による解析

手順１

„

JMPデータシートにモデル式（計算式）を含む列

を追加する．

…y-hat の列を作り，「計算式」を入力する． …「計算式」の画面で，「パラメータ」を選択し，「パラ メータの新規作成」でb₀,b₁ を作成し，適当な初期値 （たとえば， b₀ = 2, b₁= 1）を指定する． …計算式 b₀ + b₁ x を入力する．

33

JMPの「非線形回帰」 手順１

《 y-hat の列作成と計算式入力》

b0 + b1 * x 34 „

「分析」→「モデル化」→「非線形回帰」を選択

する．

…パラメータ指定画面で，y を「Y, 応答変数」に，y-hat

を「X, 予測式列」に指定する． …OK をクリックする． …必要があればパラメータの初期値をプロットを見な がら設定する． …「実行」をクリックする． …「信頼限界」をクリックする．

JMPの「非線形回帰」による解析

手順２

35

JMPの「非線形回帰」 手順２

《パラメータ指定画面》

36

JMPの「非線形回帰」 手順３

《初期値の変更》

37

JMPの「非線形回帰」による解析結果

16.726190476 SSE 6 DFE 2.7876984 MSE 1.6696402 RMSE b0 b1 パラメータ 2.8571428571 1.2261904762 推定値 1.30097317 0.25763108 近似標準誤差 -0.3262238 0.59578994 下側信頼限界 6.04050952 1.85659101 上側信頼限界 次の値で解く: 解析 NR 解 切片 x 項 2.8571429 1.2261905 推定値 1.300973 0.257631 標準誤差 2.20 4.76 t値 0.0705 0.0031 * p値(Prob>|t|) -0.326224 0.5957899 下側95% 6.0405095 1.856591 上側95% パラメータ推定値

「二変量の関係」の結果と一致する．

38

非線形モデルの解き方 本編

„ データへ非線形モデルのあてはめ „ JMP「非線形回帰」の使い方 …非線形モデル

39

非線形モデルによる非線形回帰

計量値の例

„ Emaxモデル γ γ γ EC50 x x Emax x f + ⋅ = ) ( Emax：最大反応量 x：処置量(投与量あるいは濃度) γ：係数(傾き) EC50：50%反応量(濃度)=D50 Emaxは実験の 最大投与量の 反応量ではあり ません 推定されるパラメータ „ 誤差が正規分布であると仮定した非線形モデルについてご 紹介します 40

E

_max

モデルの説明

E

_max

モデルからの式変形

γ γ γ 50 max ) ( EC x x E x f + ⋅ = _{分子、分母を x}γ_で割る γ γ x EC E x f 50 max 1 ) ( + = 指数で表わすと )) ln( ) exp(ln( 1 ln exp 1 ) ( 50 max 50 max γ γ γ γ _{EC x} E x EC E x f − + = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + =

(

)

(

ln( ) ln( )

)

exp 1 ₅₀ max x EC E − + =

γ

41

In vitro 実験例

《実験装置》

マグヌス（Magnus）装置 1μM ヒスタミン 42

In vitro 実験例

《データ》

1μM ヒスタミン 3.16μM ヒスタミン 10μM ヒスタミン 66 mm 113 mm 158 mm

43

計量値データ

使用するデータ例

„ ヒスタミンによる摘出腸管の収縮反応 „ 0.01～316μMの10用量 „ n=1 ヒスタミン濃度 平滑筋 （μM） 収縮量 （mm） 0.01 1 0.0316 3 0.1 5 0.316 23 1 66 3.16 113 10 158 31.6 171 100 171 316 165 γ γ γ EC50 x x Emax x f + ⋅ = ) ( 0 20 40 60 80 100 120 140 160 180 0 100 200 300 dose 収縮量 （ m m ） x を実数でプロット データの出典：http://www.yukms.com/biostat/takahasi/rec/004.htm 44

計量値データ

使用するデータ例

„ ヒスタミンによる摘出腸管の収縮反応 „ 0.01～316μMの10用量 „ n=1 ヒスタミン濃度 平滑筋 （μM） 収縮量 （mm） 0.01 1 0.0316 3 0.1 5 0.316 23 1 66 3.16 113 10 158 31.6 171 100 171 316 165 0 20 40 60 80 100 120 140 160 180 0.01 0.1 1 10 100 1000 dose 収縮 量（ mm ）

(

)

(

ln( ) ln( )

)

exp 1 ) ( 50 max x EC E x f − + = γ x を対数でプロット データの出典：http://www.yukms.com/biostat/takahasi/rec/004.htm

45

計量値データ

使用するデータ例

„ ヒスタミンによる摘出腸管の収縮反応 „ 0.01～316μMの10用量 „ n=1 ヒスタミン濃度 平滑筋 （μM） 収縮量 （mm） 0.01 1 0.0316 3 0.1 5 0.316 23 1 66 3.16 113 10 158 31.6 171 100 171 316 165 0 20 40 60 80 100 120 140 160 180 0.01 0.1 1 10 100 1000 dose 収縮 量（ mm ）

(

)

(

ln( ) ln( )

)

exp 1 ) ( 50 max x EC E x f − + = γ x を対数でプロット データの出典：http://www.yukms.com/biostat/takahasi/rec/004.htm 0 20 40 60 80 100 120 140 160 180 0.01 0.1 1 10 100 1000 dose 収縮 量（ m m ） S 字曲線（シグモイド曲線） 46

シグモイド曲線

0.0 0.1 0.2 0.3 0.4 0.5 -4 -3 -2 -1 0 1 2 3 4 正規分布 ﾛｼﾞｽﾃｨｯｸ分布 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 -4 -3 -2 -1 0 1 2 3 4 プロビット曲線 ﾛｼﾞｽﾃｨｯｸ曲線 正規分布 ⇒プロビット曲線 ロジスティック分布 ⇒ロジスティック曲線 „

正規分布あるいはロジスティック分布の

分布関数はシグモイド曲線となる．

プロビット曲線とロジスティック曲線は

ほとんど違いはありません

47

シグモイド曲線の比較

„

正規分布の分布関数

（プロビット曲線） „

ロジスティック分布の分布関数

（ロジスティック曲線）

∫

−∞ ⎟⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − − ⋅ = x dt t x F 2 1 exp 2 1 ) ( 2 σ μ σ π ( ) (₍ )_{) ( ( ))} x x x x F ⋅ + − + = ⋅ + + ⋅ + = 1 0 1 0 1 0 exp 1 1 exp 1 exp β β β β β β „ Emaxモデル γ γ γ D50 x x Emax x f y + ⋅ = = ( )

₍

₍

₎

₎

Emax x D50 − ⋅ + = ) ln( ) ln( exp 1 1 γ

Emaxモデルによるシグモイド曲線はロジスティック

曲線と本質的には一致している．

48 値指定 値指定範囲 EC50 = 1 (0.01～100) γ = 2 (-5～5) Emax = 8 (0～10) 0 2 4 6 8 10 0.01 0.1 1 10 100 x ( ) ( ln( ) ln()) exp 1 ) ( 50 max x EC E x f − + = γ γ γ γ 50 max ) ( EC x x E x f + ⋅ =

EXCELで確かめるパラメーター

の意味

EXCELへ

49

JMPの「非線形回帰」

《 y-hat の列作成と計算式入力》

b0 + b1 * x „ Emaxモデル 自然対数 50

JMPによる解析結果

-50 0 50 100 150 200 y 0 1 100 x Emax EC50 gamma パラメータ 171.57994341 1.5873589117 1.1677651165 推定値 85 0.795 0.6 最小値 255 2.385 1.8 最大値 プロット 121.1020042 SSE 7 DFE 17.300286 MSE 4.1593613 RMSE Emax EC50 gamma パラメータ 171.57994341 1.5873589117 1.1677651165 推定値 2.68144664 0.11320162 0.08213898 近似標準誤差 165.519251 1.34387134 1.00246152 下側信頼限界 177.952336 1.87906987 1.37033761 上側信頼限界 次の値で解く: 解析 NR 解

51

用量反応関係式の類似性

γ γ γ

D50

x

x

E

y

max

+

⋅

=

n D n F n F max B K ] L [ ] L [ B ] L [ + = h h h K S S V V 5 . 0 max ] [ ] [ + = Emax理論 受容体結合実験 酵素反応速度論 Michaelis-Menten M K S S V V max + ⋅ = ] [ ] [ h = 1 の場合 52

小まとめ

„

JMPの「非線形回帰」の手順を線形モデルで，

説明した．

„

「非線形回帰」では，最小二乗法によって，パ

ラメータの最適値を初期値から試行錯誤的に

数値計算して求める．

„

Emaxモデルを例にして，パラメータの意味，

「非線形回帰」での解析手順を示した．

„

Emaxモデルが酵素反応などの解析にも利用

できることを示した．

53

その他のモデル

„ 非線形回帰は極めて広範囲のモデルに当てはめ ることができる． „ それらのモデルの一部は，JMPのマニュアルに取 り上げられている． „ そのデータは JMP のサンプルデータに登録され ている． „ JMPを立ち上げ，「ヘルプ」から「サンプルデータ ディレクトリ」を選択すると，一覧表が表示される． „ それらの一部が 配布される「補助資料」の中で解 説されている． 54

ロバスト回帰

„

次のデータはアン

スコムが提供した外

れ値を含むデータで

ある．

„

最小二乗法による

回帰直線は点線の

ように，外れ値に

引っ張られる．

x y y-hat e 4 5.39 5.39 0.000 5 5.73 5.73 -0.005 6 6.08 6.08 0.000 7 6.42 6.42 -0.005 8 6.77 6.77 0.000 9 7.11 7.11 -0.005 10 7.46 7.46 0.000 11 7.81 7.81 0.005 12 8.15 8.15 0.000 13 12.74 8.50 4.245 14 8.84 8.84 0.000 4 6 8 10 12 14 4 6 8 10 12 14 y 線形 (y)

55

ロバスト推定

„ きたないデータ（外れ値を含むデータ）では平均値 よりも中央値の方が代表値として安定した値が得ら れることが知られている． „ 平均値は残差の2乗の和を最小にする最小二乗 推定値である，中央値は残差の絶対値の和を最小 とする推定値である． „ アンスコムのデータに，残差の絶対値の和を最小 とする回帰式を当てはめることを考える． „ また，中間を取って，残差の絶対値の 1.5 乗の和 を最小にすることも考えられる． 56

「非線形回帰」の別の使い方

„ 「非線形回帰」の標準的な使い方は， y-hat を 「X， 予測式列」に，ｙ を「Y, 応答変数」に指定すると， JMPは両者の差(残差）の2乗和が最小となるように， y-hat の中のパラメータを決める． „ もう一つの方法は，ｙ を指定せず，残差の2乗に相 当する列を準備し，それを「損失」に指定する．JMP は「損失」の合計が最小になるようにパラメータを決 める．

57 „

次のようなデータ表を準備する．

„

y-hat，loss には次の計算式を入力する．

58 „ y-hat 列 を 「X,予測式列」に， loss 列を「損失」に 指定する． „ メニューの一番下の「数値微分のみ」の左の☐ をクリックして5 を表示する．（次シート） „ 「OK」 をクリックする． „ 「非線形回帰のあてはめ」で「実行」をクリックす ると次次シートの解が得られる．

59 60 2.5 5 7.5 10 12.5 15 Y 2.5 5.0 7.5 10.0 12.5 15.0 x Y y y-hat 重ね合わせプロット 8.4461669554 損失 9 DFE 0.938463 平均損失 0.968743 平均損失平方根 b0 b1 パラメータ 3.7328054482 0.3882806004 推定値 0.82459453 0.10083384 近似標準誤差 . . 下側信頼限界 . . 上側信頼限界 次の値で解く: 数値 SR1 解

61

ちょっと，休憩

62

効力比

„

効力比の定義

„

Emaxモデルの拡張

„

データ例による効力比の算出

63

効力比の推定

„

効力比の定義：

A剤の投与量 x と薬効 y の間に y = f (x) の関係が あるものとする．ここで関数 f (x) はどのようなもの であっても構わない． B剤の投与量と薬効の関係が y = f (cx) で表わさ れるとき， c が効力比である． 0 2 4 6 8 10 12 14 16 0 2 4 6 x A B 0 2 4 6 8 10 12 14 16 -1 0 1 2 log(x) A B 1 c ln( c ) 実数 対数 ln(x) 64

実験データ

薬剤A の投与量を，0.192，0.48，1.2，3.0とほぼ

等比級数的に変化させて，4匹の動物に投与し

て，薬効を調べる．

薬剤B の効力は薬剤A の半分程度（効力比 c =

0.5）と想像されるので，薬剤A の2倍，すなわち，

0.384，0.96，2.4，6.0 と変化させて投与した．

効力比を求める実験データ 薬剤 x X=ln(x) y-bar A 0.192 0.93 1.17 1.03 1.19 -1.65 1.08 0.480 1.08 1.61 1.60 1.64 -0.73 1.48 1.200 1.61 1.76 1.75 1.70 0.18 1.71 3.000 1.73 1.63 1.66 1.84 1.10 1.72 B 0.384 0.71 1.00 0.91 0.96 -0.96 0.90 0.960 1.26 1.32 1.18 1.31 -0.04 1.27 2.400 1.53 1.68 1.82 1.65 0.88 1.67 y

65 対数目盛 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 -2.00 -1.00 0.00 1.00 2.00 A B

グラフ化

66 対数目盛 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 -2.00 -1.00 0.00 1.00 2.00 A B

グラフ化

用量反応関係を表現するモデルを当てはめるの

がよいだろう．

ただし，反応は0から始まるのではなさそう．

67

Emaxモデルの拡張

„

Emaxモデルの y の範囲は 0～Emaxの範囲

であった．

„

任意の y( = y

_min

)から最大値 y

_max

までの範囲

に拡張する必要がある．

68

拡張したモデル

拡張したモデルの X₅₀は？

Emaxモデル：y

_min

= 0， y

_max

= Emax

2 2 min max min max min y y y y y y = + − = +

となる X である．

ここで，X=ln(x) あるいは， X₅₀=ln(x₅₀)とする． )) ( exp( 1 )) ( exp( 1 ₅₀ min max min 50 min max min X X B y y y X X B y y y y − − + − + = − + − + = Emaxモデルの X₅₀は となる x であった．

2

max

E

69

ロジスティック曲線の当てはめ

) ( exp( 1 ) ( exp( 1 ₅₀ min max min 50 min max min X X B y y y X X B y y y y − − + − + = − + − + = 薬剤Aのモデル式 ここで，X=ln(x) あるいは， X₅₀=ln(x₅₀)とする． 薬剤B の効力比が c のとき，薬剤B の x₅₀は薬剤A の x₅₀ の 1/c となる（薬剤B が2倍の効力比と持つとき，薬剤A の 半分の投与で同じ効果が得られる）．これから， ))) ln( ( exp( 1 ))) / ln( ( exp( 1 50 min max min 50 min max min c X X B y y y c x X B y y y y + − − + − + = − − + − + = EXCELへ 70

JMPによる解析

0.4294967837 SSE 27 DFE 0.0159073 MSE 0.1261241 RMSE ymin ymax B X50 ρ パラメータ 0.7993648439 1.7319219573 2.2933720753 -1.234178388 -1.173473463 推定値 0.16773849 0.04853574 0.83905082 0.2142144 0.14962186 近似標準誤差 . 1.64834873 . . -1.4841826 下側信頼限界 1.02354324 . . -0.8248258 -0.8562815 上側信頼限界 次の値で解く: 解析 NR 解 ))) ln( ( exp( 1 ₅₀ min max min c X X B y y y y + − − + − + = JMPの解析では ln(c)=ρとして いる． ρ= ln(c)= -1.173 より，c = exp(ln(c))=0.309 である．

71

JMPによる解析

0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 A B A B y₅₀ = (1.73+0.80)/2 = 1.266 AのX₅₀= -1.234 BのX₅₀= -1.173-(-1.234) = -0.061 ρ= ln(c) = -1.173: c=0.31, cの95%信頼区間(0.23, 0.42) 72

小まとめ

„

Emaxモデルを拡張し，y の範囲を y

_min

から

y

_max

までにした．

„

効力比を定義した．

„

2本の曲線を「非線形回帰」で当てはめた．

…2剤の選択は If 関数で識別

73

ロジスティック回帰

„ p の範囲は0から1で，誤差は二項分布． „ 最尤法で解く． „ JMPの「ニ変量の関係」あるいは「モデルのあて はめ」では解析的に解く． „ 「非線形回帰」では損失関数に尤度関数を設定し， 数値計算で解く． 74

ロジスティック回帰のモデルの意味

„

pは0～1の範囲の値しかとらない．

„

オッズ

は0～∞の値しかとらない．

z p p p _⎟⎟= ⎠ ⎞ ⎜⎜ ⎝ ⎛ − = 1 ln ) ( logit „

対数オッズは-∞～∞の値をとることができる．

対数オッズはロジットlogit と呼ばれる． ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − p p 1 X z=β₀+β₁⋅ と定義したとき， の関数を考える．

75

pとｚとの関係

-∞～∞ 0 ～ 1 シグモイド曲線(0 ~ 1)を直線(-∞ ~ ∞)に変換して いることを意味する．（線形化） 76

ロジスティック回帰のモデルの意味

„

pは0～1の範囲の値しかとらない．

( ) (β β ) ( (β β )) ε β β ₊ ⋅ + − + = ⋅ + + ⋅ + = X X X p 1 0 1 0 1 0 exp 1 1 exp 1 exp X p p _{= + ⋅} ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − 0 1 1 ln β β このモデルから求められる p は0 ～ 1の範囲に収まる． この式を p について解くと， „

対数オッズは-∞～∞の値をとることができる．

したがって， 誤差 ε は 二項分布 に従う． ロジスティック曲線(0 ~ 1)となる．

77

モデルのバリエーション

„

ロジスティック回帰

…z = A+ Bx の場合： JMP｢ニ変量の関係」 or「 モデルのあてはめ」 …z = a +bx + cx2 _の場合： JMP｢モデルの当てはめ」 …z = B(x -x₅₀) の場合： JMP｢非線形回帰」

)

exp(

1

1

z

p

−

+

=

どのモデルでも誤差は二項分布に従う． よって，ロジスティック回帰は最尤法で解析する． 78

二項分布の確率と尤度

n 10 f .05 .10 .15 .20 .25 .30 .35 .40 .45 .50 .55 .60 .65 .70 .75 .80 .85 .90 .95 0 .60 .35 .20 .11 .06 .03 .01 .01 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 1 .32 .39 .35 .27 .19 .12 .07 .04 .02 .01 .00 .00 .00 .00 .00 .00 .00 .00 .00 2 .07 .19 .28 .30 .28 .23 .18 .12 .08 .04 .02 .01 .00 .00 .00 .00 .00 .00 .00 3 .01 .06 .13 .20 .25 .27 .25 .21 .17 .12 .07 .04 .02 .01 .00 .00 .00 .00 .00 4 .00 .01 .04 .09 .15 .20 .24 .25 .24 .21 .16 .11 .07 .04 .02 .01 .00 .00 .00 5 .00 .00 .01 .03 .06 .10 .15 .20 .23 .25 .23 .20 .15 .10 .06 .03 .01 .00 .00 6 .00 .00 .00 .01 .02 .04 .07 .11 .16 .21 .24 .25 .24 .20 .15 .09 .04 .01 .00 7 .00 .00 .00 .00 .00 .01 .02 .04 .07 .12 .17 .21 .25 .27 .25 .20 .13 .06 .01 8 .00 .00 .00 .00 .00 .00 .00 .01 .02 .04 .08 .12 .18 .23 .28 .30 .28 .19 .07 9 .00 .00 .00 .00 .00 .00 .00 .00 .00 .01 .02 .04 .07 .12 .19 .27 .35 .39 .32 10 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 .01 .01 .03 .06 .11 .20 .35 .60 π f = 3 0.0 0.1 0.2 0.3 0.0 0.2 0.4π0.6 0.8 1.0 尤度 π= 0.2 0.0 0.1 0.2 0.3 0 1 2 3 4f5 6 7 8 9 10 確率 f = 3 が得られたとき，尤度が最大となる πは0.3となる．これが最尤法．

79

最尤法

) ( ) ˆ 1 ( ˆ ) ˆ ( i i i i i f n i f i f n i C p p p L = − − 尤度は， 尤度は，Excelでは JMPでは で定義される．

(

f n p

)

BINOMDIST , , ˆ ) ( , ) exp( 1 1 ˆ z x x₅₀ z p _{i i} i i= _{+ −} =β − のとき，

(

p n f

)

y Probabilit Binomial ˆ, , „ 全体としての尤度は個々の尤度の積となり，対数をとること によって個々の尤度の和として求めることができる． „ （対数）尤度が最大となるように， を定めるパラメータ β と x50 を推定するのが最尤法である． „ -2対数尤度はχ2_{値と関連する．} i pˆ 80

の例

(

)

(

B x x

)

p − + = 50 exp 1 1 ˆ

81

列の作成

82

83

JMP解析結果

切片 x 項 -26.211489 5.00433617 推定値 6.350431 1.207517 標準誤差 17.04 17.18 カイ2乗 <.0001* <.0001* p値(Prob>ChiSq) 推定値は次の対数オッズに対するものです： 0/1 パラメータ推定値 0.50000000 確率 5.23775541 予測値 x 5.05542585 下側限界 5.41318710 上側限界 0.9500 1-Alpha 逆推定 11.077731404 損失 4 DFE 2.7694329 平均損失 1.6641613 平均損失平方根 b x50 パラメータ 5.0043403751 5.237755368 推定値 1.20752566 0.07823367 近似標準誤差 . . 下側信頼限界 . . 上側信頼限界 次の値で解く: 数値 NR 解 「ニ変量の関係」 「逆推定」 で解析した結果 「非線形回帰」での 解析結果 信頼区間の推定に おいて収束しない． 84

まとめ

85

まとめ１

„

線形と非線形の区別を明確にした．

„

非線形モデルを解く方法として，最小二乗法

を説明した．

„

非線形最小二乗法は数値計算で最適値を求

める方法であった．したがって，初期値の設

定が重要である．

„

JMPの「非線形回帰」では， Emaxモデルを例

にパラメータとその信頼区間を直接求めた．

86

まとめ２

„

Emaxモデルで，酵素反応速度論などの解析

ができることを示した．

„

Emaxモデルを拡張し，さらに，効力比につい

ても紹介した．

„

Emaxモデルとロジスティック回帰のモデル式

の類似性を示し，ロジスティック回帰における

最尤法について説明した．

„

「非線形回帰」で最尤法が解けることを示した．

87

まとめ３

„

非線形モデルの解析結果として得られた信

頼区間は統計ソフトによって異なる可能性が

あるので，使用しているソフトがどのような信

頼区間を算出しているかを確認することは重

要である．

„

非線形モデルが解けるようになったことで，適

切なモデルを作成するには，データの裏にあ

る専門分野の知識が重要となる．

88

セミナーのご案内

„

2007年9月3日(月)及び18日(火)の2日間で非

線形回帰に関するセミナーを予定．

„

Excelのソルバーを利用して，考え方を説明し，

JMPで追試する．

„

非線形回帰の解法として，最小二乗法と，最

尤法について説明する．

„

効力比を求める非線形モデルについて説明

する．

„

非線形モデル構築における考え方を経時

データを例に説明する．

89

JMPのユーザーサポート

非線形回帰

【モデルライブラリー】抜粋

90

モデルライブラリー

91

ロジスティック2p

"Bates, D.M.and Watts, D.G.(1988), Nonlinear Regression Analysis & its Applications.New York, John Wiley and Sons."

1

1+ theta1 * Exp theta2* X

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Y -5 0 5 X 92 ロジスティック3p

"Bates, D.M.and Watts, D.G.(1988), Nonlinear Regression Analysis & its Applications.New York, John Wiley and Sons."

theta1 1 + theta2 * Exp theta3 *X

0.0 0.5 1.0 1.5 2.0 Y -1.0 -0.5 .0 .5 1.0 1.5 2.0 X

93

ロジスティック4p

S. (Sylvie) Huet et al.(1996), Statistical Tools for Nonlinear Regression: A Practical Guide With S-Plus and R Examples, Springer

theta1 + theta2 - theta1 1 + Exp theta3 * X- theta4

0.0 0.5 1.0 1.5 2.0 Y 0 1 2 3 4 5 6 7 8 X ) ( exp( 1 50 min max min X X B y y y y − − + − + = 94 Michaelis Mentenモデル(2P)

Raymond H. Myers(1986), Classical and Modern Regression With Applications, Pws Pub Co theta1 * X theta2 + X 0 25 50 75 100 125 150 Y .0 2.5 5.0 7.5 10.0 12.5 15.0 X

95

モデルE(2P)

"Draper, N.and Smith, H.(1981), Applied Regression Analysis , Second Edition, New York: John Wiley and Sons, Inc."

theta1 * Xtheta2 0 1 2 3 4 5 6 7 8 9 10 Y 0 1 2 3 4 5 6 7 8 9 10 X 1 0 b

x

b

y

=

96 モデルH(Gompertz成長モデル、3P)

"Bates, D.M.and Watts, D.G.(1988), Nonlinear Regression Analysis & its Applications.New York, John Wiley and Sons."

theta1 * Exp - Exp theta2 - theta3 * X

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Y 0 1 2 3 4 5 6 7 8 9 10 X