1
医薬開発におけるモデルのあてはめ
Excel-ソルバー による非線形最小2乗法・
最尤法・
ロバスト回帰
芳賀敏郎
（
元 東京理科大学）
2007. 9. 29
2本日のテーマ
• 非線形回帰分析 は今後広い分野で活用が 期待される統計 手法である． • 本日の主題は「非線形回帰分析 」であるが，応用統計家 に どのように 説明したら良いかを考え，通常の解説とはかなり 変わったアプローチで説明する． • 極めて身近 な平均値とその標準誤差，信頼区間から入 り， 単回帰分析に拡張する． • 次いで，直線回帰における逆推定問題を取り上げて，非線 形パラメータの意味と，解析結果 が前と異なること示す． • このような準備の後 に，非線形回帰分析を取 り上げる． 3平均値とは
• n 個の観測値 ｘi が得られたとき，代表値 a を決め たい． • 代表値は観測値に全体として近くしたい． • ｘiとa の距離の総合値として，差の2乗の和Sを考 える． • a によってSがどのように 変化するかをグラフ化する． • グラフは次シートに示すように，放物線となり，平均 値で最小となる． • このような推定を最小2乗法という． S =P(xiÄ a)2 4 x 1 3 4 7 4 9 8 12 平均 6.0 S 92 n 8 自由度 7 平均平方 13.14 標準偏差 3.63 標準誤差 1.28 t(0.05) 2.36 下限 2.97 上限 9.03 S+V 105.1 S+FV 165.5S = 8a2Ä 96a + 380 = 8(a Ä 6)2+ 92 = n(a Ä x)2+ S; s e =pV =n
5 • 元のデータから，平方和S，平均平方V，平均 値の標準誤差se，区間推定値などを計算する． • 平均値±se，区間推定値に縦線を引き，曲線と の交点を通る水平線を引く． • 交点の縦軸の値はS+V，S+FVである． • 2次式の，2次の係数５は，パラメータの推定値 の標準誤差を求めるときの分母である． • この関係は，線形最小2乗推定値に拡張できる． • 非線形最小2乗推定値では，少し修正される． F = F(1; n Ä 1; 0:05) = t(n Ä 1; 0:05)2 6
回帰式
• 次の散布図に回帰直線を当てはめる． • 平均値と同様，最小2乗法が用いられる． • 次シートの散布図で，ｙ＝a＋ｂｘの a, b を試行 錯誤で変化させて， を最小とする解を求める．
S = n X i=1 (yiÄ byi)2= n X i=1 (yiÄ (a + bxi))27
•
x y y-hat -2 1 0.00 -1 3 1.50 0 4 3.00 1 7 4.50 2 4 6.00 3 9 7.50 4 8 9.00 5 12 10.50 a 3.000 b 1.500 S 20.000 0 2 4 6 8 1 0 1 2 1 4 -2 0 2 4 8ソルバー とは
• 試行錯誤を自動的に実行してくれるのがソル バーである． • 下の結果が得られる．
x y y-hat -2 1 1.25 -1 3 2.61 0 4 3.96 1 7 5.32 2 4 6.68 3 9 8.04 4 8 9.39 5 12 10.75 a 3.964 b 1.357 S 14.643 0 2 4 6 8 1 0 1 2 1 4 - 2 0 2 4 9 • ソルバーを選択すると，「パラメータ設定」画面が 表示される． • 目的セルに S のセルを，目標値を最小値 に，変 化させるセルに a, b のセルを指定して実行する． 10ソルバーの利用
• ソルバーは，目的とするセルの値を最小とする だけでなく，最大，または，ある値に等しくすること ができる． • 変化させるセルは，一塊りのセルに限らず， “,” を使って自由に組合わせることができる． • 変化させるセルの値に制約を加えることができる． • これらを上手に使うとその適用範囲は無限であ る．後にいくつかの例を紹介する． 11通常の方法
•
ソルバーを使わなくても，ExcelのLINEST関数で， a,b だけでなく，それらの標準誤差も簡単に求めら れる． • LINEST関数の出力 ソルバーの解 x const a,b 1.357 3.964 se 0.241 0.660 R^2 0.841 1.562 sd F 31.698 6 fe SR 77.357 14.643 Se a 3.964 b 1.357 S 14.643 12 • b の標準誤差 0.241 の求め方を調べる．•
b の最小2乗推定値の前後でδｂだけ変化させ て残差平方和の変化を見る． • 2次式が完全に当てはまる．2次の係数は60で ある． y = 60x2 - 4E-06x + 14.643 0 10 20 30 40 50 60 -1.0 -0.5 0.0δｂ 0.5 1.0 Se13 • これから，b の標準誤差を計算する． • この方法ではLINESTの結果 0.241 は得られな い． • 単純にbを変化させたときの，あてはまりが 悪さ を見るのは誤りである． se(by) =qVe 60 = q 14:643=6 60 = 0:201 14 • 上の計算は，下の図の赤線のように，a =4 を固 定してbを変化させて，残差平方和を求めている． • b を変えたあと，青の点線のように，残差が最小 となるように a も変化させる必要がある． -4 -2 0 2 4 6 8 10 12 14 16 - 2 0 2 4 15 • 次シートに示すような計算表を準備し，ΣSを最 小とするｂを求める．前シートの点線に相当する． • 横軸にδ，縦軸にS をとってグラフを描き，2次曲 線を当てはめる ． • 2次の係数 42 から，ｂの標準誤差を計算する． となり，LINEST関数の結果が求められた． y = 42x2 - 2E-11x + 14.643 0 10 20 30 40 50 -1.0 -0.5 0.0 δb 0.5 1.0 Se se(by) =qVe 42= q 14:643=6 42 = 0:241 16 x y -2 1 0.55 4.05 3.35 2.65 1.95 1.25 0.55 -0.15 -0.85 -1.55 -1 3 2.11 4.61 4.11 3.61 3.11 2.61 2.11 1.61 1.11 0.61 0 4 3.66 5.16 4.86 4.56 4.26 3.96 3.66 3.36 3.06 2.76 1 7 5.22 5.72 5.62 5.52 5.42 5.32 5.22 5.12 5.02 4.92 2 4 6.78 6.28 6.38 6.48 6.58 6.68 6.78 6.88 6.98 7.08 3 9 8.34 6.84 7.14 7.44 7.74 8.04 8.34 8.64 8.94 9.24 4 8 9.89 7.39 7.89 8.39 8.89 9.39 9.89 10.39 10.89 11.39 5 12 11.45 7.95 8.65 9.35 10.05 10.75 11.45 12.15 12.85 13.55 δb 0.2 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 a 3.664 5.164 4.864 4.564 4.264 3.964 3.664 3.364 3.064 2.764 b 1.557 0.557 0.757 0.957 1.157 1.357 1.557 1.757 1.957 2.157 S 16.32 41.52 29.76 21.36 16.32 14.64 16.32 21.36 29.76 41.52 ΣS y-hat 248.91 17 • bの標準誤差は，前々シートのグラフに Se+ Ve=14.643+2.440=17.083 の横線を引き，曲 線との交点から求めることができる． • 信頼区間は，Se+F Ve= 14.643+5.987*2.441 =29.25 の横線を引き，曲線との交点から求める ことができる． 0 1 0 2 0 3 0 4 0 5 0 -1.0 -0.5 0.0 δb 0.5 1.0 Se 18 • a,bの両方を変化させてSを計算し，Sの等高線 を描くと次のグラフが得られる． • 等高線は楕円になる． • 前のグラフは鉛直線，後のグラフは斜線での切 り口の断面に対応する．
-0.8 -0.6 -0.4 -0.2 0 0 . 2 0 . 4 0 . 60.8-1.6 -1.2 -0.8 -0.4 0 0 . 4 0 . 8 1 . 2 1 . 6 a b 60.0 -70.0 50.0 -60.0 40.0 -50.0 30.0 -40.0 20.0 -30.0 10.0 -20.0 0.0 -10.0
19
逆推定
• y=8 となるｘ（x8で表す）はいくらだろうか？ • 散布図に回帰直線とその信頼区間の幅をつけ たグラフを左下に示す． • y=8 の水平線と回帰直線との交点が解である． • 式で計算するとx8=(8-a)/b=(8-3.96)/1.36=2.974 - 2 0 2 4 6 8 1 0 1 2 1 4 -2 -1 0 1 2 3 4 5 x y y-hat - 2 1 1.25 - 1 3 2.61 0 4 3.96 1 7 5.32 2 4 6.68 3 9 8.04 4 8 9.39 5 1 2 10.75 y 0 8 x 8 2.974 b 1.357 S 14.643 20回帰式の変形とその解
• これまで取り扱っていたデータのモデルは であった．これを次のように 変形する． この式では，未知パラメータはｂとx8である． • Excelの計算表で，ｙの予測値を計算する式に （２）式を用い，ソルバーで解いた結果が前シート の右下に示されている． y = a + bx (1) y = 8 + b(x Ä x8) (2) 21x
8の信頼区間
• 前と同じ考えで，x8の信頼区間を求めてみよう． • 次シートに示すように，曲線は左右非対称で， 放物線とはならない． • ｂの区間推定の場合と同様に，Se=29.25 の水 平線を引き，曲線との交点を求める． • 信頼区間は 1.96∼4.68 となる．これは次シート で y=8 の水平線とy-hat の信頼区間の曲線との 交点である． 22 10 15 20 25 30 35 40 45 50 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 • • x8の標準誤差はS+Vに横線を引いて，曲線との 交点を求めるという方法では求められない． • 推定値近傍を2次式で近似して，その係数から計 算する．近似の方法で区々の結果が得られる． 23• 逆推定のモデルで，残差の等高線を描く
と，非対称で，楕円とはならない．
1 2 3 4 5 6 7 8 9S1 S2 S3 S4 S5 S6 S7 S8 S9 23.00 -24.00 22.00 -23.00 21.00 -22.00 20.00 -21.00 19.00 -20.00 18.00 -19.00 17.00 -18.00 16.00 -17.00 15.00 -16.00 14.00 -15.00 24非線形回帰分析とは
• 最初に示した２つのモデル式を比較する． （１）式はパラメータa, b に関して線形（1次式）であ るが，（２）式は，展開すると，パラメータb，x8の積 が含まれ，線形ではない． • このように，パラメータに関して線形でないとき， 上に述べたように線形とは異なる現象が現れる． • このようなモデルを当てはめるのが 非線形回帰分析である， y = a + bx (1) y = 8 + b(x Ä x8) = 8 + bx Ä bx8 (2)25
非線形回帰分析
26べき関数
•
x とy の間に下記の左の関係があるとき，両辺の対 数を取ると，パラメータに関して線形となり，通常の直 線回帰分析が適用できる． • 従来は，このように線形化を工夫して解析した．y = abx; log(y) = log(a) + log(b)x
= axb; log(y) = log(a) + alog(x)
27
•
上は abx，下は axbの曲線を表わす．• 左は，ｘとｙの関係，右は ｘとlog(y) または log(x) とlog(y) の関係を表わす．凡例は a,b
0 10 20 30 40 50 60 70 0 2 4 6 8 10 10,1.2 10,0.8 1 10 100 0 2 4 6 8 10 10,1.2 10,0.8 0 20 40 60 80 100 120 0 1 2 3 4 5 10,1.5 10,0.5 0 1 10 100 1000 0.1 1 10 10,1.5 10,0.5 28 • しかし，これらの式の右辺に+cが追加されると， この方法では線形化できない． • abxで b<1 のとき，下のような曲線を当てはめ たい場合が多い． • ｘ＝０，ｘ＝∞のときの ｙをｙ０，ｙ∞で表わす． y = y1+ (y0Ä y1)bx = y1+ (y0Ä y1) exp(Bx) 0 2 4 6 8 10 0 2 4 6 8 10 2,10,0.5 10,2,0.7 29 • ここで，誤差の等分散性について考慮が必要で ある． • 対数を取ることにより，等分散性が成立する場 合が多いが，逆に等分散性が成立しなくなる場合 もありえる． • ソルバーを用いると，パラメータに関して非線形 で，対数や指数関数などの超越関数が含まれて いても，そのままでモデルを当てはめることがで きる． • 線形化と等分散性の両方を考慮して解析方法 を決める．重みつき最小2乗法を使えば適用範囲 はさらに拡大するが，今回は触れない． 30
温度と水蒸気圧の関係
• 「
理科年表」
，岩波書店 からデータを取り
出し，散布図を描く．右は水蒸気圧 y の対
数を取ったもの．
-4 -2 0 2 4 6 8 -50 0 50 100 x ln(y) 0 200 400 600 800 -50 0 50 100 x y31 • 対数変換では直線化が行 き過ぎるときは，Box-Cox変 換が用いられる． • pの初期値を設定してy* を計算する．ｘとy* の相関 係数 ｒを計算する． • ソルバーで，ｒを最大にす るp をで求める．p=0.135 となる．右のグラフで 寄与 率は 0.9999 となり，モデル の当てはまりは 極めて良い． • しかし，こうして得られた 0.135 には物理的な意味は 考えられない．
y
É=
ypÄ 1 p R2 = 0.9999 - 4 - 2 0 2 4 6 8 1 0 1 2 -50 0 50 100 x y^0.135 32 • 物理・化学現象で広く適 用されるモデルにアレニゥ スの式がある(中 上・森川訳 ， 「医薬統計学 」サイエンティスト社）． • 温度を絶対温度Tに変換 し，横軸に絶対温度のー逆 数を取ると，右のように 直 線関係が得られる． • このモデルでの傾斜は蒸 発エネルギーと結び付く． • 良く見ると，途中で折れて いるようである．y = A exp
Ä
ÄERTÅ
R2 = 0.9982 -2 -1 0 1 2 3 -4.5 -4.0 -3.5 -3.0 -2.5 1000/T ln(y) 33 • 折点の値（x0，y0）， 左側と右側の傾斜 （bL，bR）の４つの非 線形パラメータモデ ルを当てはめると， 下の結果が得られる． • x0 を絶対温度に戻 すと，1000/3.6045 =277°Kとなる． - 4 - 2 0 2 4 6 8 -4.5 -4.0 -3.5 -3.0 -2.5 -1000/K ln(p) x 0 -3.6045 y 0 1.8793 b L 6.1317 b R 5.1756 S 0.0025 34ロジスティック曲線
• 投与量（または投与量の対数）を横軸に，効果を縦 軸にとってグラフを描くとき，下のような曲線となるこ とがある． • このようなS字型の曲線を一般に生長曲線という． 生長曲線の代表的なものには(左）ロジスティック曲 線や(右）ゴンペルツ曲線がある． 0 2 4 6 8 10 0 5 10 0 2 4 6 8 10 0 5 10 35 • ロジスティック曲線は一般に次の式で表わされる． • 説明を簡単にするため，ｙの変化範囲を ０∼ｙmax とする（ymin＝０）． • （２）式で x=ｘ5 0のとき，分母=2，y=ymax/2 となる． • 一般には（１）式が用いられるが，a には物理的意 味が無いので，（２）式の方が利用しやすい．y = ymin+ ymaxÄ ymin
1 + exp(Ä(a + bx)) y = ymax 1 + exp(Ä(a + bx)) (1) = ymax 1 + exp(Äb(x Ä x50)) (2) 36 • 薬物動態の分野で良く用いられるモデルに Emaxモデルがある． 上の式を変形すると次のように ロジスティックモデ ルの式に帰着する． y = EmaxÅx ç xç+ ECç 50 y = Em ax 1 +EC50ç xç = Emax 1 + exp(ç(ln(EC50) Ä ln(x))) 0 2 4 6 8 10 0 20 40 60 0 2 4 6 8 10 0.1 1 10 100
37 • 薬物動態の分野で，用量反応関係式としては， Emaxモデルの他に，酵素反応速度論， Michaelis-Menten，受容体結合実験などがある． これらは，記号が区々で別のものと思っている人 も多いが，変換するとすべてロジスティックモデ ルである．Michaelis-Menten は b=1 としたもの である． • これらは，すべて，ソルバーを使って簡単にモ デルを当てはめることができる． • ただし，x ではなく，ln(x) に対してロジスティック 曲線が当てはまる．誤って log(x)で計算すると 係数γの値が変わるので，注意を要する． 38
ロジスティック曲線の当てはめ
• 投与量 x を0.125 から公比2の等比級数で変化 して，次シートの効果ｙが観測された． • 横軸を対数変換したグラフが示されている． • 下の式を当てはめる． y = ymax 1 + exp(ÄB(ln(x) Ä ln(x50))) 39 x y y-hat 0.125 1 0.46 0.25 3 1.75 0.5 5 6.54 1 23 22.63 2 66 63.18 4 113 118.55 8 158 153.66 16 171 166.52 32 171 170.23 64 165 171.23 ymax 171.581 x50 2.642 B 1.939 S 121.394 0 2 0 4 0 6 0 8 0 100 120 140 160 180 0.1 1.0 10.0 100.0 y-hat=ymax /(1+EXP( -B*(LN(x)-LN(x50)))) 40効力比
• 「薬剤Bが薬剤Aのｃ倍の効力を持っている」と は，「薬剤Bを x 投与した効果は，薬剤Aをｃｘ投 与した効果に等しい」ことを表わす．この関係は次 の式で表わされる． • この関係が x の値によらず成立するとき，効力 比 x が一意に定義できる． • 効力比が定義できるときには，横軸に投与量の 対数を取るとき，グラフは水平移動で重なる． yA = f (x) = g(ln(x)) yB = f (c Çx) = g(ln(c) + ln(x)) 41 • Emaxモデルのグラフに，効力比 c=2 の薬剤の 曲線（左側 ■）を加えたグラフを下に示す． 0 2 4 6 8 10 0 20 40 60 0 2 4 6 8 10 0.1 1 10 100 42効力比の例題
• 薬剤Aの投与量ｘを0.192, 0.48, 1.2, 3 と等比級 数的に変化させ，4匹の動物に投与して薬効 y を 調べる． • 薬剤Bの薬効は薬剤Aの半分程度と想像される ので，薬剤Aの2倍，すなわち，0.384, 0.98, 2.4, 6 と変化させて投与した． • 4匹の薬効の平均値を次シートのｙに示す． • ymin，ymax，b，D50 と効力比 c をパラメータと してモデルを当てはめると，次シートの右の曲線 のようにモデルがあてはまる．43 薬剤 x y y-hat A 0.192 1.08 1.06 0.480 1.48 1.51 1.200 1.71 1.70 3.000 1.72 1.73 B 0.384 0.90 0.91 0.960 1.27 1.28 2.400 1.67 1.63 6.000 1.71 1.72 c 0.3093 ymin 0.7994 ymax 1.7319 b 2.2934 D50 0.2911 S 0.0028 0.50 1.00 1.50 2.00 0.10 1.00 10.00 A B A B
byA = ymin+ ymaxÄ ymin
1 + exp(Äb(ln(x) Ä ln(x50)) byB = ymin+ ym axÄ ymin 1 + exp(Äb(ln(cx) Ä ln(x50)) 44
薬剤の相乗・拮抗効果
• 薬剤AとBを単独および併用投与して次の結果 が得られた． • A,Bの単独の効果は5,7であるのに対して，併 用した効果は11で 5+7=12よりも小さい． • 実験計画法の考えによれば，わずかではある が交互作用があり，拮抗効果が認められる． B=0 B=1 A=0 10 17 A=2 15 21 A=0 0 7 A=2 5 11 A0,B0 との差 10 12 14 16 18 20 22 B=0 B=1 A= 0 A= 2 45 • 投与量を増やして実験して，下の結果が得られた． • A=4 または B=2単独投与と，両剤を半量ずつの 併用投与を比較する． • 単独投与の平均効果(18+20)/2=19に比べて併用 投与の効果は 22 と大きい．すなわち，相乗効果が 認められる． • A=6, B=3 の併用でも，(19+22)/2=20.5 に比べて 併用投与は 25, 25 と大きく，相乗効果がある． • 前の結果と逆になった．どのように考えるか？ B=0 B=1 B=2 B=3 A=0 10 17 20 22 A=2 15 22 25 A=4 18 25 A=6 19 46 • 以下に述べる方法は，薬剤A,Bに 効力比 のモ デルが成立することが前提とされる． • 相乗・拮抗効果がない（相加性が成立する）場 合は ｙ＝ｆ（ｘA＋ｃｘB） という一般式で表わされる．ここで，ｃは効力比， ｆは任意の関数である． • この式では表わされないときは，d xA ｘB の項 を追加することにより相乗・拮抗効果を評価する ことができる．この項は dxAｘB には限らない． • 薬剤の相乗・拮抗効果は，実験計画法の交互 作用とは別のものであるという認識が必要である． 47複数の曲線の当てはめ
• 蜂に刺されると腫れる．腫れを抑制する薬剤を 開発したので，その抑制効果を評価したい． • 無投与，低用量，中用量，高用量で，浮腫量の 時間変化がどのように変るかを調べる実験を実 施した． • 各投与量毎に，複数匹の動物に対して，刺激 を与えた直後に，薬剤を投与し，0, 1, 2, 3, 4, 5 時 間後の浮腫量を測定した． • 複数の動物の浮腫量の平均値の値と，時間的 変化を表わすグラフを示す． 4849
Excelソルバーによる解析
• ４つの群毎に次の式を当てはめる． • 4群＊3＝12個のパラメータの初期値を下の青 色のセルに入力する． • 初期値を使って，予測値と残差，残差の2乗和 を求める． • さらに，2乗和の合計を黄色のセルに求める． • 2乗和の合計が最小になるパラメータをソルバー で求める．y = yinfÄ (yinfÄ y0) exp(Bx)
50 • グラフを見ると，３つの曲線は x=0 で重なってお り，上の式の y0の値はほぼ等しい． • そこで，y0の値は，薬剤の投与量によらず一定 のモデル式を当てはめる． • パラメータ入力領域を下に示す． • y0の2番目以降のセルには「＝左のセル」が入 力されている． • 着色部分のみを変化させるセルに指定して，ソ ルバーを実行すると下の結果が得られる． 51 • モデルのパラメータを3個少なくしたが，残差平 方和 はわずかに 0.01872−0.01835=0.00037 し か増加しない． • 得られた 4組の yinf とB の散布図を描いたの が，下のグラフの■マークである． • 近似直線を引き，相関係数を求めると， R2=0.9191 で 両者には深い関連がある． 52
•
yinfとB の何れか一方が4群で共通という
モデルを当てはめることを試みる．
• その方法は前と同じである．
• 何れの場合も残差平方和がかなり大きく
なり，グラフを描いても不適切なモデルであ
る．
• そこで，前のグラフの直線上に２つのパラ
メータが乗るモデルを当てはめると，かなり
良くあてはまる．
• これらの試行錯誤の末に到達したのが次
に述べるモデルである．
53 • 次のように考える． • 刺激を与えると浮腫が進むが，薬剤の効果の 表れるには時間を要する． • 最初の間は薬剤の影響は現れず，同じ傾斜で 浮腫が進む． • 前に示した曲線を表わす式で，t=0 における傾 斜は，この式をx について微分し，x に 0 を代入 して得られる． • 初期の傾斜をc とすると，B は次の式で求めら れる． dy dx= (y0Ä y1)BeBx ) (y0Ä y1)B (x = 0) B = c y0Ä y1 54 • パラメータにｃを追加 し，B は，ｃ，y0とyinfか ら計算する． • 右のグラフに示すよう に，実測点と曲線とは 良く合っているように見 える． • Se は 0.02047 で，y0 のみを共通とした残差 平方和 0.01872 に比べ ると，パラメータを3個 減らしたにもかかわら ず，わずか 0.0175 しか 増加しない． F = 0:00175=3 0:01872=(24 Ä9)ô 0:655
ロバスト推定
• 右のデータはア
ンスコムが提供し
た外れ値を含むデー
タである．
• 最小2乗法による
回帰直線は点線
のように，外れ値
に引っ張られる．
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) 56平均値と中央値
• きたないデータ（外れ値を含むデータ）では平均値 よりも中央値の方が代表値として安定した値が得ら れることが知られている． • 平均値は残差の2乗の和を最少にする最小２乗推 定値であ，中央値は残差の絶対値の和を最小とす る推定値である． • 外れ値の影響を受け難い推定がロバスト推定 • アンソコムのデータに，残差の絶対値の和を最小 とする回帰式をソルバーで当てはめた結果が次の シートである． 57 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 a 4.01 b 0.35 S 4.265 =SUM(ABS(D18:D28)) 0 2 4 6 8 10 12 14 4 9 14 58• 得られた回帰直線は，外れ値の影響を
全く受けないロバストな回帰式が得られる．
• 残差の絶対値の代わりに絶対値のｐ乗
（
1<p<2) を用いることも考えられる．
• p＝１.５とすると，次シートのように，外れ
値も幾分考慮した回帰直線が得られる．
59 x y y-hat e 4 5.39 5.29 0.104 5 5.73 5.67 0.056 6 6.08 6.06 0.018 7 6.42 6.45 -0.031 8 6.77 6.84 -0.069 9 7.11 7.23 -0.117 10 7.46 7.62 -0.156 11 7.81 8.00 -0.194 12 8.15 8.39 -0.242 13 12.74 8.78 3.960 14 8.84 9.17 -0.329 a 3.73 b 0.39 S 8.446 =SUM(ABS(D18:D28)^1.5) 0 2 4 6 8 10 12 14 4 9 14 60 • 最小ｐ乗法（１＜ｐ＜２）を，外れ値の含まれるデー タに適用することにより，ロバストな推定値を求め ることができる． • ｐが2以外のとき，平方和の微分が不連続にな るので，理論的に展開することができず，数理統 計学者は興味を持たないようである． • 最適化にGauss-Newton法（数式微分法）を用い る統計解析プログラムではこの方法は用いられ ない．数値微分が必要． • しかし，現実のデータ解析の分野では役に立つ 方法であると思われる．61
最尤法
• 誤差が正規分布に従うときは，最小2乗法が用 いられ，パラメータに関して非線形であっても，ソ ルバーで解くことができた． • 誤差が2項分布に従うときは，最小2乗法では なく，最尤法が用いられる． • この場合にも，ソルバーで解析が可能である． • まず，尤度とは何かから説明を始める． 62 • 真の有効率πが 0.2 のとき，n=10 のサンプル中 の有効数 f=3 となる確率は，2項分布で計算され る． • ｆを0 ∼ 10 に変えたときの確率が次シートの π=0.2 の列に求められている． • それをグラフ化したのが，左下のグラフである． • ｆは整数値しか取らないので，グラフは棒グラフ となる． 63 64 • 逆に，ｆ＝ ３が得られたとき，この結果はπ=0.2 のとき どの位起こり易いかを考える． • f=3の行を横に見ると0.20 である． • π=0.5のとき，この値は 0.12 となる． • πの値を変化させたときこの値の変化を見た のが右下のグラフである．πは連続的に変化す るので，滑らかな曲線で表わされる． • このとき，縦軸の値を尤度と呼ぶ． • 尤度が最大となるπは 0.3 となる． • これが，最尤法の考え方で，得られた結果は 最尤推定値と呼ばれる． 65ロジスティック回帰分析
• 薬剤の投与量 ｘを等比級数で変化して，それ ぞれの投与量で n=10 匹に投与し，効果の見ら れた匹数 ｆを求めた． • その結果を次シートの左に示す． • グラフの横軸は投与量ｘの対数表示，縦軸は 有効率 p=f/nである． • 両者の関係はロジスティック曲線で表わされる と思われる． • ソルバーで解いた結果が次シートに示されてい る． 66 x f n p p-hat -2ln(L) p-hat -2ln(L) 101 0 10 0.0 0.042 0.868 0.042 0.868 136 2 10 0.2 0.164 2.483 0.164 2.483 183 5 10 0.5 0.465 2.854 0.465 2.854 247 8 10 0.8 0.796 2.396 0.796 2.396 333 9 10 0.9 0.946 2.223 0.946 2.224 450 10 10 1.0 0.987 0.254 0.987 0.254 a -26.211 x50 188.247 b 5.004 b 5.004 L 11.078 L 11.078 0.0 0.2 0.4 0.6 0.8 1.0 100 1000 -2ln(L)=-2*LN(BINOMDIST(f, n ,p-hat,FALSE)) p-hat=1/(1+EXP(-(a+b*LN(x)))) p-hat=1/(1+EXP(-b*(LN(x)-LN(x50))))67
