SF7 Beyond Linearity

雖然感覺光研究線性模型就能做出很多變化，不過事實是它永遠不夠（科學家也不會甘於用線性模型解釋所有世俗現象）

The truth is never linearity, yet the linearity assumption is good enough (sometimes)

這裡稍微整理一下非線性模型的做法及評估方式：

• Basis Functions

  • Special case
    • Polynomial Regression
    • Step Functions: cut the variable into distinct regions

• Regression Splines

  • keyword: knots, $df$
  • In R: lm(wage ~ bs(age, knots = c(25, 40, 60)), data = Wage)

• Smoothing Splines

  • keyword: $\lambda$, $df_\lambda$ (effective degree of freedom)
  • In R: fit = smooth.spline(age, wage, df = 16)

• Local Regression

  • In R: loess()

• Generalized Addictive Models (GAM)

  • In R: lm(y ~ f_1(x_1)+ f_2(x_2)+f_3(x_3)) ; plot.gam()
  • In R: MGCB package

Polynomial Regression

我們將一般的回歸式進行延伸：

$y_i = \beta_0+ \beta_1x_i + \beta_2x_i^2+ \beta_3x_i^3+ ... + \beta_dx_i^d + \epsilon_i$

一般來說，我們不會讓 $d$ 大於 3 或 4，因為太大的指數會使得模型過於彈性，曲線會長得非常奇怪。

根據 $x_i$ 我們已經得到實線適配的結果，接著我們想要得到該線的信賴區間（虛線）。其實利用最小平方法計算時，就可以得到係數 $\hat{\beta_j}$ 的變異量了，而估計的信賴區間可藉由 pointwise standard error ($Var[\hat{f}(x_0)]$) 就能得到。實際上，兩倍標準誤就是信賴區間了。

如果我們再敏感一點，其實會看到原始資料在 Wage 250 的位置有明顯的分群，假如我們想看各個年齡層在收入大於250的機率，可以用 logistic regression model:

$Pr(y_i > 250|x_i)=\frac{exp(\beta_0+ \beta_1x_i + \beta_2x_i^2+ \beta_3x_i^3)}{1+exp(\beta_0+ \beta_1x_i + \beta_2x_i^2+ \beta_3x_i^3)}$

在 R 的實作上：

fit = glm(I(wage>250) ~ poly(age,3), data=Wage, family="binomial")
preds = predict(fit,list(age=age.grid),se=T)
preds$se.fit
se.bands <- preds$fit + cbind(fit=0,lower=-2*preds$se.fit,upper=2*preds$se.fit)
prob.bands <- exp(se.bands)/(1+exp(se.bands))
matplot(age.grid,prob.bands,col="blue",lwd=c(2,1,1),lty=c(1,2,2),type='l',ylim=c(0,1))

就能繪出類似這樣的結果

Step function

有時也被稱作 piece-wise constant function。概念其實和 Polynomial Regression 類似，只是將資料進行人為的分類，也就是把連續變項轉成次序類別變項，用新的變數 (dummy variable) 來進行模型適配。

$I(\cdot)$ 稱為 indicator function, 當條件成立為 1 不成立 0，其實就是 dummy variable。 step function 常用的領域為生物統計及流行病學。這裡特別注意，除非這些斷點本身有意義，否則不建議用 step function 對原始資料進行切割，因為這樣可能會有資訊遺漏的風險。

Regression Splines

可以看成是 polynomial regression 和 step function 結合的進階版：

$y_i = \begin{cases} \begin{align} \beta_{01}+ \beta_{11}x_i + \beta_{21}x_i^2+ \beta_{31}x_i^3 +\epsilon & \text{ if } x_i <c \newline \beta_{02}+ \beta_{12}x_i + \beta_{22}x_i^2+ \beta_{32}x_i^3 +\epsilon & \text{ if } x_i \geq c \end{align} \end{cases}$

從上圖可看到，原本的 polynomial regression 因為後面的條件 (knot)，被分拆成兩個新的 regression models，因此我們又稱之為 piecewise polynomials。不過光這樣切，其實會有個問題，圖形會變得不連續：

好的做法是讓他在 knots 的位置能夠連續且平滑，我們必須做一些條件的限制。以此例，是 degree = 3 的函數，為了能夠使曲線連續且平滑，我們必須讓 age = 50 在一次及二次微分後的函數的值是一樣的，也就是說我們的 piecewise polynomials 會變成：

$y_i = \begin{cases} \begin{align} \color{royalblue}{\beta_{0}}+ \color{royalblue}{\beta_{1}}x_i + \color{royalblue}{\beta_{2}}x_i^2+ \beta_{31}x_i^3 +\epsilon & \text{ if } x_i <c \newline \color{royalblue}{\beta_{0}}+ \color{royalblue}{\beta_{1}}x_i + \color{royalblue}{\beta_{2}}x_i^2+ \beta_{32}x_i^3 +\epsilon & \text{ if } x_i \geq c \end{align} \end{cases}$

自由度會從原本的 $df=8$ 變成 $df=5$，繪圖的結果會變成連續平滑了：

不過到底要怎麼決定節點數及節點的位置呢？

1. 最簡單的方式就是讓電腦算給你（不負責任的做法…）
2. 直接肉眼觀察曲線的長相
3. 進行交叉驗證：90% 資料做訓練，去比對不同的節點位置或節點數的模型平均的 test RSS 大小，選擇最小 RSS 模型。

通常 regression splines 的結果會優於 polynomial regression，由於兩者的自由度都一樣，但是 splines 多了節點使得曲線更彈性。

Smoothing Splines

除了上述平滑曲線的方式，還有另外一種做法：在有限制的情況下使 $RSS$ 越小越好。如果我們只用原始的 $RSS$ 概念，$RSS \text{ (loss function)} = \sum_{i=1}^{n} (y_i - g(x_i))^{2}$ ，因為要符合最小的結果，一定能透過內差法找出最適合的 $g(x_i)$ 使 $RSS$ 為零，這會導致過度適配。 其中一種方式就是，加入一個非負整數的調節參數 $\lambda$：

$RSS=\sum_{i=1}^{n} (y_i - g(x_i))^{2}+\lambda\int g''(t)^2dt$

$\lambda\int g''(t)^2dt$ 作為懲罰項（限制 $g()$ 的變動大小），二次微分的用意就是看曲線的凹凸性，也就是斜率的變化，我們假想兩種情況，如下圖所示：

以圖 Ａ 來說，函數 $g(t)$ 附近的斜率相當平滑，斜率變化趨近於常數，反觀圖 B，函數 $g(t)$ 附近的斜率變化大，會使得 $g''(t)$ 變得很大。因此 $\lambda$ 的目的就是就是在限制 $g''(t)$ 的收縮，如果 $\lambda \rightarrow \infty$ ，也就是把吸收了大量的懲罰效果，最後的結果就是使曲線變成直線。

Choosing the parameter $\lambda$

有兩種做法可以決定 $\lambda$

• tune effective degrees of freedom ($df_{\lambda}$) rather than $\lambda$

這裡不做證明，但可以得到結論是 $\lambda \text{ increases from 0 to } \infty \text{ as }df_{\lambda} \text{ decreases from n to 2}$ $n \text{: the number of unique }x_i$

example R:

smooth.spline(x, y, df=16) # max df: length(unique(x), min df: 2

• LOOCV:

  $RSS_{cv}(\lambda)=\sum_{i=1}^{n}(y_i-\hat{g}_{\lambda}^{-i}(x_i))^2$,

  $\hat{g}_{\lambda}^{(-i)} \text{ indicates smoothing spline function }$

example R:

smooth.spline(x, y, cv=T)

下圖我們可以看到兩種方式的結果：

Local Regression

這非常有意思，我們先看圖：

可以把它看成是單點估計的回歸，作法如下：

1. $k$ 為分群數，$n$ 樣本數，決定訓練集的比例 $s = k/n$ ，找出離目標點 $x_i$ 的 $s$ 比例的觀察值作為訓練集，$s$ 取得越小，適配模型會越區域化、曲線的跳動頻繁，反之，$s$ 取得越大，適配模型會越整體化，曲線越平滑。
2. 對於這 $s$ 的每個點 $x$ 根據離 $x_i$ 的距離重新賦值 $K_{i0}=K(x_i,x_0)$ ，離 $x_i$ 最遠的點為 0

3. 利用 2. 適配 weighted leasted squares regression: $\sum_{i=1}^{n}K_{i0}(y_i-\beta_0-\beta_1x_i)^2$

4. 由 $\hat{f}(x_0)=\hat{\beta}_0+\hat{\beta}_1x_0$ 得到點估計 $\hat{x}_0$

這種過程屬於 memory-based procedure，除此之外，考量 step-size 也是蠻重要的，也就是選擇 moving fitting 的過程需要多少個目標點 $x_i$。 另外一點要注意的是，不建議在高維度（$p>3 \text{ or } 4$）的回歸進行臨近點分群回歸，因為隨著維度增加，鄰近目標點 $x_0$ 的 $x$ 會越來越少，適配模型效果會非常的差

Theoretically the same approach can be implemented in higher dimensions, using linear regressions fit to p-dimensional neighborhoods. However, local regression can perform poorly if p is much larger than about 3 or 4 because there will generally be very few training observations close to $x_0$.

Generalized Addictive Models (GAM)

GAM 概念其實蠻直覺的，原始的線性迴歸公式：

$y_i = \beta_0+ \beta_1x_i + \beta_2x_i^2+ \beta_3x_i^3+ ... + \beta_dx_i^d + \epsilon_i$

我們可以直接針對各獨變項給予非線性函數：

$y_i=\beta_0+\sum_{j=1}^{p}f_j(x_{ij})+\epsilon = \beta_0+ f_1(x_{i1})+ f_2(x_{i2})+...+f_p(x_{ip})+\epsilon_i$

但有沒有發現，變項間的交互作用項消失了，當然我們可以自己將 $f_{jk}(X_j,X_k)$ 放入模型，這個而適配交互作用項，可以透過上述的 local regression 或是 splines 做法啦！