折れ線型回帰分析。非線形回帰としてのスプライン１

線形回帰は汎用性が高く、モデルの解釈もしやすいとあってよく使われます。
しかしデータの間に直線関係=比例関係があると考えるのは最も単純な仮定です。
以前の記事では一次式でなく二次式以上の多項式を使うことによって直線ではない回帰分析を試みました。
しかし直線は維持したまま非線形回帰をする方法もあります。
たとえば折れ線を使った回帰です。

【折れ線回帰モデル定式化】

線形回帰のモデル式は次の式で表すことができました。

$\displaystyle y=\beta_0+\beta_1x$

$x$ が説明変数で $y$ が目的変数、 $\beta_0, \beta_1$ がパラメータです。
しかし分析したいデータはある変化点において前提条件が変わってしまったとしましょう。
その場合、その変化点の前後では異なる直線を使いたいと考えるでしょう。
そのためには変化点を $c$ とすればモデル式は

$y=\begin{cases}\beta_{10}+\beta_{11}x & (x\leq c)\\ \beta_{20}+\beta_{21}x & (x > c)\end{cases}$

となります。

このモデル式はステップ関数 $I(x)$

$I(x)=\begin{cases}0 & (x\leq 0)\\ 1 & (x>0)\end{cases}$

を使って

$y = (\beta_{01}+\beta_{11}x)I(c-x) + (\beta_{02}+\beta_{12}x)I(x-c)$

と書き直せます。

さらに任意の $x$ に対して成立するステップ関数についての基本的な関係式

$I(c-x)+I(x-c)=1$

を代入すると

$y = \beta_0+\beta_1x+\beta_2xI(x-c)+\beta_3I(x-c)\ \cdots\ (1)$

とも書けます。
ただし

$\beta_0=\beta_{01}\\\beta_1=\beta_{11}\\\beta_2=\beta_{12}-\beta_{11}\\\beta_3=\beta_{02}-\beta_{01}$

です。

このモデル式は自由度４（パラメータが４つ）であり、これらのパラメータをいつもの線形回帰と同じように最小二乗法によって推定することができます。
これを使って回帰分析すると次のような結果が得られます。

線が切れてしまっていますのでこれでは折れ線とは言えません。
自由度を１つ犠牲にして線を繋ぐ制約条件を加えましょう。

$\beta_{01}+\beta_{11}c=\beta_{02}+\beta_{12}c$

この制約条件を加えると、簡単な計算によりモデル式(1)は

$y = \beta_0+\beta_1x+\beta_2(x-c)I(x-c)\ \cdots\ (2)$

と整理することができます。
これが今回我々が望む折れ線回帰のモデル式です。

同じデータに当てはめると次のような図が得られます。

折れ目をもう少し前に持ってきたらもっと当てはまりがよくなるかな、という感じです。

なお、2ヶ所 $c_1,c_2$ で折れ目を付けたい場合は

$y = \beta_0+\beta_1x+\beta_2(x-c_1)I(x-c_1)+\beta_3(x-c_2)I(x-c_2)$

といくらでも拡張していくことができます。
折れ目を１つ増やすたびにモデルの自由度が１増えていくことになります。

【分析例】

先の例で出したグラフはe-Statから取ってきた家賃の１９７０以降、過去４０年以上の月次推移です。一応ここにも置いておきます。

bukka

個人的にはバブル崩壊でも家賃相場には全く影響せずに上がり続けていたのは意外でした。
バブルというのは本当に株と地価だけの話なんですねえ。投機に参加せず平凡に暮らしていれば何の影響も無かったのでしょう。
なお現在は黒田バズーカもなんのそのと人口減少に伴い単調現象中です。

#R
library(dplyr)
library(lubridate)
b0 <- read.csv("../data/bukka.csv", stringsAsFactor = FALSE) %>%
        rename(ym = Group.Item) %>%
        select(ym, Rent) %>%
        mutate(year = ym %/% 100, month = ym %% 100) %>%
        mutate(date = make_date(year, month, 1)) %>%
        select(-ym)
sp1 <- b0 %>% 
        group_by(year) %>%
        summarise(Rent = mean(Rent)) %>%
        mutate(year.p = if_else(year < 2001, 0, year - 2001)) %>%
        lm(Rent ~ ., data =.)
summary(sp1)

#Call:
#lm(formula = Rent ~ ., data = .)

#Residuals:
#    Min      1Q  Median      3Q     Max
#-7.8715 -3.3030  0.9138  3.6796  5.6153

#Coefficients:
#              Estimate Std. Error t value Pr(>|t|)
#(Intercept) -4.627e+03  1.421e+02  -32.56   <2e-16 ***
#year         2.367e+00  7.151e-02   33.10   <2e-16 ***
#year.p      -3.169e+00  2.060e-01  -15.39   <2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Residual standard error: 4.238 on 45 degrees of freedom
#Multiple R-squared:  0.9687,    Adjusted R-squared:  0.9673
#F-statistic: 696.5 on 2 and 45 DF,  p-value: < 2.2e-16

それぞれの説明変数のP値もバッチリです。
それは上の図で見て一本の直線よりも当てはまりが良いことからも明らかです。

【エクセルでの実現方法】

モデル式(2)はエクセルの数式でも簡単に実装できます。
説明変数がB列に入っていたとすれば、C列に次の式を追加するだけです。

=(B2-2000)*IF(B2-2000>0, B2 - 2000, 0)

ここで2000とは折れ目の値であり、モデル式(2)でいうところの $c$ のことです。
あとはこれを回帰分析の説明変数にC列も加えてやればOKです。

参考文献

An Introduction to Statistical Learning with Applications in R

An Introduction to Statistical Learning: with Applications in R (Springer Texts in Statistics)

posted with カエレバ

G., Witten, D., Hastie, T., Tibshirani, R. James Springer 2013-08-12

Amazonで購入

楽天市場で購入