【学習動機】
折れ線回帰や区分回帰(piecewise regression, segmented regressionといったところか)をこれまで触れたことが無かったのだが, 少し触れる機会があった. どうやって推定するものかと頭で考えていたのだけれど, よく分からなかったので, 調べてみたら理解できた. ついでに, いつものように人工データを作って, 実験してみる.
基本的には, 以下の記事を参考にした.
【学習内容】
パッケージ読み込み.
library(tidyverse)まずは, 真の構造を作成. 20のときに最大を迎える関数ftrueを作る. 関数の形は以下の図のような感じ. イメージは消化力で, 20歳のときにピークを迎えるように作った.
x <- seq(from = 0, to = 100, by = 0.1) %>% as.numeric()
ftrue <- function(x){
if(x<20){return(5*x)}
else{return(125 - 1.25*x)}
}
df <- data.frame(
x1 = as.numeric(x),
y = apply(array(x),1,ftrue)
)
plot(df$x1,df$y)
{wakefield}パッケージについては, 昔の記事参照.
ftrue関数に, \(N(0,9^2)\)の正規分布の誤差項を加えたものを生成.
後で, 折れ線回帰でなく二乗項を加えて回帰をするために, xxという変数も作成.
set.seed(271)
n <- 100
x <- wakefield::age(n, x=0:100)
df <- data.frame(
x1 = as.numeric(x)) %>% mutate(
xx = x1^2,
y = apply(array(x),1,ftrue) + rnorm(n, mean = 0, sd = 9)
)
plot(df$x1, df$y)
この例で推定するケースを考える. 真のモデルは以下の通り.
\[ y = \begin{cases} 5x + \epsilon & (x < 20) \\ 125 - 1.25x + \epsilon & (x \ge 20) \end{cases} \]
\[ \epsilon ～ N(0,9^2) \] 最初に挙げた記事に照らし合わせると, \(\beta_{10}=0, \beta_{11}=5, \beta_{20}=125, \beta_{21}=-1.25\)ということになる. これらのパラメータを推定することになる.\(c=20\)は既知とする. このように, 折れ線回帰の折れ目が分かっているということは普通なのだろうか. 折れ線回帰を考えながら, マルコフスウィッチングモデルだとか構造変化の検定だとか隠れマルコフモデルだとかを思い出したけど, 今回はこれ以上広げずまた今度にする.
推定するための折れ線回帰(折れ線としてつなぎ目あり)のモデル式は
\[ y = \beta_{0} + \beta_{1}x + \beta_{2}(x-c)I(x-c) \] となる.
このモデルが推定できるように, x2という変数を追加する. 上の\(\beta_{2}\)を係数として持つ変数の部分である.
この推定方法を見ていて, サポートベクターマシーンでカーネルトリックで次元を一つ増やしてやって平面分離するイメージが思い浮かんでいた. 次元を増やしてきれいに解決するという点ですごくスッキリした. まぁそんなことはともかく, 推定してみた.
df2 <- df %>% mutate(x2 = if_else(x1>20,x1-20,0))
res1 <- lm(y~x1+x2, data = df2)
summary(res1)##
## Call:
## lm(formula = y ~ x1 + x2, data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.8269 -5.2227 0.5048 7.2089 20.8555
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.0248 2.9980 -1.009 0.316
## x1 5.1876 0.1943 26.693 <2e-16 ***
## x2 -6.4207 0.2206 -29.112 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.547 on 97 degrees of freedom
## Multiple R-squared: 0.9122, Adjusted R-squared: 0.9104
## F-statistic: 503.7 on 2 and 97 DF, p-value: < 2.2e-16係数の対応関係に合わせて, 元のパラメータに戻すと
\[ \beta_{01} = \beta_{0} \\ \beta_{11} = \beta_{1} \\ \beta_{12} = \beta_{1} + \beta_{2} \\ \beta_{02} = \beta_{0} + \beta_{3} =\beta_{0} - c \beta_{2} \] となるので, \(\beta_{01}=\) -3.0248435, \(\beta_{11}=\) 5.1876236, \(\beta_{12}=\) -1.2330436, \(\beta_{02}=\) 125.3884994 という推定結果であったが, 真のモデルパラメータは順番に0, 5, -1.25, 125なので悪くない結果だ. 真のモデルと同じ構造をあてはめているのだから, まぁそれはそうというところだが.
ただし, 戻したパラメータの不確実性はどうなるのだろうか. 少し気になるところではある.
さて, 折れ線回帰が思いつかない場合, 思いつくモデルは二乗項を加えたモデルだろうから, 推定してみた.
res2 <- lm(y~x1+xx, data = df2)
summary(res2)##
## Call:
## lm(formula = y ~ x1 + xx, data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.548 -10.847 -2.097 10.885 53.518
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.007703 4.689733 5.546 2.53e-07 ***
## x1 2.320621 0.221435 10.480 < 2e-16 ***
## xx -0.027601 0.002166 -12.742 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.22 on 97 degrees of freedom
## Multiple R-squared: 0.6801, Adjusted R-squared: 0.6735
## F-statistic: 103.1 on 2 and 97 DF, p-value: < 2.2e-16AICも比較してみた. 当然ながら, 真のモデルと同じ構造で推定しているres1の方がAICが低い.
AIC(res1,res2)## df AIC
## res1 4 739.9789
## res2 4 869.2265【学習予定】
上で挙げた疑問の解決と, Rのパッケージで{segmented}というものがあるみたいなので使ってみたいところ.
https://cran.r-project.org/web/packages/segmented/segmented.pdf
