Excel VBAで散布図に信頼区間と予測区間を表示するマクロを作ってみた１

実は、前回まで４回をかけて行なった『マンガでわかる統計学　回帰分析編』に従って単回帰分析の手順を追ってみたシリーズは、今回のための下準備でした。

これまでの計算を全てVBAで自動化した上で、信頼区間＆予測区間付きの散布図も作っちゃえ！　というのが、本当にやりたかったことなのです。*1

あらかじめワークシート上で計算したおかげで、サクサク実装できました。
ので、いきなりソースコードから。

ソースコード 

'単回帰の信頼区間・予測区間付きで散布図を作成する
Sub Scatter_Plot()

    Const Ttl = "散布図 with 信頼区間＆予測区間"
    Dim Rng, y, x, n, a, b, xbar, Sxx, Se, F, haty, CI, PI, i
    
    If Selection.Columns.Count <> 2 Then
        MsgBox "独立変数（x）、従属変数（y）の順に２列の" & vbCrLf & _
            "セル範囲を選択してから実行してください。", vbCritical, Ttl
        Exit Sub
    ElseIf Selection(1).Row = 1 Then
        MsgBox "ラベル行は選択しないでください。", vbCritical, Ttl
        Exit Sub
    End If
    
    Application.ScreenUpdating = False
    
    Set Rng = Selection.Resize(1, 1).Offset(-1, 2)  '基準セル
    Set y = Selection.Resize(, 1).Offset(, 1)       '従属変数
    Set x = Selection.Resize(, 1)                   '独立変数
    n = x.Rows.Count                                '観測数
    ReDim haty(1 To n)
    ReDim CI(1 To n)
    ReDim PI(1 To n)
    With Application.WorksheetFunction
        a = .Slope(y, x)                        '回帰係数
        b = .Intercept(y, x)                    '切片
        xbar = .Average(x)                      'xの平均値
        Sxx = .DevSq(x)                         'xの偏差平方和
        For i = 1 To n
            haty(i) = a * x(i, 1) + b           'yの予測値
        Next
        Se = .SumXMY2(y, haty)                  '残差平方和
        
        'Ｆ検定
        F = (a ^ 2 * Sxx) / (Se / (n - 2))
        With Rng.Offset(1, 6)
            .Value = "回帰係数の検定"
            .EntireColumn.AutoFit
        End With
        Rng.Offset(1, 7).Value = .F_Dist_RT(F, 1, n - 2)
        
        F = .F_Inv_RT(0.05, 1, n - 2)   '区間推定用のＦ値
    End With
        
    '信頼区間と予測区間と標準化残差
    With Rng
        .Value = "信頼区間" & vbLf & "下限95%"
        .Offset(, 1).Value = "信頼区間" & vbLf & "上限95%"
        .Offset(, 2).Value = "予測区間" & vbLf & "下限95%"
        .Offset(, 3).Value = "予測区間" & vbLf & "上限95%"
        .Offset(, 4).Value = "標準化" & vbLf & "残差"
        For i = 1 To n
            '信頼区間の幅
            CI(i) = Sqr(F * (1 / n + (x(i, 1) - xbar) ^ 2 / Sxx) * Se / (n - 2))
            '予測区間の幅
            PI(i) = Sqr(F * (1 + 1 / n + (x(i, 1) - xbar) ^ 2 / Sxx) * Se / (n - 2))
            .Offset(i).Value = haty(i) - CI(i)       '信頼区間下限
            .Offset(i, 1).Value = haty(i) + CI(i)    '信頼区間上限
            .Offset(i, 2).Value = haty(i) - PI(i)    '予測区間下限
            .Offset(i, 3).Value = haty(i) + PI(i)    '予測区間上限
            With .Offset(i, 4)
                .Value = (y(i, 1) - haty(i)) / Sqr(Se / (n - 2)) '標準化残差
                '標準化残差の絶対値が３を超える場合は警告
                If Abs(.Value) > 3 Then .Font.Color = vbRed
            End With
        Next
        '書式を整える
        With .Resize(, 5)
            .Interior.Color = Rng.Offset(, -1).Interior.Color
            .Font.Color = Rng.Offset(, -1).Font.Color
            .HorizontalAlignment = xlCenter
            .Resize(n + 1).Borders.LineStyle = True
        End With
    End With
    
    '散布図作成
    ActiveSheet.Shapes.AddChart2(240, xlXYScatter).Select
    With ActiveChart
        .SetSourceData Source:=Range("'" & ActiveSheet.Name & "'!" & _
            Rng.Offset(, -2).Resize(n + 1, 6).Address)
        .HasTitle = False   'グラフタイトル除去
        .HasLegend = False  '凡例除去
        '横軸の最小値
        .Axes(xlCategory).MinimumScale = Round(WorksheetFunction.Min(x) - 1, 0)
        '縦軸の最小値
        .Axes(xlValue).MinimumScale = _
            Round(WorksheetFunction.Min(Rng.Offset(1, 2).Resize(n)) - 1, 0)
        '回帰直線を追加
        With .FullSeriesCollection(1)
            .Trendlines.Add
            With .Trendlines(1)
                .DisplayEquation = True '回帰式
                .DisplayRSquared = True '寄与率
                With .DataLabel
                    .Left = 29
                    .Top = 11
                End With
                '書式を整える
                With .Format.Line
                    .ForeColor.RGB = RGB(192, 0, 0)
                    .Weight = 1.25
                End With
            End With
        End With
        '信頼区間と予測期間の近似曲線を追加
        For i = 2 To 5
            With .FullSeriesCollection(i)
                'マーカー非表示
                .Format.Fill.Visible = msoFalse
                .MarkerForegroundColorIndex = xlColorIndexNone
                .Trendlines.Add
                With .Trendlines(1)
                    .Type = xlPolynomial '多項式
                    .Order = 2           '２次
                    '書式を整える
                    With .Format.Line
                        .Weight = 1.25
                        Select Case i
                        Case 3, 5
                            .ForeColor.ObjectThemeColor = i + 3
                        End Select
                    End With
                End With
            End With
        Next
        .Parent.Top = Rng.Offset(3, 6).Top      '位置調整（上端）
        .Parent.Left = Rng.Offset(3, 6).Left    '位置調整（左端）
    End With
    
    ActiveWindow.DisplayGridlines = False       'グリッド除去
    Application.ScreenUpdating = True

End Sub

コードの解説

前半の信頼区間と予測区間を出力するところは、ワークシート上の計算をＶＢＡに置き換えただけですので、説明は省かせていただきます。

前４回では触れなかったのがコチラ↓

Rng.Offset(, 4).Value = "標準化" & vbLf & "残差"
For i = 1 To n
    With Rng.Offset(i, 4)
        .Value = (y(i, 1) - haty(i)) / Sqr(Se / (n - 2)) '標準化残差
        '標準化残差の絶対値が３を超える場合は警告
        If Abs(.Value) > 3 Then .Font.Color = vbRed
    End With
Next

詳しくは、結果のところで見てみたいと思います。

それから散布図の信頼区間と予測区間ですが、母回帰が回帰直線上の ŷ_i を平均とする正規分布を仮定しているため、どちらも回帰直線をはさんで滑らかな曲線を描きます。

今回、この滑らかさを利用して、観測値から導かれる予測値に基づく２次多項式の近似曲線を引くことで、区間の表示を実現しています。

'信頼区間と予測期間の近似曲線を追加
For i = 2 To 5
    With .FullSeriesCollection(i)
    'マーカー非表示
    .Format.Fill.Visible = msoFalse
    .MarkerForegroundColorIndex = xlColorIndexNone
    .Trendlines.Add
    With .Trendlines(1)
        .Type = xlPolynomial '多項式
        .Order = 2 '２次
    End With
Next

ですので、標本データの独立変数のレンジに比べてサンプルサイズが極端に小さい（観測数が少ない）など、独立変数の観測値にスキマが多い場合は、線が歪むこともあろうかと思います。ご不便をおかけしますが、ご使用の際はこの点ご承知おきください。

*1:今回から、タイトルを「実行プログラム」→「マクロ」に変えました。