線性回歸(10)——多項式回歸

本來以為分享到第9期，線性回歸系列就算完結了，但是最後想了想，線性回歸大家族中的非線性回歸系列在臨床和科研中運用也較為常見，所以最終決定把多項式回歸，分段回歸還有限制性立方樣條也加進來。現實中，自變量和因變量的線性關係其實很少，那麼，如何處理變量之間的非線性關係便成為了一種問題。多項式回歸，分段回歸還有限制性立方樣條則是臨床及科研中較為常見的處理這種情形的方法。1.多項式回歸的理解多項式回歸通過對變量進行升維將變量之間的非線性關係轉化為線性關係。以預測變量的冪形式作為新的預測變量，其實本質上還是屬於線性回歸。這是一種通過增加自變量上的次數，而將數據映射到高維空間的方法，只要我們設定一個自變量上的次數（大於1），就可以相應地獲得數據投影在高次方的空間中的結果。
對於y=wx+b，是個簡單的一元一次方程，我們將其轉化為多項式回歸方程：

不難注意到，多項式變化後數據看起來不太一樣了：首先，數據的特徵（維度）增加了，這正符合我們希望的將數據轉換到高維空間的願望。其次，維度的增加是有一定的規律的。不難發現，如果我們本來的特徵矩陣中只有一個特徵x。假設x=1，2，3，則有：

假設這個被投影到更高空間中的數據在某個角度上看起來已經一條直線了，於是我們可以繼續使用線性回歸來進行擬合，線性回歸是會對每個特徵擬合出權重w的，由此推斷，假設多項式轉化的次數是n，則數據會被轉化成形如：

而擬合出的方程也可以被改寫成：

這個過程看起來非常簡單，只不過是將原始的x上的次方增加，並且為這些次方項都加上權重w，然後增加一列所有次方為0的列作為截距乘的x0。其實，如果把x0~xn改為Z1~Zn，那麼該方程跟普通的線性回歸就完全一樣了：

這是個妥妥的線性方程：並不存在任何高次項，只要擬合結果不錯，大概也沒有人會在意這個幾個特徵的原始數據究竟是怎麼來的，那多項式回歸不就變成一個含有部分共線性的線性方程了麼？在很多資料中，多項式回歸被稱為「線性回歸的一種特殊情況」，這就導致了很多人會產生混淆，多項式回歸究竟是線性模型還是非線性模型呢？這就需要聊到我們對」線性模型「的狹義和廣義定義了。狹義線性模型：自變量上不能有高此項，自變量與標籤之間不能存在非線性關係。

廣義線性模型：只要標籤與模型擬合出的參數之間的關係是線性的，模型就是線性的。這是說，只要生成的一系列之間沒有相乘或者相除的關係，我們就認為模型是線性的。

就多項式回歸本身的性質來說，如果我們考慮狹義線性模型的定義，那它肯定是一種非線性模型沒有錯——否則如何能夠處理非線性數據呢，並且在統計學中我們認為，特徵之間若存在精確相關關係或高度相關關係，線性模型的估計就會被「扭曲「，從而失真或難以估計準確。多項式正是利用線性回歸的這種」扭曲「，為線性模型賦予了處理非線性數據的能力。但如果我們考慮廣義線性模型的定義，多項式回歸就是一種線性模型，畢竟它的係數之間也沒有相乘或者相除。總結一下，多項式回歸通常被認為是非線性模型，但廣義上它是一種特殊的線性模型，它能夠幫助我們處理非線性數據，是線性回歸的一種進化。另外一個需要注意的點是，線性回歸進行多項式變化後被稱為多項式回歸，但這並不代表多項式變化只能夠與線性回歸連用。在現實中，多項式變化瘋狂增加數據維度的同時，也增加了過擬合的可能性，因此多項式變化多與能夠處理過擬合的線性模型如嶺回歸，Lasso等來連用，與在線性回歸上使用的效果是一致的，感興趣的話大家可以自己嘗試一下。2.多項式回歸的R實現

x <- c(4,8,12,25,32,43,58,63,69,79)y <- c(20,33,50,56,42,31,33,46,65,75)plot(x,y)
上圖給我們一個直觀的感受：x和y不是線性關係。我們可以先線性擬合一下，看下效果：
dd <- as.data.frame(x=x,y=y)lm1 <- lm(y~x,data = dd)mean((y - predict(lm1,dd))^2)
library(ggplot2)theme_set(theme_classic())ggplot(dd,aes(x,y))+  geom_point()+  geom_smooth(method = 'lm',formula = y~x,se = F)
均方誤差160.547，而且直線擬合如上圖所示，效果很差，所以我們接下來使用多項式回歸看下效果。多項式的代碼實操（以二次項為例）
lm2 <- lm(y~x+I(x^2),data = dd)summary(lm2)
lm3 <- lm(y~cbind(x,x^2),data = dd)summary(lm3)
lm4 <- lm(y~poly(x,degree = 2,raw = F))summary(lm4)
lm5 <- lm(y~poly(x,degree = 2,raw = T))summary(lm5)

當參數raw設置為False（默認）：
表示該模型使用正交矩陣。正交矩陣代表著各個預測變量互相垂直，內積為0，這表示各個變量之間不用考慮交互項。舉個例子，當 x 具有在區間 (0,1) 上的均勻分布時，x 和 x^2 具有大約0.97的相關性，而線性回歸的前提是兩個預測變量之間要互相獨立，不能存在高度相關或精確相關。
當參數raw設置為TURE：
表示該模型等價於lm2和lm3，直接把x^2看成一個獨立的變量放入線性回歸方程中。使用poly()能夠通過產生正交矩陣來防止變量之間的相關性，所以接下來的分析我們接著使用poly()函數。
mean((y - predict(lm4,dd))^2)
ggplot(dd,aes(x,y))+  geom_point()+  geom_smooth(method = 'lm',formula = y~poly(x,degree = 2,raw = F),se = F)
最佳Degree如何確定從上圖得知，二階應該不是最佳的。在多項式中，階數可以是無窮盡的，但是，我們如何知道幾階是最佳的呢？當然，可以通過臨床經驗或者基礎知識來確定。但是很多情況下，我們沒有辦法通過先驗知識來確定合理階數。所以，如何確定最佳階數成了多項式擬合的一個問題。
mse <- c()for (i in 2:9) {  lm <- lm(y~poly(x,i),dd)  mse[i] <- mean((y - predict(lm,dd))^2)}mse <- data.frame(Degree=2:9,MSE=mse[-1])
ggplot(mse,aes(Degree,MSE))+  geom_point(size=3,color='firebrick3')+  geom_line(color='grey50',size=0.7)+  scale_x_continuous(breaks = 2:9)
我們可以看到，當Degree=4的時候模型的MSE都已經趨於穩定了，為了使模型不過於複雜或發生過擬合，則我們認為4是最佳階數。一般的，採用多項式建模確實會較為明顯的提高擬合優度，如前面舉的例子，當直接把x和y擬合成簡單線性回歸時，MSE=161，而當擬合為9階多項式回歸時，MSE接近於0.泰勒公式告訴我們如果一個函數足夠光滑，那麼就可以在函數上某點的一個領域內用一個多項式來對函數無限逼近，而且隨著多項式階數的提高，這種逼近的效果也會越來越好。同理，如果確實有一條光滑的曲線可以對所有數據點都進行毫無偏差的擬合，理論上就可以找到一個多項式對這條曲線進行較為精確的擬合。
pp <- lapply(1:9,function(i){  ggplot(dd,aes(x,y))+    geom_point(size=3,color='#984ea3')+    geom_smooth(method = 'lm',formula = y~poly(x,i),se = F,color='#fdc086')+    ggtitle(paste0('Degree = ',i))+    theme(plot.title = element_text(hjust = 0.5))+    annotate('text',label=paste0('MSE = ',format(mean((y - predict(lm,dd))^2),digits = 3)),             x = 30,y=max(dd$y))})
cowplot::plot_grid(plotlist = pp,nrow = 3)
可以清楚的看到當3階時，效果還是一般，但從4階後，對於數據的擬合性就明顯優於3次項了，當Degree=8 時，曲線呈現出明顯的震蕩，這說明模型出現了過擬合現象。
降低Degree (如上推薦，D>4可能會導致over-fitting)；
增加訓練樣本；
添加正則化項
Fitting Polynomial Regression in R | Datascience+
【機器學習】多項式回歸
菜菜的Sklearn
多項式回歸- 維基百科