2017-08-06

ロジスティック回帰

教科書のデータ

表7.2 かぶと虫の死亡データ
ガス濃度はlog10CS2mgl-1
ガス濃度かぶと虫の数死亡数
xny
1.6907596
1.72426013
1.75526218
1.78425628
1.81136352
1.83695953
1.8616261
1.88396060


出典

ロジスティック回帰

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

このデータに基づいて、死亡確率をガス濃度の関数として表したいと思います。但し確率は0から1の間ですので、直線や多項式を用いると、ガス濃度によっては0から1の範囲をはみ出てしまうので不適です。そこで0から1の範囲をはみ出さないモデルとして次の3つがしばしば使われます。
表の上からi番目の行は、ガス濃度がxiであるときに、かぶと虫ni匹の中のxi匹が死んだということを表しますので、その濃度で各々のかぶと虫が死ぬ確率をpiとすると、死亡数の確率分布は二項分布B(ni,pi)になります。このpiをxiのどのような関数で表すか、という問題を考えます。二項分布の確率関数p(x)=nCxpx(1-p)n-xで、ExcelではnCx=COMBIN(n,x)です。実はこの項にはp(x)が含まれないので省いも構いません。
モデル1:ロジスティックモデル
死亡確率はp(x)=exp(α+βxi)/(1+exp(α+βxi)))である。
モデル2:プロビットモデル
死亡確率はp(x)=Φ((xi-μ)/σ)である。但しΦは標準正規分布の分布関数で、Excelでは Φ((x-μ)/σ))=NORM.DIST(x,μ,σ,1)で書けます。
モデル3:極値モデル
死亡確率はp(x)=1-exp(-exp(α+βxi))である。
まずロジスティックモデルで確率を計算しましょう。α、βは最初は分からないので、xの係数βは0にして、αはp(x)=1/2になるようにα=0にします。

次に対数尤度の計算です。二項分布B(n,p)の確率関数はnCxpx(n-p)n-xでした。対数尤度はlognCx+xilog pi+(ni-xi)log (1-pi)をiを1から8まで全部足したものになります。

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

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

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

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

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

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

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

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

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

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

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

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

レポート課題

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