静粛に、只今統計勉強中

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

Excel VBAでクロス表の残差分析ができる関数を作ってみた

前回、以下の表についてカイ二乗検定を行い、クラメールのVを計算して、学年と満足度の間には何らかの関連があるらしいことがわかりました。*1
しかし、具体的にどのような関連があるのでしょうか?

f:id:cyclo-commuter:20171209101909j:plain

このような疑問を持ったときに行なうのがクロス集計表の残差分析です。
今回は、Exel VBAで関数を作って残差分析を行います。

クロス集計表の残差分析は、表中の各セルについて調整済み標準化残差を計算して行われる由。「調整済み」「標準化」「残差」とはなんぞ?
前回までは、いきなりソースコードをベタッと貼り付けていましたが、理解を深めるため、手順を一つ一つ追ってみようと思います。以下、B列を例に。

まず、残差。これは「観測値-期待値」のこと。B3セルの残差をr1、以下r2、r3として、r1 = 10-(22*60/150) = 1.2 、r2 = 3-(22*50/150) = -4.33...、r3 = 9-(22*40/150) = 3.13...となりますね。r1 + r2 + r3 = 0 となるのを覚えておきます。

次に、標準化残差。残差を期待値の平方根で割って標準化します。はて?

f:id:cyclo-commuter:20171212102047j:plain
なぜ期待値の平方根で割ると標準化されるのでしょう? ネットで調べると「期待値の平方根は残差の標準偏差とありますが、何故そうなるかの解説が見当たりません。ここが統計学を勉強していて辛いところで、わかっている人は当然のこととして計算を進めていくので、初学者はたいてい置いてきぼりです。やむを得ないので割り切って計算を進めますが、概念としては、残差はもともとバラバラの度数を元に計算されていてセル同士の比較ができないので、比較できるように統一基準に合わせる(=標準化する)、と理解します。B3セルの標準化残差をsr1、以下sr2、sr3として、
sr1 = 1.2 / √(22*60/150) = 0.40、sr2 = -4.33 / √(22*50/150) = -1.60、sr3 = 3.13 / √ (22*40/150) = 1.29...となりますね。sr1 + sr2 + sr3 ≠ 0 になりました。

最後に、調整済み標準化残差。標準化残差をさらに[標準化残差の分散の平方根]で割ったのが、調整済み標準化残差。こうすると、2×j表の縦計が 0 になるそうです。でも、i×j表(i>2)では、縦計は 0 になりません。それでも調整はするみたいです。

f:id:cyclo-commuter:20171212132409j:plain
この式も何の脈絡もなく「こうである」と来るのですが、なぜそうなるのかはわかりません。標準化残差は近似的に正規分布に従うそうです(これはなんとなくわかります)ので、もしかすると正規分布の定義から導出できるのかもしれません。私には無理ですが。ともかく、式を完成させます。
f:id:cyclo-commuter:20171212134751j:plain
B3セルの調整済み標準化残差をasr1、以下asr2、asr3として、
asr1 = 1.2 / √(22*60/150)(1-22/150)(1-60/150) = 0.57、asr2 = -4.33 / √(22*50/150)(1-22/150)(1-50/150) = -2.12、asr3 = 3.13 / √ (22*40/150)(1-22/150)(1-40/150) = 1.64 となりました。

ところどころ式の意味がわからないのに、計算それ自体は中学までの算数・数学でできてしまうのだから、統計学ってたちが悪いです。結局、たいして理解は深まりませんでした・・・

関数のソースコードは以下のとおりです。*2

'クロス表の残差分析(調整済み標準化残差z値または両側p値を返す)
Function Resid_Anal(実測値範囲, 参照, 種類 As Boolean)

    Dim o_Row As Integer '行数
    Dim o_Clm As Integer '列数
    Dim o_R_Sum() As Long '行合計
    Dim o_C_Sum() As Long '列合計
    Dim o_Sum As Long 'N
    Dim i As Integer
    Dim j As Integer
    Dim o As Double '実測値
    Dim e As Double '期待値

    o_Row = 実測値範囲.Rows.Count
    o_Clm = 実測値範囲.Columns.Count
    ReDim o_R_Sum(1 To o_Row)
    ReDim o_C_Sum(1 To o_Clm)

    For i = 1 To o_Row
        For j = 1 To o_Clm
            o_R_Sum(i) = o_R_Sum(i) + 実測値範囲(1).Offset(i - 1, j - 1)
            o_C_Sum(j) = o_C_Sum(j) + 実測値範囲(1).Offset(i - 1, j - 1)
            o_Sum = o_Sum + 実測値範囲(1).Offset(i - 1, j - 1)
        Next
    Next

    i = 参照.Row - 実測値範囲(1).Row
    j = 参照.Column - 実測値範囲(1).Column
    o = 実測値範囲(1).Offset(i, j)
    e = o_C_Sum(j + 1) * o_R_Sum(i + 1) / o_Sum
    If 種類 = True Then
        Resid_Anal = (1 - WorksheetFunction.Norm_S_Dist(Abs( (o - e) / Sqr(e * (1 - o_C_Sum(j + 1) / o_Sum) * (1 - o_R_Sum(i + 1) / o_Sum))), True)) * 2
    Else
        Resid_Anal = (o - e) / Sqr(e * (1 - o_C_Sum(j + 1) / o_Sum) * (1 - o_R_Sum(i + 1) / o_Sum))
    End If

End Function

引数が3つになったので、わかりやすいように日本語名称にしました。
 実測値範囲:クロス表のラベルと計を除く範囲
 参照   :計算対象となるセル
 種類   :True  両側検定のp値を返す False  z値を返す
 入力例  :B13セルには、=Resid_Anal($B$3:$F$5,B3,FALSE) と入力

f:id:cyclo-commuter:20171212141703j:plain

見出しを付け忘れましたが、上(12:15行)が調整済み標準化残差の表、下(17:20行)が有意確率(両側)の表です。
セルのフォントカラーと背景色は、条件付き書式を使って有意水準0.05で色分けしました。
結果を見ると、2年生は他の学年よりも[大変満足]が少なく、[少し不満]が多い。また、3年生は他の学年よりも[だいたい満足]が多いことがわかります。

どうでしょう、だいぶ分析っぽくなってきたのではないでしょうか。
数ある統計分析ソフトは皆<コマンド→結果出力>タイプですが、Excel上で完結できるユーザー定義関数って、重宝する場面は多いと思うのです。

*1:くどいようですが、データは架空のものです。

*2:2017/12/23 コード個所をpre記法に直しました