06:35:00
T
|
rong phạm vi bài viết, mình tập tành
chút trong việc xây dựng mô hình tuyến tính bằng Bayesian Model Average (BMA).
Bài viết được áp dụng trên cơ sở giáo trình “Phân tích dữ liệu với R - Nguyễn Văn Tuấn” (trang 186-188). Theo
đó, phép tính BMA tìm tất cả các mô hình khả dĩ và trình bày kết quả của các mô
hình được xem là “tối ưu” nhất về lâu dài. Tiêu chuẩn tối ưu cũng dựa vào giá
trị Akaike Information Criterion (AIC).
Trước tiên, để sử dụng phép tính BMA,
chúng ta cần tải và cài package BMA.
> library(BMA)
Về dữ liệu:
> sk_kll
d h
dt scanh re than canh la
1 2.16 0.90 0.93 5 0.35
0.12 0.17 0.35
2 5.25 2.10 2.15 14 3.34
1.34 1.10 2.20
3 4.17 1.84 2.01 15 2.25
0.75 0.67 0.95
..........................................
..........................................
Trước tiên, chúng ta cần tạo ra một
ma trận chỉ gồm các biến độc lập. Trong data (sk_kll) chúng ta có các biến d
(đường kính gốc), h (chiều cao cây), dt (đường kính tán) và scanh (số cành dài
trên 50cm)... của cây Keo lá liềm trồng trên đất cát vùng ven biển các tỉnh Bắc
Trung bộ. Biến y (biến tiên đoán), có thể là (sinh khối rễ, sinh khối thân,
sinh khối cành và sinh khối lá). Mục đích của phép tính này, là tiên đoán sinh
khối bộ phận thân cây Keo lá liềm theo một số chỉ tiêu sinh trưởng (d, h, dt,
scanh...).
Bước tiếp theo, chúng ta tạo ra một
data frame mới gồm các biến độc lập (d, h, dt và scanh).
> biendoclap=sk_kll[,c(1:4)]
Lệnh trên có nghĩa, tạo ra data frame
mới gồm 4 biến (d, d, dt và scanh) với tất cả các dòng dữ liệu.
Đối với biến tiên lượng (áp dụng cho
từng biến tiên lượng). Cụ thể:
> re=sk_kll[,c(5)]
> than=sk_kll[,c(6)]
> canh=sk_kll[,c(7)]
> la=sk_kll[,c(8)]
Bây giờ chúng ta đã sẵn sàng phân
tích bằng phép tính BMA. Để phân tích hồi qui tuyến tính, chúng ta sử dụng hàm
bicreg trong package BMA.
> re1=bicreg(biendoclap,re, strict=FALSE, OR=20)
> summary(re1)
Call:
bicreg(x = xvars, y = re,
strict = FALSE, OR = 20)
7
models were selected
Best
5 models (cumulative posterior
probability = 0.7413 ):
p!=0 EV
SD model 1 model 2
model 3
Intercept 100.0
-1.91834 NaN -1.74414
-1.72798 -1.72285
d 61.2 0.61783
NaN 0.96480 1.02960
1.07714
h 38.8 1.22103
NaN . . -0.28197
dt 38.8
0.18728 NaN .
-0.15695 .
scanh 38.8
-0.06831 NaN .
. .
nVar 1 2
2
r2 0.999 0.999
0.999
BIC -19.62465 -18.52604
-18.52604
post prob 0.224 0.129
0.129
model 4 model 5
Intercept -1.75410
-2.61290
d 0.99922 .
h . .
dt . 5.16129
scanh -0.01084 -0.36742
nVar 2 2
r2 0.999 0.999
BIC -18.52604
-18.52604
post prob 0.129
0.129
Kết quả bên trên đưa ra 05 mô hình được
đánh giá là tối ưu nhất cho mô hình tiên đoán sinh khối rễ thân cây Keo lá liềm.
Chúng ta có thể thể hiện kết quả trên
qua biểu đồ dưới đây:
> imageplot.bma(re1)
Tương tự, cho việc tiên đoán cho các
biên tiên lượng (sinh khối thân, cành và lá).
Tiên lượng cho sinh khối thân cây
> than1=bicreg(biendoclap,than, strict=FALSE, OR=20)
> summary(than1)
Call:
bicreg(x = xvars, y = than, strict = FALSE, OR = 20)
8 models were selected
Best 5 models (cumulative posterior probability = 0.625 ):
p!=0 EV SD model 1 model 2 model 3
Intercept 100.0 -0.8168 NaN -0.68873 -0.67005 -0.67005
d 50.0 0.3396 NaN 0.62034 0.79334 0.79334
h 62.5 0.8939 NaN . -1.02619 -1.02619
dt 50.0 -0.3971 NaN -0.57119 . .
scanh 37.5 -0.0464 NaN . . .
nVar 2 2 2
r2 0.999 0.999 0.999
BIC -18.52604 -18.52604 -18.52604
post prob 0.125 0.125 0.125
model 4 model 5
Intercept -0.78377 -1.22190
d 0.50976 .
h . .
dt . 2.63306
scanh -0.03946 -0.22137
nVar 2 2
r2 0.999 0.999
BIC -18.52604 -18.52604
post prob 0.125 0.125
> imageplot.bma(than1)
Tiên lượng sinh khối cành
> canh1=bicreg(biendoclap,canh, strict=FALSE, OR=20)
> summary(canh1)
Call:
bicreg(x = xvars, y = canh, strict = FALSE, OR = 20)
7 models were selected
Best 5 models (cumulative posterior probability = 0.7143 ):
p!=0 EV SD model 1 model 2 model 3
Intercept 100.0 -0.54887 NaN -0.45181 -0.43983 -0.43983
d 57.1 0.27623 NaN 0.44565 0.55664 0.55664
h 57.1 0.38324 NaN . -0.65834 -0.65834
dt 42.9 -0.03839 NaN -0.36644 . .
scanh 42.9 -0.03740 NaN . . .
nVar 2 2 2
r2 0.999 0.999 0.999
BIC -18.52604 -18.52604 -18.52604
post prob 0.143 0.143 0.143
model 4 model 5
Intercept -0.51279 -0.83484
d 0.37471 .
h . .
dt . 1.93548
scanh -0.02532 -0.15903
nVar 2 2
r2 0.999 0.999
BIC -18.52604 -18.52604
post prob 0.143 0.143
> imageplot.bma(canh1)
Tiên lượng cho biến sinh khối lá
> summary(la1)
Call:
bicreg(x = xvars, y = la, strict = FALSE, OR = 20)
7 models were selected
Best 5 models (cumulative posterior probability = 0.7143 ):
p!=0 EV SD model 1 model 2 model 3 model 4
Intercept 100.0 -1.1493 NaN -0.7806 -0.7117 -1.1311 -2.0101
d 42.9 0.6460 NaN 1.4305 2.0686 1.0226 .
h 57.1 2.4123 NaN . -3.7850 . .
dt 57.1 -1.4977 NaN -2.1068 . . 5.2823
scanh 42.9 -0.1348 NaN . . -0.1456 -0.5105
nVar 2 2 2 2
r2 0.999 0.999 0.999 0.999
BIC -18.5260 -18.5260 -18.5260 -18.5260
post prob 0.143 0.143 0.143 0.143
model 5
Intercept -0.9351
d .
h 8.4853
dt -6.8298
scanh .
nVar 2
r2 0.999
BIC -18.5260
post prob 0.143
> imageplot.bma(la1)
Trên đây mình có tập tành chút về áp
dụng R hoặc tương tác với RStudio để xây dựng mô hình tuyến tính bằng Bayesian
Model Average (BMA). Có thể tóm tắt các bước để tìm mô hình “tối ưu” bằng tiêu
chuẩn BMA như sau:
> library(BMA)
> biendoclap=sk_kll[,c(1:4)]
> re=sk_kll[,c(5)]
> than=sk_kll[,c(6)]
> canh=sk_kll[,c(7)]
> la=sk_kll[,c(8)]
> re1=bicreg(biendoclap,re, strict=FALSE, OR=20)
> summary(re1)
> imageplot.bma(re1)
