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

前回のラストで、『統計学がわかる【回帰分析・因子分析編】』のデータが、線形よりも２次多項式回帰モデルのほうが当てはまりがよい（精度が高い）ことに気づきました。
f:id:cyclo-commuter:20180202100706j:plain

横軸に最高気温、縦軸に架空のアイスクリーム店への来店者数をとって、日ごとの観測値がプロットされているのですが、回帰式から、３４℃を頂点にそれより涼しくなるほど、または暑くなるほど来店者数が減る予測が立っています。暑くなるほどアイスが食べたい。でも暑すぎると外に出たくない、って感じでしょうか。

そこで今回は、前々回作ったマクロを実行後、分布の形状を見て近似曲線を２次多項式に変えたときに、信頼区間と予測区間を再計算するマクロを作ってみたいと思います。

計算方法
ソースコード
コードの解説
使い方と実行結果

計算方法

とは言っても、２次多項式回帰の区間推定ってどうやるの？

と思ったら、Yahoo知恵袋に答えがありました。

質問した方と回答した方に感謝です！
回答によると、

式は違っても，いずれも残差分散Ｖeや標本数n，ｘの平均，ｘの平均からの偏差二乗和 Sxx を用います。
つまり，信頼区間を求めるのに，式の形がどうであるか，は関係ないのです。
回帰式が対数であっても，その式とデータとの残差の平方和を求めてＶｅを出せば，単回帰と同様に信頼区間が求まります。

との由。
とすると、新しい回帰式に基づいて予測値を再計算して、その後は前と同じ計算を繰り返せばいいだけ、なんですね。得たり賢し。

ソースコード

'近似曲線を２次多項式にした際の再計算
Sub Bivariate_Polynomial()

    Dim Cht, Tmp, Rng, y, x, n, a, b, c, xbar, Sxx, Se, F, haty, CI, PI, i, j
    
    If Selection.Columns.Count <> 2 Then
        MsgBox "独立変数（x）、従属変数（y）の順に２列の" & vbCrLf & _
            "セル範囲を選択してから実行してください。", vbCritical
        Exit Sub
    ElseIf Selection(1).Row = 1 Then
        MsgBox "ラベル行は選択しないでください。", vbCritical
        Exit Sub
    End If
    
    Set Cht = ActiveSheet.ChartObjects(1).Chart
    With Cht.SeriesCollection(1).Trendlines(1)
        With .DataLabel
            .Left = 29
            .Top = 11
        End With
        Tmp = .DataLabel.Text   '回帰式を取得
    End With
    i = InStr(1, Tmp, "x2")
    j = InStr(i + 1, Tmp, "x") - 1
    
    If i = 0 Then
        MsgBox "散布図中の近似曲線を「多項式近似」" & vbCrLf & _
            "次数「２」にしてから実行してください。", vbCritical
        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)
    a = Val(Mid(Tmp, 4, i - 1))                     '回帰係数(x^2の)
    b = Val(Mid(Tmp, i + 3, j))                     '回帰係数(xの)
    c = Val(Mid(Tmp, j + 2, InStr(j, Tmp, "") - 1)) '切片
    With Application.WorksheetFunction
        xbar = .Average(x)                          'xの平均値
        Sxx = .DevSq(x)                             'xの偏差平方和
        For i = 1 To n
            haty(i) = a * x(i, 1) ^ 2 + b * x(i, 1) + c 'yの予測値
        Next
        Se = .SumXMY2(y, haty)                      '残差平方和
        
        'Ｆ検定
        F = ((Syy - Se) / 2) / (Se / (n - 3))
        Rng.Offset(1, 7).Value = .F_Dist_RT(F, 2, n - 3)
        
        F = .F_Inv_RT(0.05, 2, n - 3)   '区間推定用のＦ値
    End With
        
    '信頼区間と予測区間と標準化残差
    For i = 1 To n
        '信頼区間の幅
        CI(i) = Sqr(F * (1 / n + (x(i, 1) - xbar) ^ 2 / Sxx) * Se / (n - 3))
        '予測区間の幅
        PI(i) = Sqr(F * (1 + 1 / n + (x(i, 1) - xbar) ^ 2 / Sxx) * Se / (n - 3))
        With Rng
            .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 - 3)) '標準化残差
                '標準化残差の絶対値が３を超える場合は警告
                If Abs(.Value) > 3 Then .Font.Color = vbRed
            End With
        End With
     Next
    
    With Cht
        '横軸の最小値
        .Axes(xlCategory).MinimumScale = Round(WorksheetFunction.Min(x) - 1, 0)
        '縦軸の最小値
        .Axes(xlValue).MinimumScale = _
            Round(WorksheetFunction.Min(Rng.Offset(1, 2).Resize(n)) - 1, 0)
    End With

    Application.ScreenUpdating = True

End Sub

コードの解説

ほとんど前々回に作ったマクロの丸写しなのですが、再計算なのでセル書式のところとかグラフ作成の部分はカットしてあります。

新たに付け加えたり変更したのは２箇所。

ひとつは２次多項式から回帰係数と切片を取得するところです。

Set Cht = ActiveSheet.ChartObjects(1).Chart
With Cht.SeriesCollection(1).Trendlines(1)
Tmp = .DataLabel.Text '回帰式を取得
End With
i = InStr(1, Tmp, "x2")
j = InStr(i + 1, Tmp, "x") - 1
a = Val(Mid(Tmp, 4, i - 1)) '回帰係数(x^2の)
b = Val(Mid(Tmp, i + 3, j)) '回帰係数(xの)
c = Val(Mid(Tmp, j + 2, InStr(j, Tmp, "") - 1)) '切片

ご覧のとおり、グラフのデータラベルのテキストから一つずつ抜き出しています。

もうひとつは、残差分散・F 検定・F 値の計算が前と違ってます。

Syy = .DevSq(y) 'yの偏差平方和
F = ((Syy - Se) / 2) / (Se / (n - 3))
Rng.Offset(1, 7).Value = .F_Dist_RT(F, 2, n - 3)
F = .F_Inv_RT(0.05, 2, n - 3) '区間推定用のＦ値

前回は F = (a ^ 2 * Sxx) / (Se / (n - 2)) で F 値を求めましたが、回帰係数がaとbの２個になったので、(a ^ 2 * Sxx) は使えません。一瞬途方に暮れましたが、ネットでいろいろ調べるうちに、F 値は［回帰分散÷残差分散］で求められるんだから、重回帰分析に準じて計算すればいいじゃん、ということに気づきました。
そこで『マンガでわかる統計学　回帰分析編』１２４ページの数式

${\frac{S_{yy}-S_{e}}{説明変数の個数}÷\frac{S_{e}}{個体の個数-説明変数の個数-1}}$