ロジスティック回帰

教科書のデータ

表7.2 かぶと虫の死亡データ ガス濃度はlog10CS2mgl-1

ガス濃度 かぶと虫の数 死亡数 x n y 1.6907 59 6 1.7242 60 13 1.7552 62 18 1.7842 56 28 1.8113 63 52 1.8369 59 53 1.861 62 61 1.8839 60 60 出典

ロジスティック回帰

表7.2から、ガス濃度毎の死亡割合を計算して図にしました。

このデータに基づいて、死亡確率をガス濃度の関数として表したいと思います。但し確率は0から1の間ですので、直線や多項式を用いると、ガス濃度によっては0から1の範囲をはみ出てしまうので不適です。そこで0から1の範囲をはみ出さないモデルとして次の3つがしばしば使われます。
表の上からi番目の行は、ガス濃度がx_iであるときに、かぶと虫n_i匹の中のx_i匹が死んだということを表しますので、その濃度で各々のかぶと虫が死ぬ確率をp_iとすると、死亡数の確率分布は二項分布B(n_i,p_i)になります。このp_iをx_iのどのような関数で表すか、という問題を考えます。二項分布の確率関数p(x)=_nC_xp^x(1-p)^n-xで、Excelでは_nC_x=COMBIN(n,x)です。実はこの項にはp(x)が含まれないので省いも構いません。

モデル1:ロジスティックモデル

死亡確率はp(x)=exp(α+βx_i)/(1+exp(α+βx_i)))である。

モデル2:プロビットモデル

死亡確率はp(x)=Φ((x_i-μ)/σ)である。但しΦは標準正規分布の分布関数で、Excelでは Φ((x-μ)/σ))=NORM.DIST(x,μ,σ,1)で書けます。

モデル3:極値モデル

死亡確率はp(x)=1-exp(-exp(α+βx_i))である。

まずロジスティックモデルで確率を計算しましょう。α、βは最初は分からないので、xの係数βは0にして、αはp(x)=1/2になるようにα=0にします。

次に対数尤度の計算です。二項分布B(n,p)の確率関数は_nC_xp^x(n-p)^n-xでした。対数尤度はlog_nC_x+x_ilog p_i+(n_i-x_i)log (1-p_i)をiを1から8まで全部足したものになります。

F12セルの対数尤度を最大にするようなF1、F2セルの値を、ソルバーを用いて求めてください。

モデルによって推定した死亡確率と、実際の死亡割合を比較するとこのようになります。

点は実際の死亡割合、曲線がモデルです。

プロビットモデルに関しても同様にまず確率を求めて、それから対数尤度を計算してみましょう。μは最初はデータの中ごろということで1.8、σはとりあえず1にしてソルバーを使ってください。

すると、分母のσが小さくなりすぎてエラーになりました。

エラーにならないように変数の範囲に制約をつけましょう。

σ=0.02…でエラーになっているので、σ≧0.03にしましょう。制約条件(U)の追加(A)をクリックして、制約式を入力します。

もう一度ソルバーを実行してみます。

今度はエラーになりませんでしたし、σ=0.05…ですので制約の端でもありません。もしもソルバーの解が制約の端である0.03と一致していたら、ソルバーの結果が制約式に依存してしまっているので、制約式を少し広げてσ≧0.025などを試すべきですが、今回は端ではないのでこれで大丈夫です。

最後に極値モデルに関しても最尤推定量を求めてください。α、βの初期値はロジスティックモデル同様α=β=0にします。今度もソルバーでエラーが出ました。

今回は死亡確率が1になっているのは最後の行だけで、ここは60匹全員が死亡しているので尤度としては1⁶⁰で良いのですが、公式にそのまま代入すると_nC_xp^x(1-p)^n-xの中の(1-p)^n-xの部分が0⁰となり、数学的には1なのですがExcelは底が0だとエラーになってしまいます。エクセルがエラーを起こす(1-p)^n-xの部分は、最後の行ではどうせ1を掛けているだけなので削除してしまっても計算式は変わりません。

今回はソルバーで計算することが出来ました。

レポート課題

それぞれのモデルの最尤推定量を求めて、どのモデルが一番良いか選択したエクセルファイルを提出してください。