Xây dựng mô hình tuyến tính bằng Bayesian Model Average (BMA)

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)