使用 “ISwR” package 中的 cystfibr 資料 (囊腫性纖維化病人的 資料 (7–23 years old).) 作為範例, 先將資料載入
EDA
Summary
稍微看一下 summary, 可以看到, sex 應該為 Factor, 其餘不存在太多問題
## age sex height weight bmp
## Min. : 7.00 Min. :0.00 Min. :109.0 Min. :12.9 Min. :64.00
## 1st Qu.:11.00 1st Qu.:0.00 1st Qu.:139.0 1st Qu.:25.1 1st Qu.:68.00
## Median :14.00 Median :0.00 Median :156.0 Median :37.2 Median :71.00
## Mean :14.48 Mean :0.44 Mean :152.8 Mean :38.4 Mean :78.28
## 3rd Qu.:17.00 3rd Qu.:1.00 3rd Qu.:174.0 3rd Qu.:51.1 3rd Qu.:90.00
## Max. :23.00 Max. :1.00 Max. :180.0 Max. :73.8 Max. :97.00
## fev1 rv frc tlc pemax
## Min. :18.00 Min. :158.0 Min. :104.0 Min. : 81 Min. : 65.0
## 1st Qu.:26.00 1st Qu.:188.0 1st Qu.:127.0 1st Qu.:101 1st Qu.: 85.0
## Median :33.00 Median :225.0 Median :139.0 Median :113 Median : 95.0
## Mean :34.72 Mean :255.2 Mean :155.4 Mean :114 Mean :109.1
## 3rd Qu.:44.00 3rd Qu.:305.0 3rd Qu.:183.0 3rd Qu.:128 3rd Qu.:130.0
## Max. :57.00 Max. :449.0 Max. :268.0 Max. :147 Max. :195.0Scatter plot with smooth
將 response 對變數一一畫 scatter plot, straight line 以及 lowess line, 看看是否有線性相關
可以看到, 有些部分資料呈現線性, 有些則無, 有線性的變數我們就放至 linear part, 有 non - linear 趨勢的我們就放入 smoothing part
distribution checking
pemax 變數是我們的 response, 假若我們要做檢定或是信賴區間, 那就必須要對 \(\epsilon\) 有分布假設, 一般假設為 \(\epsilon \sim N(\mathbf{0}, \sigma^{2}\mathbf{I})\), 而 \(Y\) 就會繼承其分布, 也就是 \(Y \sim N(\mathbf{X}\beta, \sigma^{2}\mathbf{I})\), 所以我們現在可以先看看 \(Y\) 的 density plot, 以及 shapiro.test 是否有過
##
## Shapiro-Wilk normality test
##
## data: data$pemax
## W = 0.891, p-value = 0.01175透過圖片可以發現有明顯右偏, 檢定結果也被拒絕, 所以我們試著做 log 轉換
##
## Shapiro-Wilk normality test
##
## data: log(data$pemax)
## W = 0.94779, p-value = 0.2234做完 log 轉換就可以發現, density plot 變得更加像是 normal, 檢定也未被拒絕, 所以接下來都用 log 轉換後的 response
若需要較客觀的轉換建議, 可以做 box cox transfrom
\[ y_{i}^{(\lambda)}=\left\{\begin{array}{ll} \frac{y_{i}^{\lambda}-1}{\lambda} & \text { if } \lambda \neq 0 \\ \ln y_{i} & \text { if } \lambda=0 \end{array}\right. \]
透過信賴區間得到 \(\lambda\) 估計值, 進一步推斷建議何種轉換
可以看到, 0 落入信賴區間, 所以 log 轉換是合理的
correlation plot
接著我們來看看各個變數間的相關性
可以看到在某些變數間存在強烈的線性相關, 比較能夠解釋的是隨年齡增長, 體重以及身高都會相對比較高, 所以符合直覺, 而性別應該要與身高體重有強烈相關性, 但是在這裡卻沒有, 可能是收集到的資料剛好沒有這個現象
透過觀察 correlation plot 我們也可以試著做看看 factor analysis, 但我們先 focus 在 linear model 上
根據上面的 correlation plot 選變數做 linear model, 若相關性高, 那假若都放入模型中, 可能會導致兩個變數試圖來做同一件事情, 例如我放年齡, 身高, 體重進入變數中, \(R^{2}\) 並不會有太多的上升, 卻會造成多浪費一個自由度, 所以這裡先用 age, fev1, bmp, rv 來作為變數配適 linear model
fit model
配適線性模型 \[ log(pemax) = \beta_0 + \beta_{1}(age) + \beta_{2}(fev1)+\beta_{3}(bmp) + \beta_{4}(rv) + \epsilon \]
##
## Call:
## lm(formula = log(pemax) ~ age + fev1 + bmp + rv, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39311 -0.18781 0.03121 0.17331 0.28520
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.7169970 0.5777809 6.433 2.83e-06 ***
## age 0.0374931 0.0110201 3.402 0.00283 **
## fev1 0.0129049 0.0057433 2.247 0.03609 *
## bmp -0.0038946 0.0048992 -0.795 0.43598
## rv 0.0009746 0.0008532 1.142 0.26685
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2246 on 20 degrees of freedom
## Multiple R-squared: 0.4814, Adjusted R-squared: 0.3777
## F-statistic: 4.642 on 4 and 20 DF, p-value: 0.008174可以看到 bmp 以及 rv 並不顯著, 但會將它保留, 若捨棄此兩個變數, 會導致模型假設被違反 (you can try it)
mean structure
Component + Residual Plot
參考 Component + Residual Plot, 此 plot 為將 \(\hat{\epsilon} + X_{[,i]}\hat\beta_{i}\) 對 \(X_{[,i]}\) 作圖
\[\begin{align*} Y &= \mathbf{X}\beta + \epsilon\\ &= \mathbf{X}\hat{\beta} + \hat{\epsilon}\\ &= \mathbf{X}_{[,-i]}\hat{\beta}_{-i} + X_{[,i]}\hat{\beta}_{i} + \hat{\epsilon} \end{align*}\]
可以看到 \(X_{[,i]}\hat{\beta}_{i} + \hat{\epsilon}\) 此項主要顯示前面 fit mean structure 後, 沒 fit 好的部分, 假若看此圖中, smoothing scatter plot 有 non-linear 的趨勢, 就要想辦法對 \(X_{[,i]}\) 做轉換, 注意一次只要更動一次變數即可, 更動一個變數, 會導致其他變數的 Component + Residual Plot 改變, 所以一次做一個變數再確認, 重複此事直到對所有變數做的 Component + Residual Plot 都沒有 non-linear structure 即可罷休
可以看到, 很幸運的是, 看起來都沒有 non-linear structure
additive variable plot
error structure
diagnostic plots
來看看模型的 diagnostic plot, 由左上角至右下角
- residual 沒有明顯的 mean structure, 也沒有 non-constant variance 的傾向
- Normal QQ plot 看起來差強人意, 且 test 未過關, 但還算可以接受
- standardized 後的 residual 看起來也沒有明顯 mean structure 及 non-constant variance
- 沒有明顯的 influence points 存在
##
## Shapiro-Wilk normality test
##
## data: model$residuals
## W = 0.91708, p-value = 0.04397將 residual 對各個變數作圖, 也可以發現對各個變數沒有明顯的 mean structure, 也沒有 non-constant variance 的傾向
同樣在 residual 對 fitted value 作圖也可以看到沒有明顯的 mean structure, 也沒有 non-constant variance 的傾向
leverage
接下來來看看資料中是否有離資料中心較遠的點, 我們稱這種觀察值為 leverage 大的觀察值, leverage (\(h_{i}\)) 為 hat matrix 的對角部分, 我們現在假設 \(e_{i}\) 為第 i 個component 為 1, 其餘為零的向量, 則
\[
h_{i} = \mathbf{H}_{ii}\\
\mathbf{H} = X(X^{T}X)^{-1}X^{T}\\
\text{redefine }X = [\mathbf{1}, X - \bar{X}]_{(n,p)}\text{ (span the same space)}\\
\mathbf{1}^{T}X_{[,-1]} = 0\\
H = \begin{pmatrix}\mathbf{1}& X_{[,-1]}\end{pmatrix}
\begin{pmatrix}
\mathbf{1}^{T}\mathbf{1} & 0 \\
0 & X_{[,-1]}^{T}X_{[,-1]}
\end{pmatrix}^{-1}\begin{pmatrix}\mathbf{1}\\ X_{[,-1]}\end{pmatrix}\\
= \frac{1}{n} + X_{[,-1]}\left(X_{[,-1]}^{T}X_{[,-1]}\right)^{-1}X_{[,-1]}^{T}\\
H_{ii} = e_{i}^{T}\left( \frac{1}{n} + X_{[,-1]}\left(X_{[,-1]}^{T}X_{[,-1]}\right)^{-1}X_{[,-1]}^{T} \right)e_{i}\\
= \frac{1}{n} + X_{[i\ ,-1]}\left(X_{[,-1]}^{T}X_{[,-1]}\right)^{-1}X_{[i\ ,-1]}^{T}\\
= \frac{1}{n} + (X_{i} - \bar{X})\left(X_{[,-1]}^{T}X_{[,-1]}\right)^{-1}(X_{i} - \bar{X})^{T}\\
= \frac{1}{n} + \frac{1}{n-1}(X_{i} - \bar{X})\hat\Sigma^{-1}(X_{i} - \bar{X})^{T}\\
\]
可以看到, 假若 \(X_{i} - \bar{X}\) 較大, 也就意味著該點距離資料中心較遠, 如此一來, leverage 就會比較大, 一般來說, 當 leverage 大於 \(\frac{2p}{n}\) 我們就需要多加留意, 他可能是 outlier 或是 influence point, 接著我們來看看那些點距離資料中心較遠, 是否有 leverage 大於 \(\frac{2p}{n}\) 的狀況
可以看到沒有甚麼大問題, 各個觀察值都在 \(\frac{2p}{n}\) 之下
standardized residual
假若 \(\epsilon \sim (\mathbf{0},\sigma^{2}\mathbf{I})\), 則 \(\hat\epsilon = (\mathbf{I - H})\epsilon \sim (\mathbf{0}, \sigma^{2}(\mathbf{I - H}))\), 若我們想要透過 \(\hat{\epsilon}\) 去反推得 \(\epsilon\), 並不能直接透過反矩陣運算還原, 因為 \(\mathbf{I - H}\) 是 singular 的, 反矩陣不存在 \[ Var(e_{i}^{T}\hat\epsilon) = Var(\hat{\epsilon_{i}}) = e_{i}^{T}(\mathbf{I - H})e_{i}\sigma^{2} = (1-h_{ii})\sigma^{2}\\ Var(\frac{\hat{\epsilon}_{i}}{\sqrt{1-h_{ii}}}) = \sigma^{2} \]
這時我們定義 internal standardized residual : \(r_{i}\) 為 :
\[ r_{i} = \frac{\hat{\epsilon}_{i}}{\sqrt{1-h_{ii}}\ \hat\sigma^{2}}\\ Var(r_{i}) \approx 1 \]
現在我們將 internal standardized residual 畫出, 並計算其變異數
## Variance of internal standardized residual is : 1.020818可以發現, standardized residual 沒有明顯 mean structure 以及 non-constant variance 的跡象, 其 Variance 也在 1 附近
outliers
當一個觀察點, 和我們配適的模型不吻合, 我們稱其為 outlier, 通常他會有比較大的 residual, 透過計算 jacknife residuals (external standardized residuals) 我們可以偵測 outliers.
jacknife residual 的想法就是先將第 i 筆資料暫時移除, 用剩下的資料去做回歸, 再透過該模型去預測第 i 筆資料, 最後將移除的第 i 筆資料和預測的數值相減後除上預測標準差, 就會是 jacknife residual, 也就是 :
\[ t_{i} = \frac{y_{i} - \hat{y}_{(i)}}{\sqrt{1+x_{i}^{T}\left(X_{(i)}^{T}X_{(i)}\right)^{-1}x_{i}\ \sigma^{2}_{(i)}}} \sim t(n-p-1) \]
其中下標 (i) 代表已移除第 i 筆資料的資料, 其又可以寫成
\[ t_{i}=\frac{\hat{\varepsilon}_{i}}{\left(1-h_{i}\right)^{1 / 2} \hat{\sigma}_{(i)}}=r_{i}\left((n-p-1) /\left(n-p-r_{i}^{2}\right)\right)^{1 / 2} \]
可以看到其和 internal standardized residuals 有關係
現在我們得到 n 個 jacknife residual, 現在我們可以個別對每筆資料做檢定(\(|t_{i}| > t_{\alpha/2}(n-p-1)\)), 看是否為 outliers, 但是當我們想做的事情, 是去判對此資料集中是否有離群值, 我們就必須考慮到 multipule testing 的問題, 所以現在就必須去看 \(|t_{i}| > t_{\alpha/2n}(n-p-1)\) (bonferroni correction)
現在我們計算該資料中的各個 jacknife residual
可以看到並未有明顯的 outlier
Some problems of outliers
當 n 大時, \(\alpha/2n\) 幾乎為 0 , 導致可能是 outliers 的點都沒被標出來, 可以用其他方法去做 multipule testing (ex : Hotelling \(T^{2}\))
假設資料集中存在多筆 outlier , outlier 就會包庇其他離群值(因為把某 outlier 移除後, 模型中尚存在其他 outlier 去配適模型, 導致計算出來的 jacknife residual 很小), 或許可以一次移除較多 outlier (先看 scatter plot)
有時候直接移除 outlier 並不是一個聰明的做法, outlier 可能隱含重要資訊
influence point
在配適模型時, 我們會透過觀察值去估計 \(\beta\) 以及 \(\sigma^{2}\) 這兩個參數, 有時候加入一些觀察值, 可能造成這兩個參數的估計變動大, 這種點就叫做 influential point, 當 leverage 很高, 同時他又是一個 outlier 時, 他就會對模型的參數估計造成很大的影響, 以下提供示意圖
透過計算 Cooks distance 我們可以得知哪個點為 influence point,
\[\begin{align*} D_{i} &=\frac{\left(\hat{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}_{i}\right)^{\mathrm{T}}\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{X}\right)\left(\hat{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}_{i}\right)}{\left(p \hat{\sigma}^{2}\right)} \\ &=\frac{\left(\hat{Y}-\hat{Y}_{(i)}\right)^{\mathrm{T}}\left(\hat{Y}-\hat{Y}_{(i)}\right)}{\left(p \hat{\sigma}^{2}\right)} \\ &=(1 / p) r_{i}^{2}\left(h_{i} /\left(1-h_{i}\right)\right) \end{align*}\]
透過最下面的等號不難看出其同時與 internal standardized residual 以及 leverage 有關, 通常大於 1, 0.5 我們就需要多加留意
對此資料計算各個觀察值的 cook distance
可以發現沒有明顯 influence point