Hướng dẫn Bootstrap Thống kê & Khoảng tin cậy (CI) trong R

Bootstrap là gì?

Bootstrap là một phương pháp suy luận về tổng thể dựa trên dữ liệu mẫu. Bradley Efron lần đầu giới thiệu trong bài báo này năm 1979. Bootstrap dựa trên việc lấy mẫu có hoàn lại từ dữ liệu mẫu. Kỹ thuật này có thể dùng để ước lượng sai số chuẩn của bất kỳ thống kê nào và để thu được khoảng tin cậy (CI) cho thống kê đó. Bootstrap đặc biệt hữu ích khi CI không có dạng đóng, hoặc có dạng rất phức tạp.

Giả sử ta có một mẫu gồm n phần tử: X = {x₁, x₂, …, x_n} và ta quan tâm đến CI cho một thống kê nào đó T = t(X). Khung Bootstrap khá trực quan. Ta chỉ cần lặp lại R lần sơ đồ sau: Với lần lặp thứ i, lấy mẫu có hoàn lại n phần tử từ mẫu sẵn có (một số phần tử sẽ được chọn nhiều hơn một lần). Gọi mẫu mới này là mẫu bootstrap thứ i, Xi, và tính thống kê mong muốn Ti = t(Xi).

Kết quả, ta sẽ có R giá trị của thống kê: T₁, T₂, …, T_R. Ta gọi chúng là các hiện thực bootstrap của T hoặc phân phối bootstrap của T. Dựa trên đó, ta có thể tính CI cho T. Có vài cách để làm việc này. Lấy các phân vị có vẻ là cách đơn giản nhất.

Bootstrap trong thực tế

Hãy dùng (một lần nữa) bộ dữ liệu iris nổi tiếng. Xem vài dòng đầu:

 head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosaĐược hỗ trợ bởi

Giả sử ta muốn tìm CI cho trung vị Sepal.Length, trung vị Sepal.Width và hệ số tương quan hạng Spearman giữa hai biến này. Ta sẽ dùng gói boot của R và hàm có tên... boot. Để tận dụng sức mạnh của nó, ta phải tạo một hàm tính các thống kê từ dữ liệu được lấy mẫu lại. Hàm này cần có ít nhất hai đối số: một dataset và một vector chứa indices của các phần tử trong dataset được chọn để tạo mẫu bootstrap.

Nếu muốn tính CI cho nhiều thống kê cùng lúc, hàm của ta phải trả về chúng dưới dạng một vector duy nhất.

Với ví dụ của ta, nó có thể trông như sau:

 library(boot)

foo <- function(data, indices){
  dt<-data[indices,]
  c(
    cor(dt[,1], dt[,2], method='s'),
    median(dt[,1]),
    median(dt[,2])
  )
}Được hỗ trợ bởi

foo chọn các phần tử mong muốn (các số được lưu trong indices) từ data và tính hệ số tương quan của hai cột đầu (method='s' dùng để chọn hệ số hạng Spearman, method='p' sẽ cho hệ số Pearson) và trung vị của chúng.

Ta cũng có thể thêm vài đối số phụ, ví dụ cho phép người dùng chọn loại hệ số tương quan:

 foo <- function(data, indices, cor.type){
  dt<-data[indices,]
  c(
    cor(dt[,1], dt[,2], method=cor.type),
    median(dt[,1]),
    median(dt[,2])
  )
}Được hỗ trợ bởi

Để hướng dẫn này tổng quát hơn, tôi sẽ dùng phiên bản sau của foo.

Giờ ta có thể dùng hàm boot. Ta phải cung cấp tên dataset, hàm vừa tạo, số lần lặp (R) và bất kỳ đối số bổ sung nào của hàm (như cor.type). Bên dưới, tôi dùng set.seed để có thể tái lập ví dụ.

 set.seed(12345)
myBootstrap <- boot(iris, foo, R=1000, cor.type='s')Được hỗ trợ bởi

Hàm boot trả về một đối tượng có lớp gọi là... (vâng, bạn đoán đúng!) boot. Nó có hai phần tử thú vị. $t chứa R giá trị của các thống kê được tạo bởi quy trình bootstrap (các hiện thực bootstrap của T):

 head(myBootstrap$t)

##             [,1] [,2] [,3]
## [1,] -0.26405188 5.70    3
## [2,] -0.12973299 5.80    3
## [3,] -0.07972066 5.75    3
## [4,] -0.16122705 6.00    3
## [5,] -0.20664808 6.00    3
## [6,] -0.12221170 5.80    3Được hỗ trợ bởi

$t0 chứa giá trị của các thống kê trong dataset gốc, đầy đủ:

 myBootstrap$t0

## [1] -0.1667777  5.8000000  3.0000000Được hỗ trợ bởi

In đối tượng boot ra console cho thêm một số thông tin:

 myBootstrap

##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = iris, statistic = foo, R = 1000, cor.type = "s")
##
##
## Bootstrap Statistics :
##       original       bias    std. error
## t1* -0.1667777  0.002546391  0.07573983
## t2*  5.8000000 -0.013350000  0.10295571
## t3*  3.0000000  0.007900000  0.02726414Được hỗ trợ bởi

original giống với $t0. bias là hiệu giữa trung bình các hiện thực bootstrap (những giá trị trong $t), gọi là ước lượng bootstrap của T và giá trị trong dataset gốc (giá trị trong $t0).

 colMeans(myBootstrap$t)-myBootstrap$t0

## [1]  0.002546391 -0.013350000  0.007900000Được hỗ trợ bởi

std. error là sai số chuẩn của ước lượng bootstrap, bằng với độ lệch chuẩn của các hiện thực bootstrap.

 apply(myBootstrap$t,2,sd)

## [1] 0.07573983 0.10295571 0.02726414Được hỗ trợ bởi

Các loại CI Bootstrap khác nhau

Trước khi bắt đầu với CI, luôn đáng để xem qua phân phối của các hiện thực bootstrap. Ta có thể dùng hàm plot, với index cho biết ta muốn xem thống kê nào trong số những thống kê tính trong foo. Ở đây index=1 là hệ số tương quan Spearman giữa chiều dài và chiều rộng đài hoa, index=2 là trung vị chiều dài đài hoa, và index=3 là trung vị chiều rộng đài hoa.

plot(myBootstrap, index=1)
Được hỗ trợ bởi

Phân phối của các hệ số tương quan bootstrap có vẻ khá giống chuẩn. Hãy tìm CI cho nó. Ta có thể dùng boot.ci. Mặc định là CI 95%, nhưng có thể thay đổi bằng tham số conf.

 boot.ci(myBootstrap, index=1)

## Warning in boot.ci(myBootstrap, index = 1): bootstrap variances needed for
## studentized intervals

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = myBootstrap, index = 1)
##
## Intervals :
## Level      Normal              Basic         
## 95%   (-0.3178, -0.0209 )   (-0.3212, -0.0329 )  
##
## Level     Percentile            BCa          
## 95%   (-0.3007, -0.0124 )   (-0.3005, -0.0075 )  
## Calculations and Intervals on Original ScaleĐược hỗ trợ bởi

boot.ci cung cấp 5 loại CI bootstrap. Một trong số đó, khoảng studentized, là đặc biệt. Nó cần ước lượng phương sai bootstrap. Ta không cung cấp, nên R in cảnh báo: bootstrap variances needed for studentized intervals. Có thể thu được ước lượng phương sai bằng bootstrap cấp hai hoặc (dễ hơn) bằng kỹ thuật jackknife. Điều này hơi vượt ra ngoài phạm vi hướng dẫn này, nên hãy tập trung vào bốn loại CI bootstrap còn lại.

Nếu không muốn xem tất cả, ta có thể chọn những loại phù hợp trong đối số type. Các giá trị có thể là norm, basic, stud, perc, bca hoặc một vector gồm các giá trị này.

 boot.ci(myBootstrap, index=1, type=c('basic','perc'))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = myBootstrap, type = c("basic", "perc"), index = 1)
##
## Intervals :
## Level      Basic              Percentile     
## 95%   (-0.3212, -0.0329 )   (-0.3007, -0.0124 )  
## Calculations and Intervals on Original ScaleĐược hỗ trợ bởi

Hàm boot.ci tạo đối tượng có lớp... (đúng vậy, bạn đoán rồi!) bootci. Các phần tử của nó được đặt tên giống với các loại CI dùng trong đối số type. $norm là một vector 3 phần tử chứa mức tin cậy và hai cận của CI.

 boot.ci(myBootstrap, index=1, type='norm')$norm

##      conf                       
## [1,] 0.95 -0.3177714 -0.02087672Được hỗ trợ bởi

$basic, $stud, $perc và $bca là các vector 5 phần tử cũng chứa các phân vị dùng để tính CI (ta sẽ quay lại điều này sau):

 boot.ci(myBootstrap, index=1, type='basic')$basic

##      conf                                    
## [1,] 0.95 975.98 25.03 -0.3211981 -0.03285178Được hỗ trợ bởi

Một chút ký hiệu (xin lỗi!)

Để hiểu cảm giác về các loại CI khác nhau, ta cần giới thiệu một số ký hiệu. Gọi:

t^⋆ là ước lượng bootstrap (trung bình các hiện thực bootstrap),
t₀ là giá trị thống kê trong dataset gốc,
se^⋆ là sai số chuẩn của ước lượng bootstrap,
b là độ chệch của ước lượng bootstrap b = t^⋆ − t₀
α là mức tin cậy, thường α = 0.95,
z_α là bách phân vị $1-\frac \alpha 2$ của phân phối chuẩn chuẩn hóa,
θ_α là một α-phân vị của phân phối các hiện thực bootstrap.

CI theo phân vị

Với ký hiệu trên, CI theo phân vị là:

(θ_{(1 − α)/2}, θ_{1 − (1 − α)/2})

nên về cơ bản chỉ là lấy các phân vị tương ứng. Không hơn.

CI chuẩn

CI Wald điển hình sẽ là:

t₀ ± z_α ⋅ se^⋆

nhưng trong trường hợp bootstrap, ta nên hiệu chỉnh độ chệch. Vậy nó trở thành:

$$ t_0 - b \pm z_\alpha \cdot se^\star \\ 2t_0 - t^\star \pm z_\alpha \cdot se^\star$$

CI cơ bản

CI theo phân vị thường không được khuyến nghị vì hoạt động kém với các phân phối có đuôi “kỳ lạ”. CI cơ bản (còn gọi là pivotal hoặc CI empirical) vững hơn nhiều trước trường hợp này. Lý do là tính hiệu giữa mỗi lần lặp bootstrap và t₀ rồi dùng các phân vị của phân phối các hiệu đó. Chi tiết đầy đủ có thể xem, ví dụ trong sách của L. Wasserman All of Statistics

Công thức cuối cùng cho CI cơ bản là:

(2t₀ − θ_{1 − (1 − α)/2}, 2t₀ − θ_{(1 − α)/2})

BC_α CI

BC_α xuất phát từ bias-corrected, accelerated (hiệu chỉnh độ chệch, tăng tốc). Công thức không quá phức tạp nhưng khá khó trực giác, nên tôi sẽ lược bỏ. Xem bài viết của Thomas J. DiCiccio và Bradley Efron nếu bạn quan tâm chi tiết.

Yếu tố “tăng tốc” trong tên phương pháp yêu cầu dùng các phân vị cụ thể của các hiện thực bootstrap. Đôi khi có thể xảy ra việc các phân vị đó là các giá trị cực đoan, có thể là ngoại lệ. BC_α có thể không ổn định trong các trường hợp như vậy.

Hãy xem CI BC_α cho trung vị chiều rộng cánh hoa. Trong dataset gốc, trung vị này đúng bằng 3.

 boot.ci(myBootstrap, index=3)

## Warning in boot.ci(myBootstrap, index = 3): bootstrap variances needed for
## studentized intervals

## Warning in norm.inter(t, adj.alpha): extreme order statistics used as
## endpoints

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = myBootstrap, index = 3)
##
## Intervals :
## Level      Normal              Basic         
## 95%   ( 2.939,  3.046 )   ( 2.900,  3.000 )  
##
## Level     Percentile            BCa          
## 95%   ( 3.0,  3.1 )   ( 2.9,  2.9 )  
## Calculations and Intervals on Original Scale
## Warning : BCa Intervals used Extreme Quantiles
## Some BCa intervals may be unstableĐược hỗ trợ bởi

Ta nhận được CI BC_α (2.9, 2.9). Điều đó khá kỳ lạ, nhưng may là R đã cảnh báo rằng extreme order statistics used as endpoints. Hãy xem chính xác điều gì xảy ra ở đây:

 plot(myBootstrap, index=3)Được hỗ trợ bởi

Phân phối của các hiện thực bootstrap là bất thường. Phần lớn (hơn 90%) đều bằng 3.

 table(myBootstrap$t[,3])

##
##  2.9 2.95    3 3.05  3.1 3.15  3.2
##    1    1  908   24   63    1    2Được hỗ trợ bởi

Trong tình huống như vậy, “tăng tốc” trong phương pháp BC_α dễ “nhảy” đến các giá trị cực đoan. Ở đây, CI theo phân vị, vốn là cơ sở cho CI BC_α, là (3, 3.1). Giờ khi ta mở rộng nó “về bên trái” ta phải dùng phân vị 0.002. Phần mở rộng “về bên phải” cũng đi đến cực trị: phân vị 0.997.

Tái lập kết quả

Đôi khi ta cần tái tạo các lần lặp bootstrap. Nếu có thể dùng R cho việc này thì không thành vấn đề. Hàm set.seed giải quyết điều đó. Tái tạo các lần lặp bootstrap trong phần mềm khác sẽ phức tạp hơn nhiều. Hơn nữa, với R lớn, việc tính lại trong R cũng có thể không khả thi (ví dụ do thiếu thời gian).

Ta có thể xử lý vấn đề này bằng cách lưu các chỉ số của các phần tử trong dataset gốc đã tạo nên mỗi mẫu bootstrap. Đây là điều mà hàm boot.array (với đối số indices=T) thực hiện.

 tableOfIndices<-boot.array(myBootstrap, indices=T)Được hỗ trợ bởi

Mỗi hàng là một mẫu bootstrap. Ví dụ, mẫu đầu tiên của ta gồm các phần tử sau:

 tableOfIndices[1,]

##   [1] 109  12 143  65  41  28 105  62  23 102  37  34  55  53  38   3  11
##  [18] 146  31 122  85 141 118 116 117  78 100  13  98  49  59  88  24  56
##  [35] 136   4  59  67 126 118 104 101  70  13  12  19  21 149  73 133  67
##  [52]  86   6  88 131 105 121 145 121  83  70  68  71 111  76  73 122 116
##  [69]  85 144 114 139  89 124  64  27  81  78  58  39  17  50  37 117  76
##  [86]  97 114  64  17  93 141  65 124  80 137  54  37  57  70 128  10  55
## [103]  13  53  59  45 116  14  77 118 108 138  50  78  49 104  49   1  85
## [120]  28  43  82  74  64  33  55  32  59  62 112   8  11  96  30  14  30
## [137]  38  85  66  85  97 107   4  18  76  35  31 133  27  69Được hỗ trợ bởi

Đặt indices=F (mặc định), ta sẽ có câu trả lời cho câu hỏi “Mỗi phần tử của dataset gốc xuất hiện bao nhiêu lần trong mỗi mẫu bootstrap?”. Ví dụ, trong mẫu đầu tiên: phần tử thứ 1 của dataset xuất hiện một lần, phần tử thứ 2 không xuất hiện, phần tử thứ 3 xuất hiện một lần, phần tử thứ 4 hai lần, v.v.

 tableOfAppearances<-boot.array(myBootstrap)
tableOfAppearances[1,]

##   [1] 1 0 1 2 0 1 0 1 0 1 2 2 3 2 0 0 2 1 1 0 1 0 1 1 0 0 2 2 0 2 2 1 1 1 1
##  [36] 0 3 2 1 0 1 0 1 0 1 0 0 0 3 2 0 0 2 1 3 1 1 1 4 0 0 2 0 3 2 1 2 1 1 3
##  [71] 1 0 2 1 0 3 1 3 0 1 1 1 1 0 5 1 0 2 1 0 0 0 1 0 0 1 2 1 0 1 1 1 0 2 2
## [106] 0 1 1 1 0 1 1 0 2 0 3 2 3 0 0 2 2 0 2 0 1 0 1 0 0 1 0 2 0 0 1 1 1 1 0
## [141] 2 0 1 1 1 1 0 0 1 0Được hỗ trợ bởi

Một bảng như vậy cho phép ta tái tạo các hiện thực bootstrap ngoài R. Hoặc ngay trong R khi ta không muốn dùng set.seed và làm lại toàn bộ tính toán.

 onceAgain<-apply(tableOfIndices, 1, foo, data=iris, cor.type='s')Được hỗ trợ bởi

Hãy kiểm tra xem kết quả có giống nhau không:

 head(t(onceAgain))

##             [,1] [,2] [,3]
## [1,] -0.26405188 5.70    3
## [2,] -0.12973299 5.80    3
## [3,] -0.07972066 5.75    3
## [4,] -0.16122705 6.00    3
## [5,] -0.20664808 6.00    3
## [6,] -0.12221170 5.80    3

head(myBootstrap$t)

##             [,1] [,2] [,3]
## [1,] -0.26405188 5.70    3
## [2,] -0.12973299 5.80    3
## [3,] -0.07972066 5.75    3
## [4,] -0.16122705 6.00    3
## [5,] -0.20664808 6.00    3
## [6,] -0.12221170 5.80    3

all(t(onceAgain)==myBootstrap$t)

## [1] TRUEĐược hỗ trợ bởi

Vâng, giống hệt!

Nếu bạn muốn tìm hiểu thêm về Machine Learning trong R, hãy học khóa Machine Learning Toolbox của DataCamp và xem hướng dẫn Machine Learning trong R cho người mới bắt đầu.

Chủ đề

R Khoa học Dữ liệu

Chủ đề

Công nghệ

Hướng dẫn Bootstrap trong R

Khám phá cùng AI

Bootstrap là gì?

Bootstrap trong thực tế

Các loại CI Bootstrap khác nhau

Một chút ký hiệu (xin lỗi!)

CI theo phân vị

CI chuẩn

CI cơ bản

BC_α CI

Tái lập kết quả

Chủ đề

Tìm hiểu thêm về R và Machine Learning

Machine Learning với caret trong R

Machine Learning cho Phân tích Marketing bằng R

Học không giám sát với R

Liên quan

Claude Opus 4.6: Tính năng, Điểm chuẩn, Bài kiểm tra thực hành và hơn thế nữa

Chủ đề

Công nghệ

Tìm hiểu thêm về R và Machine Learning

Machine Learning với caret trong R

Machine Learning cho Phân tích Marketing bằng R

Học không giám sát với R

Bootstrap là gì?

Bootstrap trong thực tế

Các loại CI Bootstrap khác nhau

Một chút ký hiệu (xin lỗi!)

CI theo phân vị

CI chuẩn

CI cơ bản

BCα CI

Tái lập kết quả

Chủ đề

Tìm hiểu thêm về R và Machine Learning

Liên quan

Tìm hiểu thêm về R và Machine Learning

BC_α CI