Linear Regression | 線性迴歸模型 | using AirQuality Dataset

Linear Regression 線性迴歸模型是用來預測連續型目標變數與預測變數間的線性關係，並存在許多資料符合常態分佈與線性關係等基本假設。預測變數可以是數值或類別型態，投入模型進行時類別變數都會被轉換成虛擬變數。對於有遺失值的觀測值則會直接忽略，且易受極端值影響，因此在使用模型前必須做好資料前處理。

Linear Regression Introduction

Linear Regression 線性迴歸模型主要是使用最適線性迴歸線來建立依變數Y和一個或多個自變數X間的關係，並可以此線性迴歸方程式，使用給定自變數X的值來預測依變數Y值（適用於目標變數為連續變數時，若目標變數為二元變數，則適用羅吉斯回歸）。（＊線性迴歸方程式的建立則是使用最小平方法Least Sqaures Method來找尋X和Y之間的關係與趨勢）

Linear Regression 線性迴歸方程式：\(Y=\beta_{1}+\beta_{2}*X+\varepsilon\)

Example Problem

使用資料集為R dataset套件中的airquality。此資料集主要搜集自1973年5月到9月，每日空氣品質相關的衡量指標，包括臭氧濃度、太陽輻射、平均風速、最高溫度，以及資料月份與資料日。

分析目標：找出影響空氣中臭氧(Ozone)濃度的重要影響因子與關係。

目標變數（依變數）：連續變數(臭氧濃度，Ozone)
預測變數（自變數）：連續/類別變數(太陽輻射Solar.R、平均風速Wind、最高溫度Temp、資料月份Month、資料日Day)(*此資料集剛好都是連續變數)

為了達成分析目標，主要的資料處理與分析架構會包括以下步驟：

1. 資料載入與檢視
2. 資料探索(Exploratory Data Analysis, EDA)：視覺化分析
   1. 視覺化預測變數與目標變數間的關係：相關係數矩陣、散佈圖scatter plot
   2. 偵測離群值：盒須圖box plot
   3. 檢查預測變數分佈是否為常態分佈(鐘型)而沒有左右偏移：機率密度圖density plot
3. 資料前處理
   1. 離群值處理
   2. 標準化-使變數間關係衡量不受單位規模數值影響
   3. 變數轉換-使變數接近常態分佈，以符合模型假設
   4. 遺失值填補-使用MICE()套件來預測遺失值
4. 建立預測線性模型
   1. Pooling – 綜合多組遺失值填補後的datasets建立回歸模型，並摘要模型配適結果
   2. CompleteData – 挑選其中一組遺失填補後的dataset建立回歸模型，並摘要模型配適結果
5. 線性回歸模型適合度診斷
   1. 資料是否符合線性迴歸模型假設
   2. 使用gvlma套件檢測
   3. 使用shapiro test 來看變數是否成常態分佈
6. 模型驗證
   1. 檢視(比較)模型配適能力常用指標
   2. 80% 訓練 20%驗證
      •  真實值與預測值相關性(Correlation)
      • 正確率(min_max_accuracy)
      • 錯誤率(MAPE ,Mean absolute percentage error)
   3. 加強版模型效果驗證：K-fold Cross Validation
      1. mean squared error
7. 結論

1. 資料載入與檢視

載入資料

檢視資料

2. 資料探索(Exploratory Data Analysis, EDA)

基本敘述統計

使用summary()函數查看基本敘述統計。可以看到此資料集共有153列，以及6個特徵變數的最大最小值、第一與第三分位數以及遺失值數量。

視覺化分析

視覺化預測變數與目標變數間的關係

使用GGally套件中ggcorr()來查看各變數間的相關係數強弱。

由上圖可以觀察到

(1)具有較強正相關性(相關係數>=0.4)的變數關係有：

• 氣溫(Temp)&月份(Month)
• 臭氧濃度(Ozone)與氣溫(Temp)

(2)具有較強負相關性(相關係數<= -0.4)的變數關係有：

• 風速(Wind)與氣溫(Temp)
• 臭氧濃度(Ozone)與風速(Wind)

但僅參考相關係數，並無法知道兩者是否符合線性關係。

散佈圖scatter plot ＆ 機率密度圖

由上圖可以觀察到三種數據型態，分別為：

1. 對角線上各個變數的機率密度分佈（用來檢視變數分佈是否為常態）
   • 除了Temp稍微看似常態，其餘變數Ozone(右偏嚴重)、Solar.R、Wind都不太滿足常態分佈。
2. 左上方三角區塊變數間的相關係數數值大小
   • 可以與前面ggcorr()畫的圖做比對
3. 左下方三角區塊變數間的散佈圖
   • 可以發現此散佈圖在繪製時，已自動使用各變數的最大最小值將數值標準化，所以變數間的散佈圖將不受尺度大小所影響，都剛好在方形的繪圖範圍中。

雖然ggpar可以一次畫出全部，但如果要進一步加上趨勢輔助線，看變數x和y是否為線性關係，可以使用下列語法檢視。

我們將剛剛得知有較強相關性的變數關係畫成散佈圖如下。(Ozone, Temp, Wind, Month)

由上面四張散佈圖，可大略看到線性關係如下：

1. 雖然氣溫和部分月份數值成正向關，但根據常識，這並不是我們希望的正向關。
2. 臭氧濃度會隨著氣溫上升而上升。
3. 臭氧濃度會隨著風速上升而下降。
4. 氣溫會隨著風速上升而下降。

離群值偵測

經過上述變數關係探索後，我們將僅針對聚焦後的變數(Ozone,Temp,Wind)，進行偵測離群值。

我們將使用合鬚圖(box-plot)法來偵測離群值。

 

由盒鬚圖可知，臭氧濃度(Ozone)跟風速(Wind)都有一些離群值。

由上述發現可知，針對聚焦後的變數(Ozone, Temp, Wind)，接下來需要處理：

1. 離群值處理：Ozone(2), Wind(3)。
2. 變數標準化: 使不衡量尺度（單位）間的變數關係不受數值大小所影響。
3. 變數轉換：盡可能讓變數分佈符合常態，以符合線性迴歸基本假設。
4. 空值填補：Ozone(37)。

3. 資料前處理

3-1. 離群值處理

針對聚焦的三個變數(Ozone,Temp,Wind)，本範例會將離群值填補為NA，再之後跟其他空值一併使用MICE方法預測。

偵測離群值(回傳的是離群數值，而非資料列index)。

將離群值填補為NA。

檢視離群值處理後的結果。

檢視離群值處理後變數分佈的變化。

可以看到：

1. Ozone skewness有減少(skewness = 1.21 -> 0.94)，但視覺上還是沒有非常接近常態分佈。
2. Wind的skewness有從0.34減少到0.06，視覺上已接近常態。
3. Temp的部份因為透過box-plot沒有偵測到明顯離群值，故在此沒有變化。

用新的資料偵測離群值，可發現使用box-plot法，已無離群值。

3-2. 標準化

將所有變數標準化(0 to 1)。

3-3. 變數轉換(\(X_{transform} = \log_{10}X\))

透過變數log10轉換，來使變數更近似常態分配(僅針對Ozone)。

檢視變數轉換後Ozone的機率密度分佈圖。

可以發現經log10轉換後的Ozone_log10的skewness有降低，但視覺上分佈還是不是那麼常態。(skewness = 1.21 -> 0.94 ->0.69)

3-4. 空值填補(MICE, Multivariate Imputation by Chained Equations)

當一個特徵變數空值資料占比超過5%就有必要做空值填補之處理。我們建立一個函數用來檢查每個特徵變數遺失資料列佔比。可以發現臭氧Ozone的遺失比例高達25%(原始資料遺失＋離群值NA)。

在此我們選擇使用MICE Package來處理Ozone的遺失值。（更多有關遺失值處理請參考：「遺失值處理方法」）

MICE代表： Multivariate Imputations by Chained Equations (MICE)，具有以下幾個特點：

• 主要可用來產生多元變數(multivariate)的多組空值填補值(multiple imputations)（可透過mice()函式中的參數m，number of mutiple imputations，來指定要產生幾組填補值解，預設為5組）。
• 透過Gibbs sampling法來對多元變數產稱多組空值填補值。
• MICE填補方法是採用Fully Conditional Specification，即每一個不完整的變數都是分別用獨立的模型來預測空值的。
• MICE可處理連續型變數(continuous)、二元變數(binary)、無順序類別型變數(unordered categorical)、有序類別變數(ordered categorical)。
• MICE在預測目標欄位時(target column)，預設會使用其他非目標欄位來作為預測變數(predictors)。如遇預測變數本身也不完整的情況，最新產生的一組填補值會被用來完整預測變數，再進行目標變數的填補預測。
• 我們可以替不同欄位指定各自的單變數填補模型(univariate imputation model)，使用方法範例：mice(data, meth=c(‘sample’,’pmm’,’logreg’,’norm’)))。
• 單變數填補模型可以使用內建的方法或是自訂方法(mice.impute.myfunc)(mice(method = c(“xxx”, “myfunc“))。
• 常見的幾種內建填補模型包括：
  • pmm (Predictive mean matching) : 任何變數型態
  • cart (Classification and regression trees) : 任何變數型態
  • rf (Random forest imputations) : 任何變數型態
  •  logreg (Logistic regression) : 二元類別型態 (binary)

MICE package套件中的主要函式說明：

• md.pattern() : inspect the missing data pattern
• mice(): impute the missing values *m* times
• with(): analyzed completed data sets => Performs a computation of each of imputed datasets in data. with (data, expr = formula, …)。將產生的m組填補值套用到計算公式。
• pool(): combine parameter estimates => Pool the results of the repeated analyses。綜合m組填補數據所產生的模型的估計係數(estimates)、標準差(std.error)、和p-value。
• pool.compare() : pool.compare(fit1, fit0, method = c(“wald”, “likelihood”), data = NULL)。使用方法”wald”或”likelihood”來比較兩種模型公式套用在所有填補數據的綜合結果。
• complete(): export imputed data 儲存和輸出其中一組完整填補數據
• ampute(): generate missing data => Generate simulated incomplete data 產生模擬的不完整數據

首先，因為MICE處理遺失值的假設為，資料為隨機遺失(Missing Completely At Random, MCAR)，所以在使用套件前，我們先使用mice套件中的md.pattern()函數來檢視資料遺失分佈的狀況。我們可以看到有106列資料是完整的，37列資料缺Ozone，5列資料缺Solar.R，2列資料同時缺Ozone&Solar.R。3列資料缺Wind。Ozone總遺失數為39，Solar.R總遺失數為7，Wind總遺失數為3。

我們使用VIM(Visualization and Imputation of Missing Values)套件中的aggr()函數將上述結果用更有效的視覺化圖表來檢視。

• 左圖：各特徵變數的遺失資料比例
• 右圖：遺失狀況分佈pattern。由圖可知，有近7成資料是完整的，單單遺失Ozone的資料列高達近25%。

 

而如何檢視資料的遺失分佈是屬於隨機(Missing Completely At Random, MCAR)，就是使用VIM套件的marginplot()函數。關於marginplot()的一些特性如下：

• 一次只能畫出兩維度(變量A&B)遺失分佈關係。在此，因為 Ozone和Solar.R的遺失比例最高，我們將檢視此二變數遺失分佈是否隨機
• 圖的兩軸最外側盒鬚圖：變量A中，對應變量B完整（藍色）/遺失（紅色）的分佈。
• 圖的兩軸內側散佈圖：單一變量A在的變量B發生遺失值的散佈圖。
• 圖的最左下角為變量A與變量B同時遺失的列數。其上方與右方數字則分別表單變量遺失的列數。
• 圖中央則為兩變數A&B的散佈圖。
• 如果IQR區間重疊，且盒鬚圖分佈相似，則變量A與變量B的遺失值分佈為隨機，符合MCAR假設。

我們從下圖可觀察到：

• 圖左側紅色盒鬚圖為Ozone中Solar.R遺失資料點的分佈，藍色則為其餘Ozone資料點分佈。
• 圖下側紅色盒鬚圖為Solar.R中Ozone遺失資料點的分佈，藍色則為其餘Solar.R資料點分佈。
• 從圖兩側我們可以發現，IQR區間是重疊的，且紅色藍色盒鬚圖分佈相似，因此Var1(Solar.R)與Var2(Ozone)的遺失分佈為隨機的，符合MCAR假設。

 

使用mice()函數進行空值填補，幾個重要參數如下：

• 參數m表示要產生幾組填補數據集。預設為產生5組填補數據集。我們將它設為50組。
• 參數method (or meth)表示填補方法。在本範例中，因為變數皆為數值，我將使用預設的pmm(predictive mean matching)法。其他方法包括：
  • PMM (Predictive Mean Matching) – For any type variables
  • logreg(Logistic Regression) – For binary variables( with 2 levels)
  • polyreg(Polytomous logistic regression) – For unordered categorical/factor variables (>= 2 levels)
  • polr (Proportional odds model) –  For ordered categorical variables ( >= 2 levels)
• 參數maxit表示用於計算遺失值的迭代次數，預設為5。本範例將迭代次數設為50。

使用summary()來看一下填補的結果。可以看到各特徵變數所使用的補值法分別為何，以及填補各目標變數(y)所投入的預測變數(x)。

查看特定特徵值變數填補後的數值可用$imp$var。（列數為該變數有遺失值的列數，行數為各組填補數據集數，本範例為50組）。以Wind為例，遺失的三筆資料，填補後的50組數據集如下：

4. 建立線性模型 linear regression model

最後我們來使用處理好的數據來建立回歸模型。

4-1. Pooling: 使用剛剛填補的50組的資料建模，並摘要綜合50組模型的效果。建模結果顯示預測變數皆為顯著。(p-value < 0.05)

4-2. 使用complete()隨機挑選一組填補後的dataset來建模(預設是第一組，但也可以透過參數調整)。發現模型配適結果也是顯著的(p-value < 0.05)。

評估係數(Coefficients)指標說明：

1. Estimate: 表示各變數的估計係數。
2. Std. Error (SE) of Estimate: 估計標準誤，用來衡量估計回歸式的精準度。表示觀測值與回歸線間的分散或離散程度。
   \(\beta_{0}, \beta_{1}\)計算方式分別如下：
   $$SE(\beta_{0})^2 = \sigma^2[\frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^n (x_{i} – \bar{x})^2}],\space \space SE(\beta_{1})^2 = \frac{\sigma^2}{\sum_{i=1}^n (x_{i} – \bar{x})^2}$$
   其中，\(\sigma^2 = Var(\epsilon)\)表示隨機誤差的變異數。
   我們可以透過標準誤SE來計算95%信心水準係數區間如下：
   
   confint(object = modelFit2,level = 0.95) # 2.5 % 97.5 % # (Intercept) 0.03497014 0.10267304 # Temp 0.15044504 0.22094442 # Wind -0.13912079 -0.06470414

   1 2 3 4 5

   confint(object = modelFit2,level = 0.95) # 2.5 % 97.5 % # (Intercept) 0.03497014 0.10267304 # Temp 0.15044504 0.22094442 # Wind -0.13912079 -0.06470414

   
   從以上數據來說，在95%信心水準下\(\beta{1}(Temp)\)的係數區間[0.15,0.22]並沒有包括0，也因此我們可以推論說，當溫度升高100度時，log(臭氧濃度)會增加15-22個單位。
3. t-value: 衡量\(\beta{1}\)距離0有多少個標準差。t-value越大，p-value越小，表示該變數與目標變數關係顯著。t-value計算公式如下：
   $$t = \frac{\beta_{1}-0}{SE(\beta_{1})}$$
4. Pr(>|t|): 兩個預測變數p-value都小於0.05(95%信心水準)，表示與目標變數有顯著關係。

評估模型正確性(Model Accuracy)指標說明：

我們可以透過以下幾個數值指標來了解模型配適資料的適合度(goodness-of-fit)。

1. Residual standard error (RSE)：估計誤差標準差(estimate the standard deviation of error \(\epsilon\))。此模型誤差標準差為0.0458 (越小越好)，表示平均來說，實際觀測到的臭氧濃度會偏離回歸預測值約0.0458個單位(這是取log轉換後的數值）。
   誤差標準差衡量的是平均觀測值偏離回歸的程度。計算公式如下：
   $$RSE = \sqrt{\frac{1}{n-2} \sum_{i=1}^n (y_{i} – \hat{y}_{i})^2}$$
   RSE除了顯示在sumary(modelFit2)的底端，亦可透過以下sigma()函數取得。
   
   sigma(modelFit2) # [1] 0.04580728

   1 2	sigma(modelFit2) # [1] 0.04580728

   
   但畢竟RSE評估模型配適資料合適度是用實際觀測值偏離回歸預測值多少個Y單位，無從得知構成模型配適
2. R-squared (判定係數): 衡量的是預測變數共解釋了多少變異。\(R^2\)是一個介於0到1的數值，越大越好，且數值不受Y單位所影響。\(R^2\)計算公式如下：
   $$R^2 = 1 – \frac{SSR}{SST} = 1 – \frac{\sum_{i=1}^n (y_{i} – \hat{y}_{i})^2}{\sum_{i=1}^n (y_{i} – \bar{y}_{i})^2}$$
   除了可以從sumary(modelFit2)最後幾行得知\(R^2\)的資訊，我們亦可透過以下modelr套件中的rsqaure()函數取得。
   
   rsquare(modelFit2, data = completedData) # [1] 0.604058

   1 2	rsquare(modelFit2, data = completedData) # [1] 0.604058

   
   此模型的預測變數(Temp & Wind)共解釋了60%目標變數(Ozone)的變異。
3. Adjusted R-squared: 調整後的\(R^2\)有綜合考量變數總數(模型複雜度)的懲罰效果。此模型adj.\(R^2\)為 0.599(會比\(R^2\)來的小)，近似60%。
4. F-statistic: F-statistic test是用來檢測是否模型中至少有一預測變數係數不為0(當模型為複回歸時，此檢測結果特別重要)。F統計量的計算公式如下：
   $$F = \frac{TSS – RSS/p}{RSS/n-p-1}$$
   F統計量越大，p-value越顯著(p<0.05)。
   此模型的F統計量為114(越大越好)，對應的p-value為 p < 2.2e-16。

5. linear regression 線性回歸模型適合度診斷

5-1. 基本假設檢定

檢視模型配適度一個重要的步驟，就是透過視覺化殘差分佈，來檢視模型基本假設是否皆符合。

由下圖可以發現：

• 左上：殘差變異數不符合均一性。(理想上，希望紅色線條要接近水平線並接近0)
• 右上：殘差分佈不符合常態。（理想上，希望殘差值分佈不在45度線上）
• 左下：標準化殘差變異數不符合均一性。(理想上，希望紅色線條要接近水平線並接近0)
• 右下：檢視是否有對模型影響力過高的資料點。（若是有，資料點將會超過cook’s distance虛線外）

綜合以上，會發現線性模型似乎不適合用來配適此組資料。

5-2. gvlma package

另外，也可以使用gvlma套件來檢定模型是否符合假設。可以發現有部分線性迴歸模型的基礎假設是不成立的。(不過要注意的是，gvlma套件結果相關說明的資料不是很充足。)

5-3. shapiro test常態分佈檢定

另外，我們亦使用shapiro test 來看鎖定變數是否成常態分佈(H0: 分佈等於常態分佈)。可以發現除了Wind符合常態分配，另外兩個變數都離常態分佈有些距離。因此除非進一步將變數有效轉換以符合常態，否則不適合用線性迴規模性來進行資料配適。

6. 模型驗證Validation

6-1. 常見模型效果評估指標

STATISTIC	CRITERION
R-Squared	Higher the better (> 0.70)。表示解釋變異量佔總變異量的百分比。
Adj R-Squared	Higher the better。表示調整後(對模型複雜度加入懲罰)解釋變異量佔總變異量的百分比。
F-Statistic	Higher the better。
Std. Error (標準誤)	Closer to zero the better。 (表抽樣結果與母體的偏差)
t-statistic	Should be greater 1.96 for p-value to be less than 0.05
AIC	Lower the better
BIC	Lower the better
Mallows cp	Should be close to the number of predictors in model
MAPE (Mean absolute percentage error)	Lower the better
MSE (Mean squared error)	Lower the better
Min_Max Accuracy => mean(min(actual, predicted)/ max(actual, predicted))	Higher the better

6-2. 模型驗證(80%Training20%Testing)

將資料隨機區分為訓練資料集(80%)與測試資料集(20%)

使用訓練及資料建置模型 & 使用測試集資料預測結果。並摘要訓練的模型效果summary(mod)。

檢視模型預測的效果。

計算真實數據與預測數據的相關係數：越高越好。

計算AIC, BIC => 越低越好。

計算模型預測正確率min_max_accuracy： 70.4% = > 越高越好。

計算模型預測錯誤率 MAPE (Mean absolute percentage error)： 56% => 越低越好。

6-3. 加強版模型驗證: K-fold Cross Validation

使用k-fold交叉驗證。

查看交叉驗證的平均平方誤差(MSE):越低越好。

交叉5次模型配適結果圖如下。小圖示為模型預測值，大圖示則為真實資料值。可以發現第五次模型擬和的結果與其他四次有較大差異。

7. 結論

我們可以發現，由於airquality資料型態不接近理想的常態分佈，故使用Linear Regression 線性迴歸模型做資料配適不是很合適(根據「回歸模型適合度診斷」)。如果遇到資料分配非常態時，建議可以考慮其他模型如決策樹（不受限母體分配）。

---

更多線性回歸模型應用：

1. 開發一個免費App能賺多少錢? 靠 AdMod 廣告月收3萬實例分享

更多統計模型學習筆記連結：

1. Regularized Regression | 正規化迴歸 – Ridge, Lasso, Elastic Net | R語言
2. Logistic Regression 羅吉斯迴歸 | part1 – 資料探勘與處理 | 統計 R語言
3. Logistic Regression 羅吉斯迴歸 | part2 – 模型建置、診斷與比較 | R語言
4. Decision Tree 決策樹 | CART, Conditional Inference Tree, Random Forest
5. Regression Tree | 迴歸樹, Bagging, Bootstrap Aggregation | R語言
6. Random Forests 隨機森林 | randomForest, ranger, h2o | R語言
7. Gradient Boosting Machines GBM | gbm, xgboost, h2o | R語言
8. Hierarchical Clustering 階層式分群 | Clustering 資料分群 | R統計
9. Partitional Clustering | 切割式分群 | Kmeans, Kmedoid | Clustering 資料分群
10. Principal Components Analysis (PCA) | 主成份分析 | R 統計

---

參考文章連結：

1. Linear Regression
2. Linear Regression using AirQuality Dataset
3. Imputing missing data with R; MICE package
4. Assumptions of Linear Regression
5. Chapter 12: Regression: Basics, Assumptions, & Diagnostics
6. 複回歸分析
7. Linear Regression (2)