7.3 回帰モデルの審査

観測値\(y\)と対応する当てはめ値\(\hat{y}\)との差は、訓練セット誤差、あるいは「残差」で、以下のように定義されます。 \[\begin{align*} e_t &= y_t - \hat{y}_t \\ &= y_t - \hat\beta_{0} - \hat\beta_{1} x_{1,t} - \hat\beta_{2} x_{2,t} - \cdots - \hat\beta_{k} x_{k,t} \end{align*}\] ただし、\(t=1,\dots,T\)の各時点です。各残差は対応する観測値のうち予測不能な成分です。

残差には、以下の2つを含む、いくつか有益な性質があります。 \[ \sum_{t=1}^{T}{e_t}=0 \quad\text{かつ}\quad \sum_{t=1}^{T}{x_{k,t}e_t}=0\qquad\text{全ての$k$で} \] これらの性質の結果として、残差の平均はゼロで、残差と予測変数の観測値の相関係数もゼロであることは明白です。(モデルから切片を省いた場合はこの限りではありません。)

回帰変数を選んで回帰モデルを当てはめた後は、モデルの想定が満たされているかチェックするため、残差をプロットすることが必要です。当てはめたモデルと根底にある想定の、さまざまな面をチェックするために生成すべき一連のプロットがあります。では、一つ一つ議論していきましょう。

残差のACFプロット

時系列データでは、ある変数の現時点での観測値が、その1つ前時点の観測値と、あるいはさらに以前の観測値とさえも、似ていることは、よくあることです。ですから、時系列データに回帰モデルを当てはめると、残差に自己相関が見つかるのは普通にあることです。この場合、推計されたモデルは誤差の無相関の想定を破ることになり、予測は非効率ということかもしれません。つまり、より良い予測を得るためにモデルが説明すべき情報が残っている、ということです。それでも、誤差に自己相関があるモデルからの予測にバイアスはないので、「間違っている」わけではありません。ただ、通常、必要以上に区間予測は大きくなってしまいます。ですから、いつも残差のACFプロットを見るべきです。

残差のヒストグラム

残差が正規分布かチェックすることは、いつも良い考えです。先に説明したように、予測に欠かせないわけではありませんが、正規分布だと区間予測の計算がずっと容易になります。

事例

5.3節で紹介したgg_tsresiduals()関数を使って、上記の全ての役立つ残差診断を得ることができます。

fit_consMR %>% gg_tsresiduals()
米国四半期個人消費の回帰モデルの残差分析

図 7.8: 米国四半期個人消費の回帰モデルの残差分析

augment(fit_consMR) %>%
  features(.innov, ljung_box, lag = 10, dof = 5)
#> # A tibble: 1 × 3
#>   .model lb_stat lb_pvalue
#>   <chr>    <dbl>     <dbl>
#> 1 tslm      18.9   0.00204

時間プロットは変動が時を経るにつれいくらか変化していることを示していますが、それ以外に目立ったものはありません。この不均一分散のせいで、被覆確率ごとの区間予測が不正確になるおそれがあります。

ヒストグラムは残差の分布がやや歪んでいることを示しており、これもまた、区間予測の被覆確率に影響しかねません。

自己相関プロットはラグ7で顕著に飛び出ており、Ljung-Box検定では5%水準で有意です。しかし、自己相関は特に大きいわけではなく、ラグ7では、予測や区間予測にさして影響があることはなさそうです。10章で、残差に残った情報をより上手く取り込むために使う動学的回帰モデルを議論します。

残差と予測変数の散布図

残差は体系的なパターンがなく、ランダムに散らばっていて欲しいと思っています。そのことを簡単に素早くチェックする方法は、残差と各予測変数の散布図を吟味することです。これらの散布図にパターンがあるなら、関係は非線形かもしれず、モデルもそのように修正する必要が出てきます。非線形回帰の議論は7.7節を参照してください。

残差はモデルに含めなかった予測変数に対しても散布図を描く必要があります。もしそれらにパターンがあるのなら、その予測変数を(おそらく非線形な形で)モデルに加える必要があるかもしれません。

事例

7.9は、米国個人消費を予測する重回帰モデルの残差と各予測変数の散布図です。残差はランダムに散らばっているようです。ですから、このケースでは、これで満足です。

us_change %>%
  left_join(residuals(fit_consMR), by = "Quarter") %>%
  pivot_longer(Income:Unemployment,
               names_to = "regressor", values_to = "x") %>%
  mutate(regressor = recode(regressor,
                            "Income" = "可処分所得", "Production" = "鉱工業生産",
                            "Savings" = "個人貯蓄", "Unemployment" = "失業率")) %>% 
  ggplot(aes(x = x, y = .resid)) +
  geom_point() +
  facet_wrap(. ~ regressor, scales = "free_x") +
  labs(y = "残差", x = "")
残差と各予測変数の散布図

図 7.9: 残差と各予測変数の散布図

残差と当てはめ値の散布図

残差と当てはめ値の散布図にも、パターンがあってはいけません。もしパターンがあるなら、残差の分散が一定でないことを意味する、誤差の「不均一分散」があるのかもしれません。この問題が生じたら、予測対象変数を対数や平方根にするなどの変換が必要かもしれません。(3.1節、参照).

事例

先の事例を続けます。図7.10は、残差と当てはめ値の散布図です。ランダムに散らばっているので、誤差は均一分散のようです。

augment(fit_consMR) %>%
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_point() + labs(x = "当てはめ値", y = "残差")
残差と当てはめ値の散布図

図 7.10: 残差と当てはめ値の散布図

外れ値と影響の大きい観測値

データの大半と比べて極端な値を取る観測値を外れ値と言います。回帰モデル係数の推計値に大きな影響を与える観測値を影響の大きい観測値と言います。通常、影響の大きい観測値は\(x\)方向に極端な値を取る外れ値でもあります。

外れ値と影響の大きい観測値を探知する公式な手法がありますが、本教科書の範囲を超えます。2章の始めにお薦めしたように、分析を行う前にデータに慣れ親しむことは、とても重要です。\(y\)と各\(x\)の散布図は、回帰分析のいつも役立つスターティング・ポイントであり、しばしば異常値を特定する助けになります。

外れ値の源の一つは、データ入力誤りです。データの簡単な記述的統計で、意味をなさない最小値、最大値を特定できます。そうした観測値を特定し、それが誤った記録なら、直ちに修正するか、標本から除去すべきです。

外れ値は、単にいくつかの観測値が異なっている場合でも生じます。このケースでは、それらの観測値を除去することは賢明ではないかもしれません。ある観測値を外れ値らしいと特定したら、それを調べ、その背後のあり得る理由を分析することが大切です。観測値を除去するか、保持するかの決定は、(特に、外れ値が影響の大きい観測値である場合)困難かもしれません。そうした観測値は除去した場合と、除去しなかった場合の両方の結果を報告するのが賢明です。

事例

7.11は、(7.1節で紹介した事例で)米国個人消費を可処分所得で回帰する際の、1つの外れ値の影響を強調したものです。左のパネルでは、個人消費の伸び率(%)が誤って-4%と記録されたため、外れ値は\(y\)方向だけに極端な値になっています。オレンジ色の線は、外れ値を含むデータに当てはめた回帰線で、黒色の線は外れ値を除去したデータに当てはめた線です。右のパネルでは、可処分所得が6%増で個人消費は4%減と、外れ値は\(x\)方向にも極端な値になっています。このケースでは、オレンジ色の線は黒色の線と大きく違っており、外れ値は極端に影響が大きくなっています。

外れ値と影響の大きい観測値の回帰への影響

図 7.11: 外れ値と影響の大きい観測値の回帰への影響

見せかけの回帰

時系列は「非定常」なこと、つまり、時系列の値が一定の平均付近を一定の分散で変動するわけではないこと、の方が多いです。時系列の定常性については、9章でより詳細に扱いますが、ここでは、非定常なデータが回帰モデルにもたらし得る影響について述べておく必要があるでしょう。

例えば、図7.12にプロットした2つの変数を考えてみましょう。両者とも同様に上昇トレンドがあるというだけで、関係があるように見えます。しかし、オーストラリアの航空旅客数はギニアの米生産量と何の関係もありません。

オーストラリアの航空旅客数をギニアの米生産量で回帰したこの例が示すように、トレンドのある時系列データ同士は関係があるように見える

図 7.12: オーストラリアの航空旅客数をギニアの米生産量で回帰したこの例が示すように、トレンドのある時系列データ同士は関係があるように見える

非定常時系列を回帰すると見せかけの回帰につながりかねません。図7.13は、オーストラリアの航空旅客数をギニアの米生産量で回帰した出力です。高い\(R^2\)と残差の強い自己相関は、見せかけの回帰の兆候かもしれません。以下の出力で、これらの特徴量を見てください。非定常データや見せかけの回帰を取り巻く問題については、10章でより詳細に議論します。

見せかけの回帰のケースでは、短期予測はもっともらしく見えるかもしれませんが、将来になるほど機能しなくなるのが一般的です。

fit <- aus_airpassengers %>%
  filter(Year <= 2011) %>%
  left_join(guinea_rice, by = "Year") %>%
  model(TSLM(Passengers ~ Production))
report(fit)
#> Series: Passengers 
#> Model: TSLM 
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -5.945 -1.892 -0.327  1.862 10.421 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)    -7.49       1.20   -6.23  2.3e-07 ***
#> Production     40.29       1.34   30.13  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.24 on 40 degrees of freedom
#> Multiple R-squared: 0.958,   Adjusted R-squared: 0.957
#> F-statistic:  908 on 1 and 40 DF, p-value: <2e-16
fit %>% gg_tsresiduals()
見せかけの回帰からの残差

図 7.13: 見せかけの回帰からの残差