Correlation & Regression Analysis (I)
gcchen
2019/6/11
相關與迴歸分析 I (相關篇)
* 隨機變數 \(X,Y\) 的 共變異係數 定義為 \[Cov(X,Y) = E(XY) - E(X)E(Y)\]
舉例而言(請觀察與思考):
- 以下兩組資料：\(X_1=(1,2,3,4,5,6,7,8,9,10), X_1=(1,2,3,4,5,6,7,8,9,10)\)，計算可得母體共變異數 \(Cov(X_1,X_1)= \sigma _{X_1X_1}= \frac{1+4+9+\cdots+81+100}{10} - \frac{55}{10} * \frac{55}{10} = 8.25\) (以本例而言，即為母體變異數)，或樣本共變異數 \(Cov(X_1,X_1)= s_{X_1 X_1} = \frac{\displaystyle \sum_{i=1}^{n} x_i^2}{n-1} -\left(\frac{n}{n-1} \right) * \bar{x}^2 = \frac{1+4+9+\cdots+81+100}{9} - \frac{10}{9} * 5.5^2 \approx 9.166667\) (以本例而言，即為樣本變異數)，可透過以下R程式計算：
X_1 <- c(1,2,3,4,5,6,7,8,9,10)
cov(X_1,X_1)## [1] 9.166667將兩組資料視為10個點的 \(XY\) 座標繪製散佈圖
plot(X_1,X_1)
- 以下兩組資料：\(X_1=(1,2,3,4,5,6,7,8,9,10), Y_1=(2,4,6,8,10,12,14,16,18,20)\)，計算可得 \(Cov(X_1,Y_1)= \sigma_{X_1Y_1}= \frac{2+8+18+\cdots+162+200}{10} - \frac{55}{10} * \frac{110}{10} = 16.50\)，或樣本共變異數 \(Cov(X_1,Y_1)= s_{X_1 Y_1} = \frac{\displaystyle \sum_{i=1}^{n} x_i y_i}{n-1} -\left(\frac{n}{n-1} \right) * \bar{x} \bar{y} = \frac{2+8+18+\cdots+162+200}{9} - \frac{10}{9} * 5.5* 11 \approx 18.33333\)，可透過以下R程式計算：
# (已定義) X_1 <- c(1,2,3,4,5,6,7,8,9,10)
Y_1 <- c(2,4,6,8,10,12,14,16,18,20)
cov(X_1,Y_1)## [1] 18.33333將兩組資料視為10個點的 \(XY\) 座標繪製散佈圖
plot(X_1,Y_1)
- 以下兩組資料：\(X_1=(1,2,3,4,5,6,7,8,9,10), Y_2=(10,9,8,7,6,5,4,3,2,1)\)，計算可得 \(Cov(X_1,Y_2) = \sigma_{X_1Y_2} = \frac{10+18+24+\cdots+24+18+10}{10} - \frac{55}{10} * \frac{55}{10} = -8.25\)，或樣本共變異數 \(Cov(X_1,Y_2)= s_{X_1 Y_2} = \frac{\displaystyle \sum_{i=1}^{n} x_i y_i}{n-1} -\left(\frac{n}{n-1} \right) * \bar{x} \bar{y} = \frac{10+18+24+\cdots+18+10}{9} - \frac{10}{9} * 5.5*5.5 \approx -9.166667\)，可透過以下R程式計算：
# (已定義) X_1 <- c(1,2,3,4,5,6,7,8,9,10)
Y_2 <- c(10,9,8,7,6,5,4,3,2,1)
cov(X_1,Y_2)## [1] -9.166667將兩組資料視為10個點的 \(XY\) 座標繪製散佈圖
plot(X_1,Y_2)
- 以下兩組資料：\(X_1=(1,2,3,4,5,6,7,8,9,10), Y_3=(10,10,10,10,10,10,10,10,10,10)\)，計算可得 \(Cov(X_1,Y_3) = \sigma_{X_1 Y_3} = \frac{10+20+30+\cdots+80+90+100}{10} - \frac{55}{10} * \frac{100}{10} = 0\)，或樣本共變異數 \(Cov(X_1,Y_3)= s_{X_1 Y_3} = \frac{\displaystyle \sum_{i=1}^{n} x_i y_i}{n-1} -\left(\frac{n}{n-1} \right) * \bar{x} \bar{y} = \frac{10+20+30+\cdots+80+90+100}{9} - \frac{10}{9} * 5.5*10.0 = 0\)，可透過以下R程式計算：
# (已定義) X_1 <- c(1,2,3,4,5,6,7,8,9,10)
Y_3 <- c(10,10,10,10,10,10,10,10,10,10)
cov(X_1,Y_3)## [1] 0將兩組資料視為10個點的 \(XY\) 座標繪製散佈圖
plot(X_1,Y_3)
* 隨機變數 \(X,Y\) 的 相關係數 定義為：
母體相關係數：
\[\rho_{XY} = \frac{\sigma_{XY}}{\sigma_{X} \sigma_{Y}} = \frac{Cov(X,Y)}{\sigma_{X} \sigma_{Y}} = \displaystyle \frac{ \displaystyle \frac{ \displaystyle \sum_{i=1}^{N} (X_i - \mu_X)(Y_i - \mu_Y)}{N}}{\displaystyle \sqrt{\frac{\displaystyle \sum_{i=1}^{N} (X_i - \mu_X)^2}{N}} \sqrt{\frac{ \displaystyle \sum_{i=1}^{N} (Y_i - \mu_Y)^2}{N}} }= \frac{E(XY)-E(X)E(Y)}{\sigma_{X} \sigma_{Y}} = \frac{E[(X-\mu_X)(Y-\mu_Y)]}{\sigma_{X} \sigma_{Y}} = E \left[ \left( \frac{X-\mu_X}{\sigma_{X}} \right) \left( \frac{Y-\mu_Y}{\sigma_{Y}} \right) \right]\]
樣本相關係數：
\[r_{XY} = \frac{S_{XY}}{S_{X} S_{Y}} = \displaystyle \frac{ \displaystyle \frac{ \displaystyle \sum_{i=1}^{n} (X_i - \bar{X})(Y_i - \bar{Y})}{n-1}}{\displaystyle \sqrt{\frac{\displaystyle \sum_{i=1}^{n} (X_i - \bar{X})^2}{n-1}} \sqrt{\frac{ \displaystyle \sum_{i=1}^{n} (Y_i - \bar{Y})^2}{n-1}} } = \displaystyle \frac{ \displaystyle \sum_{i=1}^{n} (X_i - \bar{X})(Y_i - \bar{Y})}{ \sqrt{\displaystyle \sum_{i=1}^{n} (X_i - \bar{X})^2} \sqrt{ \displaystyle \sum_{i=1}^{n} (Y_i - \bar{Y})^2}} = \frac{1}{n-1} \sum_{i=1}^{n} \left( \frac{X_i-\bar{X}}{S_{X}} \right) \left( \frac{Y_i-\bar{Y}}{S_{Y}} \right) \]
下表為西元 1951 與 1980 年的國內各業生產毛額(資料來源:行政院主計處，單位百萬元)
| 業別 | 1951 | 1980 |
|---|---|---|
| 1.農業 | 3,997 | 119,086 |
| 2.礦業 | 162 | 15,810 |
| 3.製造業 | 1,701 | 525,623 |
| 4.水電燃氣業 | 137 | 36,357 |
| 5.營造業 | 447 | 87,963 |
| 6.商業 | 1,573 | 189,156 |
| 7.運輸倉儲及通信業 | 493 | 87,709 |
| 8.金融保險不動產及工商服務業 | 1,237 | 154,789 |
| 9.社會服務及個人服務業 | 532 | 65,483 |
利用1951(x1)與1980(x2)的國內各業生產毛額資料繪製散佈圖，並計算相關係數：
x1 <- c(3997,162,1701,137,447,1573,493,1237,532)
x2 <- c(119086,15810,525623,36357,87963,189156,87709,154789,65483)
plot(x1,x2)
cor(x1,x2)## [1] 0.346681將農業(衰退?)與製造業(蓬勃?)資料移除，再計算相關，得到極高的相關係數!
x1 <- x1[-c(1,3)]
x2 <- x2[-c(1,3)]
plot(x1,x2)
cor(x1,x2)## [1] 0.9735761新新高中學生每年都會舉行越野賽跑,繞新竹十八尖山一圈,小劉在正式測驗前一個月每天下午放學後都去跑一圈做練習,並記錄時間(單位:分鐘)及心跳數(每分鐘心跳次數),得到 23 筆資料如下:
| 時間 | 心跳數 | 時間 | 心跳數 |
|---|---|---|---|
| 26.60 | 136 | 27.37 | 148 |
| 26.00 | 148 | 27.57 | 144 |
| 26.35 | 148 | 27.43 | 136 |
| 27.62 | 132 | 27.57 | 144 |
| 27.68 | 124 | 28.17 | 136 |
| 27.28 | 132 | 27.72 | 128 |
| 27.97 | 139 | 26.12 | 152 |
| 25.93 | 156 | 35.72 | 124 |
| 26.75 | 140 | 26.72 | 140 |
| 26.70 | 144 | 26.05 | 152 |
| 26.85 | 148 | 26.13 | 146 |
| 28.05 | 124 |
畫出一個散佈圖。以 X 軸表示時間,Y 軸表示心跳數。並計算相關係數:
pulse<-read.csv("pulse.csv",header=T)
attach(pulse)
plot(時間,心跳數)
cor(時間,心跳數)## [1] -0.5953083將大於等於30分鐘的資料移除(小劉當天身體狀況有可能比較特殊，或時間有延誤?)，再計算相關，可得到較高的相關係數
cor(時間[時間<30],心跳數[時間<30])## [1] -0.7387117若將時間的單位改成秒，並不會改變相關係數的大小
cor(時間*60,心跳數)## [1] -0.5953083以下是透過 最小平方法 決定出來的 簡單迴歸方程式：
summary(lm(心跳數~時間))##
## Call:
## lm(formula = 心跳數 ~ 時間)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.252 -5.177 1.912 5.622 11.683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 219.3602 23.4177 9.367 6.01e-09 ***
## 時間 -2.8941 0.8524 -3.395 0.00273 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.788 on 21 degrees of freedom
## Multiple R-squared: 0.3544, Adjusted R-squared: 0.3236
## F-statistic: 11.53 on 1 and 21 DF, p-value: 0.002729亦即，y = 219.3602-2.8941 x ，其中 y 代表心跳數， x 代表時間(以分鐘為單位)