7.7 非線形回帰

本章では、ここまで線形関係を想定してきました。それで適切なことが多いのですが、非線形関数式の方が適しているケースも数多くあります。本節では単純化のため、予測変数は$x$1個だけを想定します。

非線形関係をモデル化する最も単純なやり方は、回帰モデル推計の前に、予測対象変数\(y\)、および/または、予測変数\(x\)を変換しておくことです。これで\(y\)\(x\)の関係は非線形になりますが、モデル式は線形のままで、そこにパラメータが付きます。最もよく使われる変換は(自然)対数(3.1節、参照)です。

log-log式は、以下のように指定できます。 \[ \log y=\beta_0+\beta_1 \log x +\varepsilon. \] このモデルでは、傾き\(\beta_1\)は弾力性、つまり、\(x\)の1%の変化が平均すると\(y\)\(\beta_1\)伸び率(%)をもたらす、と解釈できます。他にも、役に立つ式があります。log-linear式では、予測対象変数だけを変換します。linear-log式では、予測変数を変換します。

変数に対数変換を行うためには、全ての観測値がゼロ超でなければならないことを思い起こしましょう。変数\(x\)にゼロがあるケースでは、\(\log(x+1)\)変換を使います。つまり、変数に1を加え、それから対数を取るのです。これで、対数を取るのと似た効果が得られると同時に、ゼロの問題を回避できます。元の尺度でゼロのものは変換後の尺度でもゼロになるので、すっきりするという副作用もあります。

単にデータを変換するだけでは上手く行かず、モデルをより一般的に指定する必要があるケースもあります。そのケースでは、以下のモデル式を使います。 \[ y=f(x) +\varepsilon \] ただし、\(f\)は非線形関数です。標準的な(線形)回帰では、\(f(x)=\beta_{0} + \beta_{1} x\)です。非線形回帰では、単純な対数やその他の変換よりも柔軟な\(x\)の非線形関数を\(f\)に指定できます。

最も単純な指定の一つは、\(f\)区分線形にすることです。つまり、\(f\)の傾きが変化できる点をいくつか導入します。これらの点はノット(結び目)と呼ばれます。\(x_{1,t}=x\)として、以下のような変数\(x_{2,t}\)を導入することで指定できます。 \[\begin{align*} x_{2,t} = (x-c)_+ &= \left\{ \begin{array}{ll} 0 & x < c \text{の場合}\\ x-c & x \ge c \text{の場合} \end{array}\right. \end{align*}\] \((x-c)_+\)の表記は、それが正なら\(x-c\)、それ以外なら0、という意味です。これで、傾きは点\(c\)で変化できます。上記式の変数を追加することで、関係に含める屈折を追加できます。

こうしたやり方で構築される区分線形関係は、回帰スプラインの特殊ケースです。一般的に、線形回帰スプラインは以下のようにして得られます。 \[ x_{1}= x \quad x_{2} = (x-c_{1})_+ \quad\dots\quad x_{k} = (x-c_{k-1})_+ \] ただし、\(c_{1},\dots,c_{k-1}\)はノット(線形が屈折できる点)です。ノット数(\(k-1\))とノットを置く位置の選択は難しく、恣意的になり得ます。ノット選択を自動化するアルゴリズムにいくつか利用可能なものがあるのですが、広く使われてはいません。

非線形トレンドによる予測

7.4節では、\(x=t\)と設定することで時系列に線形トレンドを当てはめることを紹介しました。非線形トレンドを当てはめる最も単純なやり方は、以下のように指定して、2次、あるいはより高次のトレンドを使うことです。 \[ x_{1,t} =t,\quad x_{2,t}=t^2,\quad \dots \] しかし、2次、あるいはより高次のトレンドを予測に使うのは、お薦めしません。外挿するとしばしば非現実的な予測になるからです。

より良いアプローチは、前節で紹介した区分指定を使って、ある時点で屈折する区分線形トレンドを当てはめることです。これは、区分線形から構築された非線形トレンドと考えることができます。もし\(\tau\)時点でトレンドが屈折するなら、前節の予測変数の設定において、モデル内の\(x=t\)\(c=\tau\)を入れ替えるだけです。 \[\begin{align*} x_{1,t} & = t \\ x_{2,t} &= (t-\tau)_+ = \left\{ \begin{array}{ll} 0 & t < \tau \text{の場合}\\ t-\tau & t \ge \tau \text{の場合} \end{array}\right. \end{align*}\] \(x_{1,t}\)\(x_{2,t}\)の係数をそれぞれ\(\beta_1\)\(\beta_2\)とすると、\(\beta_1\)\(\tau\)時点以前のトレンドの傾き、\(\beta_1+\beta_2\)\(\tau\)時点以降のトレンドの傾きとなります。\(\tau\)を線が屈折する時点「ノット」として、\((t-\tau)_+\)式の変数を追加することで、関係に屈折を追加できます。

事例: ボストンマラソン勝利時間

ボストンマラソン男子の勝利時間にトレンドモデルをいくつか当てはめてみましょう。まず、男子のデータを抽出し、勝利時間を数値に変換します。1924年にコースの長さが(24.5マイルから26.2マイルに)変更されたため、勝利時間は大幅に長くなりました。ですから、それ以降のデータだけを考えることにします。

boston_men <- boston_marathon %>%
  filter(Year >= 1924) %>%
  filter(Event == "Men's open division") %>%
  mutate(Minutes = as.numeric(Time)/60)

7.20の上のパネルは、1924年以降の勝利時間です。勝利時間は年を経るごとに改善されてきたので、時系列はおおよそ下降トレンドを示しています。下のパネルは、データに線形トレンドを当てはめた残差です。線形トレンドでは取り込めなかった非線形パターンが、プロットには明らかにあります。

ボストンマラソン勝利時間に線形トレンドを当てはめるのは不適切

図 7.20: ボストンマラソン勝利時間に線形トレンドを当てはめるのは不適切

以下のように\(y\)変数を変換することで、指数トレンドをデータに当てはめることができます。(log-linear回帰と同値です) \[ \log y_t=\beta_0+\beta_1 t +\varepsilon_t \]7.21は、当てはめた線と予測を示しています。指数トレンドは、線形トレンドよりもデータに当てはまっているように見えませんが、おそらくより意味のある予測になっています。将来の勝利時間は減少しますが、そのペースは線形の固定ペースではなく、改善が衰えるペースになるはずですから。

勝利時間のプロットから3つの異なる時期があったことが分かります。1950年頃までは勝利時間のボラティリティが大きく、ほとんど減少していません。1950年以降は明確に減少しています。しかし、1980年代以降は横ばいになり、標本の最後に向けて増加の兆しさえあります。これらの変化を考慮するため、1950年と1980年をノットとして指定します。ここで警告しておきますが、ノットを主観的に特定することは過剰な当てはまりにつながり、モデルの予測性能を阻害しかねません。ですから、注意深く行うべきです。

fit_trends <- boston_men %>%
  model(
    "線形トレンド" = TSLM(Minutes ~ trend()),
    "指数トレンド" = TSLM(log(Minutes) ~ trend()),
    "区分線形トレンド" = TSLM(Minutes ~ trend(knots = c(1950, 1980)))
  )
fc_trends <- fit_trends %>% forecast(h = 10)

boston_men %>%
  autoplot(Minutes) +
  geom_line(data = fitted(fit_trends),
            aes(y = .fitted, colour = .model)) +
  autolayer(fc_trends, alpha = 0.5, level = 95) +
  labs(y = "分", level = "区間予測",
       title = "ボストンマラソン勝利時間")
 線形、指数、区分線形の各トレンドによるボストンマラソン勝利時間の予測

図 7.21: 線形、指数、区分線形の各トレンドによるボストンマラソン勝利時間の予測

7.21は、線形、指数、区分線形の各トレンドによる当てはめた線と予測を示しています。区分線形トレンドによる予測が最良に見えます。