5.8 非线性回归

在此之前,我们都假设预测变量和被预测变量之间为线性关系,但是在某些情况下,非线性函数形式可能更加符合真实情况。为了简化本节内容,我们假设只有一个预测变量 \(x\)

一般情况下,可以通过变换预测变量 \(x\) 或被预测变量 \(y\) 来构造非线性回归模型。虽然此时为非线性函数形式,但是模型中的参数认为线性的。最常见的变换方式是对变量做自然对数变换(详见3.2节)。

双对数模型形式可写为: \[ \log y=\beta_0+\beta_1 \log x +\varepsilon. \] 在此模型中,斜率 \(\beta_1\) 可以理解为弹性:当 \(x\) 变动 \(1\%\) 时, \(y\) 的平均变动程度。此外还有其他的对数模型形式: 对数-线性 模型,此类模型只对被预测变量做对数变换处理; 线性-对数 模型,此类模型只对预测变量进行对数变换处理。

需要注意的是,进行对数变换的前提是,所有的观测值都必须大于0。当观测值中存在0值时,我们人为的对所有的观测值加1,然后进行对数变换,即取 \(\log(x+1)\) 。这种变换方式具备一个良好的特点:0的原始比例的整齐效果,即取值为0的观测值在对数变换之后仍为0。

有些情况下,简单的将数据做对数处理是不够的,需要对数据做进一步处理。模型如下所示: \[ y=f(x) +\varepsilon \] 其中, \(f\) 是非线性函数。在标准线性模型中, \(f(x)=\beta_{0}+\beta_{1}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\\ (x-c) & x \ge c \end{array}\right. \end{align*}\] 上述表达式表明,当 \(x-c\) 大于0时取其真实值,当 \(x-c\) 小于0时取值为0。这就使得斜率在 \(c\) 点发生变化,进而展现出非线性形式。

\(x=t\) ,以下通过分段拟合的方式来拟合该时间序列。

以这种方式构造的分段线性关系是回归样条的特例。通常情况下,回归样条如下所示: \[ 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\))的数量及其放置位置,某些软件中可以自动选取数量及放置位置,但由于一些其他原因,这些算法并未被广泛使用。

使用三次回归样条可以得到相较于分段线性回归更为平滑的回归结果。三次回归样条中的约束是连续且平滑的(没有分段线性样条回归中,斜率突然改变的情况)。三次回归样条通常可写为: \[ x_{1}= x \quad x_{2}=x^2 \quad x_3=x^3 \quad x_4 = (x-c_{1})_+ \quad\dots\quad x_{k} = (x-c_{k-3})_+. \] 三次回归样条通常可以更好的拟合数据。但是,当 \(x\) 超出历史数据范围时, \(y\) 的预测值变得不可靠。

非线性趋势预测

5.4节中,我们介绍了通过令 \(x=t\) 拟合存在线性趋势的时间序列数据。拟合非线性趋势的最简单方法是采用二次或更高阶趋势: \[ x_{1,t} =t,\quad x_{2,t}=t^2,\quad \ldots. \] 但是,在预测问题中,推断二阶甚至更高阶趋势是不现实的,因此不建议采用二阶甚至更高阶趋势。

采用上述介绍的分段回归,并拟合出分段的线性趋势是较好的方法。我们可将非线性趋势视为线性趋势的组合。假如在时间\(\tau\)处趋势发生变化,则在模型中可以简单替换上述的 \(x=t\)\(c=\tau\)\[\begin{align*} x_{1,t} & = t \\ x_{2,t} = (t-\tau)_+ &= \left\{ \begin{array}{ll} 0 & t < \tau\\ (t-\tau) & t \ge \tau \end{array}\right. \end{align*}\] 假设 \(x_{1,t}\)\(x_{2,t}\) 的系数分别为 \(\beta_1\)\(\beta_2\) ,那么在时间点 \(\tau\) 之前,时间序列的趋势为 \(\beta_1\) ,而在 \(\tau\) 之后,时间序列的趋势则变为 \(\beta_1+\beta_2\) 。通过增加节点个数,可使得时间序列趋势变化更加平滑。

示例:波士顿马拉松的获胜时间

5.20中的上图为自1897年以来波士顿马拉松赛的获胜时间(以分钟为单位)。从图中可以看出,由于获胜者的马拉松耗时一直在减少,因此时间序列呈现总体下降的趋势。但下图为将数据拟合为线性模型之后的残差图,该残差图表现出明显的非线性形态。并且存在异方差性,即随着时间推移,残差的波动越来越小。

利用线性模型拟合波士顿马拉松获胜者所用时间。

图 5.20: 利用线性模型拟合波士顿马拉松获胜者所用时间。

通过对变量 \(y\) 做对数变换,可以实现对数据的指数趋势拟合: \[ \log y_t=\beta_0+\beta_1 t +\varepsilon_t. \] 这同时也解决了异方差的问题。指数趋势拟合和预测结果如图5.21所示。虽然指数趋势的拟合效果不如线性趋势的拟合效果,但它给出了更合理的预测,即获胜者所用时间以某一衰减速率减少而不是固定线性速率减少。

获胜者所用时间图显示了三个截然不同的时期。在1940年左右,获胜者所用时间出现了大幅度的波动,其中在1920左右,获胜者所用时间大幅降低之后又大幅上升。在1940年之后,获胜者所用时间几乎呈直线下降趋势,在1980年左右,获胜者所用时间趋于平稳。为了解释时间衰减速率的变化,我们可以在1940年和1980年分别加入节点。需要注意的是,若节点过多,可能会导致模型过度拟合,反而不利于模型的预测。

h <- 10
fit.lin <- tslm(marathon ~ trend)
fcasts.lin <- forecast(fit.lin, h = h)
fit.exp <- tslm(marathon ~ trend, lambda = 0)
fcasts.exp <- forecast(fit.exp, h = h)

t <- time(marathon)
t.break1 <- 1940
t.break2 <- 1980
tb1 <- ts(pmax(0, t - t.break1), start = 1897)
tb2 <- ts(pmax(0, t - t.break2), start = 1897)

fit.pw <- tslm(marathon ~ t + tb1 + tb2)
t.new <- t[length(t)] + seq(h)
tb1.new <- tb1[length(tb1)] + seq(h)
tb2.new <- tb2[length(tb2)] + seq(h)

newdata <- cbind(t = t.new, tb1 = tb1.new, tb2 = tb2.new) %>%
  as.data.frame()
fcasts.pw <- forecast(fit.pw, newdata = newdata)

fit.spline <- tslm(marathon ~ t + I(t^2) + I(t^3) +
  I(tb1^3) + I(tb2^3))
fcasts.spline <- forecast(fit.spline, newdata = newdata)

autoplot(marathon) +
  autolayer(fitted(fit.lin), series = "线性") +
  autolayer(fitted(fit.exp), series = "指数型") +
  autolayer(fitted(fit.pw), series = "分段回归法") +
  autolayer(fitted(fit.spline), series = "三次样条法") +
  autolayer(fcasts.pw, series = "分段回归法") +
  autolayer(fcasts.lin$mean, series = "线性") +
  autolayer(fcasts.exp$mean, series = "指数型") +
  autolayer(fcasts.spline$mean, series = "三次样条法") +
  xlab("年份") + ylab("获胜者所用时间(分钟)") +
  ggtitle("波士顿马拉松") +
  guides(colour = guide_legend(title = " "))+
  theme(text = element_text(family = "STHeiti"))+
  theme(plot.title = element_text(hjust = 0.5))
分别利用线性方法、指数方法、分段线性趋势方法和三次样条方法对波士顿马拉松获胜者所用时间进行预测。

图 5.21: 分别利用线性方法、指数方法、分段线性趋势方法和三次样条方法对波士顿马拉松获胜者所用时间进行预测。

5.21分别利用线性方法、指数方法、分段线性趋势方法和三次样条方法对波士顿马拉松获胜者所用时间进行拟合及预测。从图中可以看出,分段线性趋势的预测结果最好,而三次样条回归对数据的拟合效果非常好,但是预测结果较差。

自然三次平滑样条可以很好的解决这个问题,该方法增加了一些额外的约束,导致样条函数在末尾是线性的,该方法既可以很好的拟合数据,又可以产生较好的预测效果。在图5.22中,我们采用 splinef() 函数生成三次样条方法的预测结果。该函数采用了比图5.21中更多的节点,但是通过约束系数来防止模型过拟合,同时尾部也是平滑的。该方法的优点之一是节点选择是非主观的。此外,我们利用对数变换(lambda=0)来解决异方差性。

marathon %>%
  splinef(lambda=0) %>%
  autoplot()+
  xlab('年份')+
  ggtitle("")+
  theme(text = element_text(family = "STHeiti"))+
  theme(plot.title = element_text(hjust = 0.5))
采用自然三次平滑样条方法对马拉松数据进行拟合。预测值是观测数据末尾的线性预测。

图 5.22: 采用自然三次平滑样条方法对马拉松数据进行拟合。预测值是观测数据末尾的线性预测。

5.23为残差图,虽然该图仍表现出一定程度的异方差性,但是可以看出,模型已经较大程度的捕捉到了数据的趋势。预测区间的宽度较大反映了历史数据的波动性。

marathon %>%
  splinef(lambda=0) %>%
  checkresiduals()
三次样条趋势回归结果的残差。

图 5.23: 三次样条趋势回归结果的残差。