Excel VBAで『マンガでわかる統計学 回帰分析編』のロジスティック回帰分析の手順をマクロにしてみた３
前回までの２回でロジスティック回帰分析のマクロのソースコードについて説明しました。今回は、ソースコード全文を掲載して、使い方と実行結果について説明したいと思います。
ソースコード
'ロジスティック回帰分析を行なう
Sub Logistic_Regression()
Application.ScreenUpdating = False
Dim a, b(1 To 2), c, i, j, k, L, n, x
Dim eq, n1, n2
Dim buf1, buf2, buf3
Dim Rng As Range, y As Range, haty As Range
c = Selection.Columns.Count - 1 '説明変数の数
n = Selection.Rows.Count - 1 '観測数
ReDim a(1 To 2, 1 To c) '回帰係数のアドレス取得用
Set Rng = Selection(1) '基準セル
With Rng
Set y = .Offset(1, c).Resize(n) '目的変数
Columns(.Column + c).Insert 'ダミーの列を追加
.Offset(1, c).Resize(n).Value = 1 'ダミー列に「１」を入力
x = .Offset(1).Resize(n, c + 1) '説明変数+ダミー
Columns(.Column + c).Delete 'ダミーの列を削除
.Offset(, c + 1).ColumnWidth = 1
.Offset(, c + 3).ColumnWidth = 1
.Offset(, c + 5).ColumnWidth = 2
Columns(.Column).Copy Columns(.Column + c + 2)
Columns(.Column).Copy Columns(.Column + c + 4)
.Offset(1, c + 6).Value = "尤度関数" '尤度関数のラベル
.Offset(1, c + 7).FormulaR1C1 = "=PRODUCT(C[-3])" '尤度関数の数式
.Offset(2, c + 6).Value = "対数尤度関数" '対数尤度関数のラベル
With .Offset(2, c + 7)
.FormulaR1C1 = "=LN(R[-1]C)" '対数尤度関数の数式
L = .Address '対数尤度関数のアドレス
End With
With .Offset(1, c + 6).Resize(2, 2) '罫線
.Borders.LineStyle = True
.Borders(xlInsideHorizontal).LineStyle = False
End With
.Offset(4, c + 7).Value = "回帰係数"
.Offset(4, c + 8).Value = "オッズ比"
.Offset(4, c + 9).Value = "Wald検定"
For i = 1 To c
.Offset(i + 4, c + 6).Value = .Offset(, i - 1).Value '説明変数のラベル
a(1, i) = .Offset(i + 4, c + 7).Address(ReferenceStyle:=xlR1C1) '回帰係数のアドレス
a(2, i) = .Offset(i + 4, c + 7).Address '回帰係数のアドレス
'回帰式の組み立て
eq = eq & a(1, i) & "*" & .Offset(c + i + 11, c + 7).Address(ReferenceStyle:=xlR1C1) & "+"
Next
.Offset(c + 5, c + 6).Value = "定数" '定数項のラベル
With .Offset(4, c + 6).Resize(c + 2, 4) '罫線
.Borders.LineStyle = True
.Borders(xlInsideHorizontal).LineStyle = False
.Resize(1).Borders.LineStyle = True
End With
b(1) = .Offset(c + 5, c + 7).Address(ReferenceStyle:=xlR1C1) '定数項のアドレス
b(2) = .Offset(c + 5, c + 7).Address '定数項のアドレス
eq = eq & b(1) '回帰式の組み立て
.Offset(c + 7, c + 6).Value = "寄与率"
.Offset(c + 8, c + 6).Value = "判別的中率"
.Offset(c + 9, c + 6).Value = "尤度比検定"
With .Offset(c + 7, c + 6).Resize(3, 2) '罫線
.Borders.LineStyle = True
.Borders(xlInsideHorizontal).LineStyle = False
End With
.Offset(c + 11, c + 6).Value = "予測"
For i = 1 To c
.Offset(c + i + 11, c + 6).Value = .Offset(, i - 1).Value '説明変数のラベル
Next
.Offset(c * 2 + 12, c + 6).Value = "確率"
With .Offset(c * 2 + 12, c + 7)
.FormulaR1C1 = "=1/(1+EXP(-(" & eq & ")))" '予測用の回帰式
.NumberFormatLocal = "0.0%"
End With
With .Offset(c + 12, c + 6).Resize(c, 2) '罫線
.Borders.LineStyle = True
.Borders(xlInsideHorizontal).LineStyle = False
With .Offset(c).Resize(1).Borders
.LineStyle = True
.Weight = xlMedium
End With
.Resize(c + 1).Borders(xlInsideVertical).Weight = xlThin
End With
.Offset(, c + 6).Resize(, 4).EntireColumn.AutoFit
.Offset(, c + 2).Value = "予測値"
eq = "" '回帰式をリセット
For i = 1 To c
eq = eq & a(1, i) & "*RC[-" & c - i + 3 & "]+" '回帰式の組み立て
Next
eq = eq & b(1) '回帰式の組み立て
Set haty = .Offset(1, c + 2).Resize(n) '予測値のセル範囲
haty.FormulaR1C1 = "=1/(1+EXP(-(" & eq & ")))" '回帰式を代入
.Offset(, c + 4).Value = "（尤度関数の計算過程）"
.Offset(1, c + 4).Resize(n).FormulaR1C1 = "=IF(RC[-4]=1,RC[-2],1-RC[-2])"
End With
'ソルバーで対数尤度関数の最大値を求める
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
'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
End With
ActiveWindow.DisplayGridlines = False
Application.Calculation = xlAutomatic
Application.ScreenUpdating = True
End Sub
使い方と実行結果
使い方は、毎度おなじみ、ソースコード全行をＶＢＥで標準モジュールにコピペしたら、データ範囲を選択して、［開発］→［マクロ］→［Logistic_Regression］を選択→【実行】をクリックするだけ、です。
使用上の注意が４つあります。
- 必ず表にラベル行を作っておいてください。
- 必ず説明変数を左に、目的変数を一番右に配置してください。
- 必ずラベル行からセル範囲を選択してください。
- あらかじめソルバーを使える状態にしておいてください。
クリック後１～数秒で以下のように結果が出力されます。
出力について一つずつ確認してみましょう。
回帰係数と定数です。この値を使って予測値と予測の確率を計算しています。
試しに、日曜日の最高気温２３℃でノルンスペシャルが売れる確率を計算してみると、
44.1％と出ました。日曜日であっても、最高気温２３℃では売れないようですね。
オッズ比について、『マンガでわかる統計学 回帰分析編』には計算の仕方は書かれていたのですが、値の意味の説明はありませんでした。
なので、例によってグーグル先生に頼ってみたら、
『オッズ比』について – 原因と結果の関係性，傾向を知る方法funmatu.wordpress.com
この ２サイトの説明がとても分かりやすかったです。ただし、順序尺度と言っても
説明変数のデータ単位が全て同じ場合はオッズ比は寄与順位に適用できます。
データ単位が異なる場合、重回帰分析同様、回帰係数の単純比較はできません。
（重回帰の場合は標準回帰係数で寄与順位をします。）
Wald統計量は（回帰係数÷標準誤差）の２乗で検定統計量です。
標準誤算で割ることによって基準化され、寄与順位に適用できます。
だそうですので、そのままじゃ使えない。正直オッズ比、扱いづらいです。
それに引き換え、マルチに使えるのがWald検定のｐ値。
（寄与率は）精度が高いほど１に近づき、そうじゃないときほど０に近づくのよ。
ロジスティック回帰式の寄与率は大きくなりづらい傾向があるといわれているから、参考程度にしておいたほうがいいかもね。
『マンガでわかる統計学 回帰分析編』１７１ページ
との由。で、もう一つの方法として誤判別率についての説明があるのですが、私の方はこれを裏返して判別的中率としてみました。
尤度比検定は、「回帰係数を包括的に検討する検定」（１７５ページ）だそうです。こちらも、お好みの有意水準でどうぞ。
以上で説明終了です。
WEBで検索していると、まだまだアレコレ足せそうですが、今回はここまで。
