UNIT09C：多重線性回歸 Multiple Linear Regression
中山大學管理學院 卓雍然
2024-04-17 15:09:48
pacman::p_load(dplyr, ggplot2, corrplot)[A] 多元線性回歸 - lm()
§ A1 讀進資料
全球暖化資料集
D = read.csv('data/climate_change.csv')
head(D) Year Month MEI CO2 CH4 N2O CFC.11 CFC.12 TSI Aerosols Temp
1 1983 5 2.556 345.96 1638.6 303.68 191.32 350.11 1366.1 0.0863 0.109
2 1983 6 2.167 345.52 1633.7 303.75 192.06 351.85 1366.1 0.0794 0.118
3 1983 7 1.741 344.15 1633.2 303.80 192.82 353.73 1366.3 0.0731 0.137
4 1983 8 1.130 342.25 1631.3 303.84 193.60 355.63 1366.4 0.0673 0.176
5 1983 9 0.428 340.17 1648.4 303.90 194.39 357.46 1366.2 0.0619 0.149
6 1983 10 0.002 340.30 1663.8 303.97 195.17 359.17 1366.1 0.0569 0.093現代的資料分析過程之中，我們通常不會將所有的資料拿來做模型，而是會先保留一部分的資料下來，用來估計模型的準確性。所以一開始我們會將全部資料(D)分割為訓練資料(TR)與測試資料(TS)
TR = subset(D, Year <= 2006) # Train Data
TS = subset(D, Year > 2006) # Test Data§ A2 建立模型
使用訓練資料(TR)做模型
m1 = lm(Temp~MEI+CO2+CH4+N2O+CFC.11+CFC.12+TSI+Aerosols,TR)§ A3 模型摘要
options(digits=4, scipen=10)
summary(m1)
Call:
lm(formula = Temp ~ MEI + CO2 + CH4 + N2O + CFC.11 + CFC.12 +
TSI + Aerosols, data = TR)
Residuals:
Min 1Q Median 3Q Max
-0.2589 -0.0591 -0.0008 0.0565 0.3243
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -124.594260 19.886800 -6.27 0.0000000014310 ***
MEI 0.064205 0.006470 9.92 < 2e-16 ***
CO2 0.006457 0.002285 2.83 0.00505 **
CH4 0.000124 0.000516 0.24 0.81015
N2O -0.016528 0.008565 -1.93 0.05467 .
CFC.11 -0.006630 0.001626 -4.08 0.0000595728765 ***
CFC.12 0.003808 0.001014 3.76 0.00021 ***
TSI 0.093141 0.014755 6.31 0.0000000010959 ***
Aerosols -1.537613 0.213252 -7.21 0.0000000000054 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0917 on 275 degrees of freedom
Multiple R-squared: 0.751, Adjusted R-squared: 0.744
F-statistic: 104 on 8 and 275 DF, p-value: <2e-16\(y\) 的 預測值：
- \(\hat{y} = b_0 + b_1 x_1 + b_2 x_2 + b_3 x_3 + b_4 x_4 + ...\)
簡單講，對每一個自變數而言：
- \(b_i\) : 係數的估計值代表 \(x_i\) 對 \(y\) 的邊際效果
- \(\sigma_i\) : 係數的標準差代表該係數的不確定性
- \(p_i\) : \(p\)-value 大致上代表 \(x_i\) 與 \(y_i\) 之間沒有關係的機率， 這個機率小於顯著水準，我們就說 \(x_i\) 對 \(y\) 的關係是顯著的
嚴格來講：
- 問題很多 … 整個管統、多變量、和計量經濟大概都在處理這些問題
幸好：
- 在商業數據分析裡面通常我們比較關心目標變數的預測值(\(\hat{y}\))，這部分比較沒有問題
[B] 係數的分布(機率密度函數)
vx = c(3:5,7); cx = summary(m1)$coef
par(cex=0.8)
plot(0,0,xlim=c(-0.04,0.02),ylim=c(0,800),pch=20,
xlab="coefficients: the estimated value of b's",
ylab='probability density',
main='Probability Density Function of Coefficients')
abline(h=seq(0,800,100),v=seq(-0.04,0.03,0.0025),col='lightgray',lty=3)
for(i in vx) curve(dnorm(x,cx[i,1],cx[i,2]),add=T,col=i,lwd=2,n=1000)
abline(v=0,col='red')
legend("topleft",col=vx,lwd=2,legend=paste0(
"b",vx," ( ",rownames(cx)[vx]," )", c("**","","。","***") ) )
summary(m1)$coef[vx,] Estimate Std. Error t value Pr(>|t|)
CO2 0.006457 0.0022846 2.8264 0.0050525
CH4 0.000124 0.0005158 0.2405 0.8101456
N2O -0.016528 0.0085649 -1.9297 0.0546693
CFC.12 0.003808 0.0010135 3.7573 0.0002097
🗿 討論問題：
1.依據這個模型，在CO2、CH4、N2O、CFC12之中：
■ 哪一種氣體的效果(關係)估計值最大？ 它的效果是顯著的嗎？
■
哪一種氣體的效果最顯著？ 它的效果估計值是最大的嗎？
■
除了係數的估計值之外，還有什麼因素會影響自變數的顯著性？
2.如果樣本變大的話：
■
係數的估計值會(a)變大、(b)變小、或(c)不變？
■
係數的估計標準差會(a)變大、(b)變小、或(c)不變？
■ 自變數的\(p\)-value會(a)變大、(b)變小、或(c)不變？
■ 自變數的顯著性會(a)變大、(b)變小、或(c)不變？
[C] 自變數之間的相關性、複回歸的共線性問題
🌻 自變數之間的相關性過高，可能會使係數的估計值有很大的偏差
cor(TR[3:10]) %>% round(2) MEI CO2 CH4 N2O CFC.11 CFC.12 TSI Aerosols
MEI 1.00 -0.04 -0.03 -0.05 0.07 0.01 -0.15 0.34
CO2 -0.04 1.00 0.88 0.98 0.51 0.85 0.18 -0.36
CH4 -0.03 0.88 1.00 0.90 0.78 0.96 0.25 -0.27
N2O -0.05 0.98 0.90 1.00 0.52 0.87 0.20 -0.34
CFC.11 0.07 0.51 0.78 0.52 1.00 0.87 0.27 -0.04
CFC.12 0.01 0.85 0.96 0.87 0.87 1.00 0.26 -0.23
TSI -0.15 0.18 0.25 0.20 0.27 0.26 1.00 0.05
Aerosols 0.34 -0.36 -0.27 -0.34 -0.04 -0.23 0.05 1.00cor(TR[3:10]) %>% corrplot.mixed(tl.cex=0.7, tl.col='black')
[D] 模型選擇、挑選變數
🌻 所以我們常把相關性太高的自變數從模型中移除
m2 = lm(Temp~MEI+N2O+TSI+Aerosols,TR) # 手動挑選自變數
summary(m2)
Call:
lm(formula = Temp ~ MEI + N2O + TSI + Aerosols, data = TR)
Residuals:
Min 1Q Median 3Q Max
-0.2792 -0.0598 -0.0060 0.0567 0.3419
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -116.22686 20.22303 -5.75 0.00000002373584 ***
MEI 0.06419 0.00665 9.65 < 2e-16 ***
N2O 0.02532 0.00131 19.31 < 2e-16 ***
TSI 0.07949 0.01488 5.34 0.00000018937322 ***
Aerosols -1.70174 0.21800 -7.81 0.00000000000012 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0955 on 279 degrees of freedom
Multiple R-squared: 0.726, Adjusted R-squared: 0.722
F-statistic: 185 on 4 and 279 DF, p-value: <2e-16🌻
R語言裡面也有一個自動挑選變數的內建函數step()
m3 = step(m1) # 自動挑選自變數Start: AIC=-1348
Temp ~ MEI + CO2 + CH4 + N2O + CFC.11 + CFC.12 + TSI + Aerosols
Df Sum of Sq RSS AIC
- CH4 1 0.000 2.31 -1350
<none> 2.31 -1348
- N2O 1 0.031 2.34 -1346
- CO2 1 0.067 2.38 -1342
- CFC.12 1 0.119 2.43 -1336
- CFC.11 1 0.140 2.45 -1333
- TSI 1 0.335 2.65 -1312
- Aerosols 1 0.437 2.75 -1301
- MEI 1 0.828 3.14 -1263
Step: AIC=-1350
Temp ~ MEI + CO2 + N2O + CFC.11 + CFC.12 + TSI + Aerosols
Df Sum of Sq RSS AIC
<none> 2.31 -1350
- N2O 1 0.031 2.35 -1348
- CO2 1 0.067 2.38 -1344
- CFC.12 1 0.130 2.44 -1337
- CFC.11 1 0.139 2.45 -1335
- TSI 1 0.335 2.65 -1314
- Aerosols 1 0.440 2.75 -1303
- MEI 1 0.831 3.14 -1265summary(m3)
Call:
lm(formula = Temp ~ MEI + CO2 + N2O + CFC.11 + CFC.12 + TSI +
Aerosols, data = TR)
Residuals:
Min 1Q Median 3Q Max
-0.2577 -0.0599 -0.0010 0.0559 0.3220
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -124.515178 19.850113 -6.27 0.0000000013651 ***
MEI 0.064068 0.006434 9.96 < 2e-16 ***
CO2 0.006401 0.002269 2.82 0.0051 **
N2O -0.016021 0.008287 -1.93 0.0542 .
CFC.11 -0.006609 0.001621 -4.08 0.0000595362634 ***
CFC.12 0.003868 0.000981 3.94 0.0001 ***
TSI 0.093116 0.014729 6.32 0.0000000010355 ***
Aerosols -1.540206 0.212616 -7.24 0.0000000000044 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0916 on 276 degrees of freedom
Multiple R-squared: 0.751, Adjusted R-squared: 0.745
F-statistic: 119 on 7 and 276 DF, p-value: <2e-16[E] 誤差與準確性指標
基礎模型的差方和(SST)：\(\sum_i (y_i-\bar{y})^2\)
回歸模型的差方和(SSE)：\(\sum_i (\hat{y_i}-y_i)^2\)
判定係數(\(R^2\))：\(1 - \frac{SSE}{SST}\)
平均誤差(RMSE)：\(\sqrt{SSE/n}\)
§ E1 使用 訓練資料 做預測
pred = predict(m3) # 使用訓練資料做預測
SSE = sum((pred - TR$Temp)^2) # SSE.train
SST = sum((mean(TR$Temp) - TR$Temp)^2) # SST.train
R2 = 1 - SSE/SST # R2.train
RMSE = sqrt(SSE/nrow(TR)) # RMSE.train
c(SSE, SST, R2, RMSE)[1] 2.31351 9.28526 0.75084 0.09026§ E2 使用 測試資料 做預測
pred = predict(m3, TS) # 使用測試資料做預測
SSE = sum((pred - TS$Temp)^2) # SSE.test
SST = sum((mean(TR$Temp) - TS$Temp)^2) # SST.test
R2 = 1 - SSE/SST # R2.test
RMSE = sqrt(SSE/nrow(TS)) # RMSE.test
c(SSE, SST, R2, RMSE)[1] 0.21764 0.58602 0.62861 0.09523注意一下，我們計算 測試資料 的差方和是使用 訓練資料 的平均值做基準
\[ SST_{test} = \sum_{i=1}^n(\bar{y}_{train} - y_{i,test})^2\] 所以 \(R_{test}^2\) 有可能會小於零。
🗿 討論問題：
3.在商務數據分析的角度，我們應該關心的是自變數的：
■
係數估計值
■ \(p\)-value
■ 顯著性
■ 其他，請說明為甚麼？