静粛に、只今統計勉強中

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

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

『マンガでわかる統計学 回帰分析編』の最終章はロジスティック回帰分析について解説されています。ロジスティック回帰分析は、ある事象が起こる「確率」を予測するための分析手法で、天気予報の発雷確率を求めるためなんかにも使われているようです。

やれ指数だ対数だ、オッズだロジットだと、いよいよもって難解になってます。
ページ数の関係からか、説明も若干飛ばし気味の印象。
エクセルExcelで計算過程を追っていくにもこの本だけでは間に合わず、ちょいちょいグーグル先生のお世話にならなくてはなりませんでした。が、なんとか計算できたのでVBAでマクロを作ってみました。

分析手順のおさらいを兼ねて、順を追ってソースコードの説明をしてみたいと思います。

 

c = Selection.Columns.Count - 1 '説明変数の数
n = Selection.Rows.Count - 1 '観測数
ReDim a(1 To 2, 1 To c) '回帰係数のアドレス取得用
Set Rng = Selection(1) '基準セル

1行目の選択セル範囲の列数-1の「-1」は、選択範囲に目的変数の列があるからです。
2行目の「-1」はラベル行の分を減じています。
3行目の配列変数 a は、
 ・ソルバーで回帰係数のセルのアドレスを指定するとき
 ・予測用の回帰式を数式でセルに入力するとき
 ・予測値計算用の回帰式を数式でセルに入力するとき
に使います。
4行目で選択範囲の左上のセルを、以後の処理でoffsetするための基準にしています。

 

With Rng
Set y = .Offset(1, c).Resize(n) '目的変数
Columns(.Column + c).Insert 'ダミーの列を追加
.Offset(1, c).Resize(n).Value = 1 'ダミー列に「1」を入力
x = .Offset(1).Resize(n, c + 1) '説明変数+ダミー
Columns(.Column + c).Delete 'ダミーの列を削除

2行目、目的変数のセル範囲を変数 y に代入します。Rangeオブジェクト変数にしたのは、 後でCountIf関数の引数にするためです。
3~6行目で、Wald検定に必要な「説明変数+ダミーの観測値」を配列変数 x に代入しています。ループを回すのが面倒だったのでこんな具合にしてみましたが、邪道かも。

 

.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

ソルバーを使うための数式をセルに入力しています。
尤度関数は、観測数をn、事象が起きる確率をP、事象が起きた回数をaとしたとき、
Pa(1-P)n-a で表されるようです。*1 ここでは、後述する(尤度関数の計算過程)の値をPRODUCT関数で掛け合わせて計算しているんですね。
対数尤度関数は log{ Pa(1-P)n-a } なので、LN関数で尤度関数の自然対数を求めています。

 

.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) '回帰式の組み立て

 回帰係数、オッズ比、Wald検定の値は説明変数ごとに求められるのでひとまとめにして、まだ値は求まらないので、ここでは入れ物だけ用意しています。
ついでに変数abにアドレス(文字列)を代入して、変数eqで予測用の数式を組み立てています。

 

.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 + 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])"

 ソルバー用の(尤度関数の計算過程)の数式と、そのための予測値を求める数式を入力しています。(尤度関数の計算過程)の数式は観測値の事象が起きていればP、起きていなければ1-Pが戻りますので、これをかけ合わせれば尤度関数が求まります。

 

少々長くなったので、続きは次回。

*1:『マンガでわかる統計学 回帰分析編』157ページの説明を一般化してみました。