静粛に、只今統計勉強中

仕事でデータ分析をすることになったバリバリ文系アラフィフのおっさんが、独学で統計の勉強を始めました。

Excel VBAで『マンガでわかる統計学 回帰分析編』のロジスティック回帰分析の手順をマクロにしてみた2

『マンガでわかる統計学 回帰分析編』で解説されているロジスティック回帰分析の手順を自動化するマクロを作ったので、前回からその説明をしています。

今回はその続きです。

 

前回うっかり触れませんでしたが、ロジスティック回帰分析の回帰式は、

{y=\dfrac{1}{1+e^{-(a_{1}x_{1}+a_{2}x_{2}+・・・+a_{p}x_{p}+b)}} }

だそうです。*1
もう、見た瞬間に文系脳が受付を拒否してしまうんですが、とりあえず分母と分子の「1」に着目しておきましょう。1を1+α(α≧0)で割ることで解が0と1の間に収まるようになっているんですね。で、eネイピア数 2.718・・・だから、これを何乗しても0以上になる、と。だから、確率が求まるんですね。
式をグラフにすると、
f:id:cyclo-commuter:20180214093947g:plain
こんなふうになるのだそうです。実はこれ、正規分布の累積分布関数のグラフなんですが、よく似てます。似てることにも意味あるのかも知れませんが、とりあえず置いておいてソースコードの説明に戻りましょう。 

 

'ソルバーで対数尤度関数の最大値を求める
SolverOk SetCell:=L, MaxMinVal:=1, ValueOf:=0, ByChange:=a(2, 1) & ":" & b(2), _
Engine:=1, EngineDesc:="GRG Nonlinear"
SolverSolve userfinish:=True

 目的セルに対数尤度関数のアドレスを、変数セルに回帰係数から定数までのアドレスを指定して、ソルバーに最尤推定値を求めてもらいます。

 

With Application.WorksheetFunction
n1 = .CountIf(y, 1) '事象が起きた回数
n2 = .CountIf(y, 0) '事象が起きなかった回数
k = n1 * .Ln(n1) + n2 * .Ln(n2) - (n) * .Ln(n)
'寄与率
Rng.Offset(c + 7, c + 7).Value = 1 - Range(L).Value / k
'判別的中率
Rng.Offset(c + 8, c + 7).Value = (.CountIfs(y, 1, haty, ">=0.5") + .CountIfs(y, 0, haty, "<0.5")) / n
'尤度比検定
Rng.Offset(c + 9, c + 7).Value = .ChiSq_Dist_RT(2 * (Range(L).Value - k), c)
'オッズ比
For i = 1 To c
Rng.Offset(i + 4, c + 8).Value = Exp(Rng.Offset(i + 4, c + 7).Value)
Next

寄与率

{1-\dfrac{対数尤度関数Lの最大値}{k}}

判別的中率は、予測値0.5以上を「事象が起きる」として、観測値と予測値が一致した回数を観測数で割ったもの。

尤度比検定は、統計検定量2×{対数尤度関数Lの最大値-k}を自由度:説明変数の数でカイ二乗検定した上側確率(p)を求めています。

オッズ比は、ネイピア数を回帰係数乗した「調整されたオッズ比」です。*2 

 

'Wald検定
ReDim buf1(1 To c + 1, 1 To n)
For j = 1 To n
buf1(c + 1, j) = haty(j) * (1 - haty(j))
For i = 1 To c
buf1(i, j) = x(j, i) * buf1(c + 1, j)
Next
Next
buf2 = .MMult(buf1, x)
buf3 = .MInverse(buf2)
For i = 1 To c + 1
Rng.Offset(i + 4, c + 9).Value = .ChiSq_Dist_RT(Rng.Offset(i + 4, c + 7) ^ 2 / buf3(i, i), 1)
Next

実は、ここで行っている行列のかけ算は『マンガでわかる統計学 回帰分析編』に従ったものではありません。

こちらのサイトの計算方法のほうがシンプルで扱いやすかったからです。
ただし、Transposeの手間を省くため、x1×y', x2×y', ... , xc, y'=ŷ×(1-ŷ) の行列を最初からタテ・ヨコ逆に変数buf1に代入しています。
それから、(株)アイスタットさんの説明では、y'=ŷ×(1-ŷ) って具合にマイナスが付いてますが、試してみたらマイナスにする必要はありませんでした。標準誤差を求める数式が (-1*C26)^0.5 ってなってるんですが、正負を戻すために -1* をするなら最初から y' を負の値にしなければいいのに、と思います。
もうひとつ、(株)アイスタットさんの手順では標準誤差を使って (回帰係数÷標準誤差)2 という数式で検定統計量を求めていますが、逆行列の時点で『マンガでわかる統計学 回帰分析編』の手順とピッタリ同じ値ですので、標準誤差は求めず 回帰係数2÷逆行列の対角要素 で検定統計量を計算しています。

 

以上、細かな書式設定の部分は割愛させていただきましたが、ソースコードの説明は終わりです。次回はいよいよ使い方と実行結果を説明したいと思います。 

*1:『マンガでわかる統計学 回帰分析編』154ページ

*2:『マンガでわかる統計学 回帰分析編』190ページ