Transcript Phân tích hồi quy tuyến tính đa biến

Phân tích hồi quy tuyến tính đa
biến sử dụng phần mềm R
Tuan V. Nguyen
Garvan Institute of Medical Research,
Sydney, Australia
Số liệu
Tác động của nhiệt độ (x1), thời gian (x2) lên CO2 (y)
Id
y
X1
X2
X3
1
36.98
5.1
400
51.37
2
13.74
26.4
400
3
10.08
23.8
4
8.53
5
X4
X5
X6
X7
4.24
1484.83
2227.25
2.06
72.33
30.87
289.94
434.90
1.33
400
71.44
33.01
320.79
481.19
0.97
46.4
400
79.15
44.61
164.76
247.14
0.62
36.42
7.0
450
80.47
33.84
1097.26
1645.89
0.22
6
26.59
12.6
450
89.90
41.26
605.06
907.59
0.76
7
19.07
18.9
450
91.48
41.88
405.37
608.05
1.71
8
5.96
30.2
450
98.60
70.79
253.70
380.55
3.93
9
15.52
53.8
450
98.05
66.82
142.27
213.40
1.97
…
…
…
…
…
…
…
…
…
25
11.19
11.5
450
77.88
25.20
663.09
994.63
1.61
26
75.62
5.2
470
75.50
8.66
1464.11
2196.17
4.78
27
36.03
10.6
470
83.15
22.39
720.07
1080.11
5.88
Chú thích: y = sản lượng CO2; X1 = thời gian (phút); X2 = nhiệt độ (C);
X3 = phần trăm hòa tan; X4 = lượng dầu (g/100g); X5 = lượng than đá; X6
= tổng số lượng hòa tan; X7 = số hydrogen tiêu thụ.
Nhập số liệu vào R
y <- c(36.98,13.74,10.08, 8.53,36.42,26.59,19.07, 5.96,15.52,56.61,
26.72,20.80, 6.99,45.93,43.09,15.79,21.60,35.19,26.14, 8.60,
11.63, 9.59, 4.42,38.89,11.19,75.62,36.03)
x1 <- c(5.1,26.4,23.8,46.4, 7.0,12.6,18.9,30.2,53.8,5.6,15.1,20.3,48.4,
5.8,11.2,27.9,5.1,11.7,16.7,24.8,24.9,39.5,29.0, 5.5, 11.5,
5.2,10.6)
x2 <- c(400,400, 400, 400, 450, 450, 450, 450, 450, 400, 400, 400,
400, 425, 425, 425, 450, 450, 450, 450, 450, 450, 450, 460,
450, 470, 470)
x3 <- c(51.37,72.33,71.44,79.15,80.47,89.90,91.48,98.60,98.05,55.69,
66.29,58.94,74.74,63.71,67.14,77.65,67.22,81.48,83.88,89.38,
79.77,87.93,79.50,72.73,77.88,75.50,83.15)
x4 <- c(4.24,30.87,33.01,44.61,33.84,41.26,41.88,70.79,66.82,
8.92,17.98,17.79,33.94,11.95,14.73,34.49,14.48,29.69,26.33,
37.98,25.66,22.36,31.52,17.86,25.20, 8.66,22.39)
x5 <- c(1484.83, 289.94, 320.79, 164.76, 1097.26, 605.06, 405.37,
253.70, 142.27,1362.24, 507.65, 377.60, 158.05, 130.66,
682.59, 274.20, 1496.51, 652.43, 458.42, 312.25, 307.08,
193.61, 155.96,1392.08, 663.09,1464.11, 720.07)
x6 <- c(2227.25, 434.90, 481.19, 247.14,1645.89, 907.59, 608.05,
380.55, 213.40,2043.36, 761.48, 566.40, 237.08,1961.49,1023.89,
411.30,2244.77, 978.64, 687.62, 468.38, 460.62, 290.42,
233.95,2088.12, 994.63,2196.17,1080.11)
x7 <- c(2.06,1.33,0.97,0.62,0.22,0.76,1.71,3.93,1.97,5.08,0.60,0.90,
0.63,2.04,1.57,2.38,0.32,0.44,8.82,0.02,1.72,1.88,1.43,
1.35,1.61,4.78,5.88)
REGdata <- data.frame(y, x1,x2,x3,x4,x5,x6,x7)
Muốn lưu giữ file, cần phải dùng lệnh
Write.table(REGdata,”địa chỉ cần lưu giữ file/data.csv”)
Hoặc gọi file từ R
Nhập số liệu từ bảng biểu trình bày trên vào file
xcel. Sau đó save as dưới dạng file csv. Cách
này có thể lưu giữ số liệu gốc.
Rồi gọi file số liệu trưc tiếp vào R:
REGdata <- read.csv(“c:/đường dẫn đến
nơi lưu file/ data.csv, header=T)
Kết quả phân tích
reg <- lm(y ~ x1+x2+x3+x4+x5+x6+x7, data=REGdata)
summary(reg)
Residuals:
Min
1Q Median
3Q
Max
-20.035 -4.681 -1.144
4.072 21.214
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 53.937016 57.428952
0.939
0.3594
x1
-0.127653
0.281498 -0.453
0.6553
x2
-0.229179
0.232643 -0.985
0.3370
x3
0.824853
0.765271
1.078
0.2946
x4
-0.438222
0.358551 -1.222
0.2366
x5
-0.001937
0.009654 -0.201
0.8431
x6
0.019886
0.008088
2.459
0.0237 *
x7
1.993486
1.089701
1.829
0.0831 .
--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.61 on 19 degrees of freedom
Multiple R-Squared: 0.728,
Adjusted R-squared: 0.6278
F-statistic: 7.264 on 7 and 19 DF, p-value: 0.0002674
Tương quan chéo
30
50
50
70
90
200
1000
0
4
8
10 40 70
10
10 30 50
y
400 440
x1
50 70 90
x2
10 40 70
x3
200 1000
x4
2000
x5
8
500
x6
0
4
x7
10
40
70
400
440
10
40
70
500
2000
Hồi quy loại từng bước (stepwise)
reg <- lm(y ~ ., data=REGdata)
step(reg, direction=”both”)
Hồi quy loại từng bước
Start: AIC= 134.07
y ~ x1 + x2 + x3 +
Df Sum of Sq
- x5
1
4.54
- x1
1
23.17
- x2
1
109.34
- x3
1
130.90
<none>
- x4
1
168.31
- x7
1
377.09
- x6
1
681.09
x4 + x5 + x6 + x7
RSS
AIC
2145.37 132.13
2164.00 132.36
2250.18 133.42
2271.74 133.68
2140.83 134.07
2309.14 134.12
2517.92 136.45
2821.92 139.53
Step 2: AIC= 130.42
y ~ x2 + x3 + x4 + x6 + x7
Df Sum of Sq
RSS
- x2
1
96.8 2264.9
- x3
1
122.0 2290.0
<none>
2168.1
- x4
1
187.4 2355.5
+ x1
1
22.7 2145.4
+ x5
1
4.1 2164.0
- x7
1
385.0 2553.1
- x6
1
1526.2 3694.3
AIC
129.6
129.9
130.4
130.7
132.1
132.4
132.8
142.8
Step 1: AIC= 132.13
y ~ x1 + x2 + x3 + x4 + x6
Df Sum of Sq
RSS
- x1
1
22.7 2168.1
- x2
1
113.8 2259.1
- x3
1
133.5 2278.9
<none>
2145.4
- x4
1
170.8 2316.2
+ x5
1
4.5 2140.8
- x7
1
375.7 2521.1
- x6
1
1058.5 3203.8
+ x7
AIC
130.4
131.5
131.8
132.1
132.2
134.1
134.5
141.0
Step 3: AIC= 129.59
y ~ x3 + x4 + x6 + x7
Df Sum of Sq
RSS
- x3
1
25.4 2290.3
- x4
1
90.9 2355.8
<none>
2264.9
+ x2
1
96.8 2168.1
+ x5
1
8.3 2256.5
+ x1
1
5.7 2259.1
- x7
1
384.9 2649.7
- x6
1
2015.6 4280.5
AIC
127.9
128.7
129.6
130.4
131.5
131.5
131.8
144.8
Hồi quy loại từng bước
Step 4: AIC= 127.9
y ~ x4 + x6 + x7
Df Sum of Sq
- x4
1
73.5
<none>
+ x3
1
25.4
+ x1
1
11.3
+ x5
1
6.3
+ x2
1
0.3
- x7
1
486.6
- x6
1
1993.8
RSS
2363.8
2290.3
2264.9
2279.0
2284.0
2290.0
2776.9
4284.1
AIC
126.7
127.9
129.6
129.8
129.8
129.9
131.1
142.8
Call:
lm(formula = y ~ x6 + x7, data =
REGdata)
Coefficients:
(Intercept)
x6
x7
2.52646
0.01852
2.18575
Step 5: AIC= 126.75
y ~ x6 + x7
Df Sum of Sq
RSS
<none>
2363.8
+ x4
1
73.5 2290.3
+ x1
1
33.4 2330.4
+ x3
1
8.1 2355.8
+ x5
1
7.7 2356.1
+ x2
1
7.3 2356.6
- x7
1
497.3 2861.2
- x6
1
4477.0 6840.8
AIC
126.7
127.9
128.4
128.7
128.7
128.7
129.9
153.4
Phương pháp lựa chọn mô hình tối ưu
bằng Bayesian Model Average (BMA)
library(BMA)
xvars <- REGdata[,-1]
co2 <- REGdata[,1]
bma <- bicreg(xvars, co2, strict=FALSE, OR=20)
summary(bma)
Bayesian Model Average (BMA) analysis
16 models were selected
Best 5 models (cumulative posterior probability = 0.6599 ):
p!=0
EV
SD
model 1
model 2
model 3
Intercept 100.0
5.75672 14.6244
2.5264
6.1441
8.6120
x1
12.4 -0.01807
0.1008
.
.
.
x2
10.4 -0.00075
0.0282
.
.
.
x3
10.7
0.00011
0.0791
.
.
.
x4
20.2 -0.03059
0.1020
.
.
-0.1419
x5
10.5 -0.00023
0.0030
.
.
.
x6
100.0
0.01815
0.0040
0.0185
0.0193
0.0164
x7
73.7
1.60766
1.2821
2.1857
.
2.1628
model 4
7.5936
-0.1393
.
.
.
.
0.0162
2.1233
model 5
7.3537
.
.
-0.0572
.
.
0.0179
2.2382
nVar
r2
BIC
post prob
3
0.704
-22.9721
0.072
3
0.701
-22.6801
0.063
2
0.700
-25.8832
0.311
1
0.636
-24.0238
0.123
3
0.709
-23.4412
0.092
Bayesian Model Average (BMA) analysis
imageplot.bma(bma)
Models selected by BMA
x1
x2
x3
x4
x5
x6
x7
1
2
3
Model #
4
5
6
7
8
10
13
Lời Cảm tạ
• Chúng tôi xin chân thành cám ơn
Công ty Dược phẩm Bridge
Healthcare, Australia đã tài trợ cho
chuyến đi.