Tutorial Statistik Bootstrap & Interval Kepercayaan (CI) di R

Apa itu Bootstrap?

Bootstrap adalah metode inferensi tentang suatu populasi menggunakan data sampel. Bradley Efron pertama kali memperkenalkannya dalam makalah ini pada tahun 1979. Bootstrap bergantung pada pengambilan sampel dengan pengembalian dari data sampel. Teknik ini dapat digunakan untuk memperkirakan galat baku dari statistik apa pun dan untuk memperoleh interval kepercayaan (CI) untuk statistik tersebut. Bootstrap sangat berguna ketika CI tidak memiliki bentuk tertutup, atau memiliki bentuk yang sangat rumit.

Misalkan kita memiliki sampel berukuran n elemen: X = {x₁, x₂, …, x_n} dan kita tertarik pada CI untuk suatu statistik T = t(X). Kerangka bootstrap cukup lurus. Kita cukup mengulangi sebanyak R kali skema berikut: Untuk pengulangan ke-i, ambil sampel dengan pengembalian sebanyak n elemen dari sampel yang tersedia (beberapa elemen akan terambil lebih dari sekali). Sebut sampel baru ini sebagai sampel bootstrap ke-i, Xi, dan hitung statistik yang diinginkan Ti = t(Xi).

Hasilnya, kita akan mendapatkan R nilai statistik kita: T₁, T₂, …, T_R. Kita menyebutnya realisasi bootstrap dari T atau distribusi bootstrap dari T. Berdasarkan ini, kita dapat menghitung CI untuk T. Ada beberapa cara untuk melakukannya. Mengambil persentil tampaknya yang paling mudah.

Bootstrap dalam Praktik

Mari gunakan (sekali lagi) dataset iris yang terkenal. Lihat beberapa baris pertamanya:

 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  setosaDidukung oleh

Misalkan kita ingin mencari CI untuk median Sepal.Length, median Sepal.Width dan koefisien korelasi peringkat Spearman antara keduanya. Kita akan menggunakan paket boot di R dan sebuah fungsi bernama... boot. Untuk memanfaatkannya kita harus membuat fungsi yang menghitung statistik kita dari data hasil resampling. Fungsi tersebut harus memiliki setidaknya dua argumen: sebuah dataset dan sebuah vektor berisi indices elemen dari dataset yang terpilih untuk membentuk sampel bootstrap.

Jika kita ingin menghitung CI untuk lebih dari satu statistik sekaligus, fungsi kita harus mengembalikannya sebagai satu vektor.

Untuk contoh kita, fungsinya bisa seperti ini:

 library(boot)

foo <- function(data, indices){
  dt<-data[indices,]
  c(
    cor(dt[,1], dt[,2], method='s'),
    median(dt[,1]),
    median(dt[,2])
  )
}Didukung oleh

foo memilih elemen yang diinginkan (nomornya disimpan di indices) dari data dan menghitung koefisien korelasi dari dua kolom pertama (method='s' digunakan untuk memilih koefisien peringkat Spearman, method='p' akan menghasilkan koefisien Pearson) serta median keduanya.

Kita juga bisa menambahkan argumen ekstra, misalnya biarkan pengguna memilih jenis koefisien korelasi:

 foo <- function(data, indices, cor.type){
  dt<-data[indices,]
  c(
    cor(dt[,1], dt[,2], method=cor.type),
    median(dt[,1]),
    median(dt[,2])
  )
}Didukung oleh

Untuk membuat tutorial ini lebih umum, saya akan menggunakan versi foo yang terakhir.

Sekarang, kita dapat menggunakan fungsi boot. Kita harus memberinya nama dataset, fungsi yang baru saja kita buat, jumlah pengulangan (R) dan argumen tambahan apa pun dari fungsi kita (seperti cor.type). Di bawah ini, saya menggunakan set.seed agar contoh ini dapat direplikasi.

 set.seed(12345)
myBootstrap <- boot(iris, foo, R=1000, cor.type='s')Didukung oleh

Fungsi boot mengembalikan sebuah objek dengan kelas bernama... (ya, benar!) boot. Objek ini memiliki dua elemen yang menarik. $t berisi R nilai statistik kita yang dihasilkan oleh prosedur bootstrap (realisasi bootstrap dari 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    3Didukung oleh

$t0 berisi nilai statistik kita pada dataset asli, lengkap:

 myBootstrap$t0

## [1] -0.1667777  5.8000000  3.0000000Didukung oleh

Mencetak objek boot ke konsol memberikan informasi lebih lanjut:

 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.02726414Didukung oleh

original sama dengan $t0. bias adalah selisih antara rerata realisasi bootstrap (yang ada di $t), yang disebut estimasi bootstrap dari T dan nilai pada dataset asli (yang ada di $t0).

 colMeans(myBootstrap$t)-myBootstrap$t0

## [1]  0.002546391 -0.013350000  0.007900000Didukung oleh

std. error adalah galat baku dari estimasi bootstrap, yang sama dengan simpangan baku dari realisasi bootstrap.

 apply(myBootstrap$t,2,sd)

## [1] 0.07573983 0.10295571 0.02726414Didukung oleh

Berbagai Jenis CI Bootstrap

Sebelum mulai dengan CI, selalu ada baiknya melihat distribusi realisasi bootstrap. Kita bisa menggunakan fungsi plot, dengan index yang menunjukkan statistik mana yang dihitung dalam foo yang ingin kita lihat. Di sini index=1 adalah koefisien korelasi Spearman antara panjang dan lebar sepal, index=2 adalah median panjang sepal, dan index=3 adalah median lebar sepal.

plot(myBootstrap, index=1)
Didukung oleh

Distribusi koefisien korelasi bootstrap terlihat cukup mendekati normal. Mari cari CI-nya. Kita bisa menggunakan boot.ci. Secara bawaan ia menggunakan CI 95 persen, tetapi dapat diubah dengan parameter 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 ScaleDidukung oleh

boot.ci menyediakan 5 jenis CI bootstrap. Salah satunya, interval studentized, bersifat unik. Ia memerlukan estimasi varians bootstrap. Kita tidak menyediakannya, sehingga R mencetak peringatan: bootstrap variances needed for studentized intervals. Estimasi varians dapat diperoleh dengan bootstrap tingkat kedua atau (lebih mudah) dengan teknik jackknife. Ini agak di luar cakupan tutorial ini, jadi mari fokus pada empat jenis CI bootstrap lainnya.

Jika kita tidak ingin melihat semuanya, kita dapat memilih yang relevan di argumen type. Nilai yang mungkin adalah norm, basic, stud, perc, bca atau vektor dari nilai-nilai ini.

 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 ScaleDidukung oleh

Fungsi boot.ci membuat objek dengan kelas... (ya, Anda menebaknya!) bootci. Elemennya dinamai sama seperti jenis CI yang digunakan di argumen type. $norm adalah vektor 3 elemen yang berisi tingkat kepercayaan dan batas CI.

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

##      conf                       
## [1,] 0.95 -0.3177714 -0.02087672Didukung oleh

$basic, $stud, $perc dan $bca adalah vektor 5 elemen yang juga berisi persentil yang digunakan untuk menghitung CI (kita akan kembali ke ini nanti):

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

##      conf                                    
## [1,] 0.95 975.98 25.03 -0.3211981 -0.03285178Didukung oleh

Sedikit notasi (Maaf!)

Untuk memahami seperti apa berbagai jenis CI, kita perlu memperkenalkan beberapa notasi. Jadi, misalkan:

t^⋆ adalah estimasi bootstrap (rerata realisasi bootstrap),
t₀ adalah nilai statistik kita pada dataset asli,
se^⋆ adalah galat baku dari estimasi bootstrap,
b adalah bias dari estimasi bootstrap b = t^⋆ − t₀
α adalah tingkat kepercayaan, umumnya α = 0.95,
z_α adalah kuantil $1-\frac \alpha 2$ dari distribusi normal baku,
θ_α adalah α-persentil dari distribusi realisasi bootstrap.

CI Persentil

Dengan notasi di atas, CI persentil adalah:

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

jadi ini hanya mengambil persentil yang relevan. Tidak lebih.

CI Normal

CI Wald yang khas adalah:

t₀ ± z_α ⋅ se^⋆

tetapi dalam kasus bootstrap, kita perlu mengoreksinya terhadap bias. Jadi menjadi:

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

CI Basic

CI persentil umumnya tidak direkomendasikan karena berkinerja buruk pada distribusi dengan ekor yang aneh. CI basic (juga disebut CI pivotal atau empiris) jauh lebih tangguh terhadap hal ini. Rasionalnya adalah menghitung selisih antara setiap replikasi bootstrap dan t₀ dan menggunakan persentil dari distribusinya. Rincian lengkapnya dapat ditemukan, misalnya, di buku L. Wasserman All of Statistics

Rumus akhir untuk CI basic adalah:

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

BC_α CI

BC_α berasal dari bias-corrected, accelerated. Rumusnya tidak terlalu rumit tetapi agak tidak intuitif, jadi akan saya lewati. Lihat artikel oleh Thomas J. DiCiccio dan Bradley Efron jika Anda tertarik pada detailnya.

Percepatan yang disebutkan dalam nama metodenya memerlukan penggunaan persentil tertentu dari realisasi bootstrap. Terkadang bisa terjadi bahwa ini adalah persentil ekstrem, yang berpotensi outlier. BC_α dapat menjadi tidak stabil dalam kasus seperti itu.

Mari lihat CI BC_α untuk median lebar petal. Pada dataset asli, median ini tepat 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 unstableDidukung oleh

Kita mendapatkan CI BC_α (2,9, 2,9). Itu aneh, tetapi untungnya R sudah memperingatkan bahwa extreme order statistics used as endpoints. Mari lihat apa yang sebenarnya terjadi di sini:

 plot(myBootstrap, index=3)Didukung oleh

Distribusi realisasi bootstrap tidak biasa. Sebagian besar di antaranya (lebih dari 90%) bernilai 3.

 table(myBootstrap$t[,3])

##
##  2.9 2.95    3 3.05  3.1 3.15  3.2
##    1    1  908   24   63    1    2Didukung oleh

Dalam situasi seperti ini, "percepatan" dalam metode BC_α mudah "melompat" ke nilai ekstrem. Di sini, CI persentil yang menjadi dasar CI BC_α adalah (3, 3,1). Sekarang, ketika kita memperluasnya "ke kiri" kita harus menggunakan persentil 0,002. Perluasan "ke kanan" juga menuju ekstrem: persentil 0,997.

Replikasi Hasil

Terkadang, kita perlu membuat ulang replikasi bootstrap. Jika kita bisa menggunakan R untuk itu, ini bukan masalah. Fungsi set.seed menyelesaikan masalah tersebut. Membuat ulang replikasi bootstrap di perangkat lunak lain akan jauh lebih rumit. Selain itu, untuk R yang besar, perhitungan ulang di R juga bisa bukan pilihan (misalnya karena keterbatasan waktu).

Kita bisa mengatasi masalah ini dengan menyimpan indeks elemen dari dataset asli yang membentuk setiap sampel bootstrap. Inilah yang dilakukan fungsi boot.array (dengan argumen indices=T).

 tableOfIndices<-boot.array(myBootstrap, indices=T)Didukung oleh

Setiap baris adalah satu sampel bootstrap. Misalnya, sampel pertama kita memuat elemen-elemen berikut:

 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  69Didukung oleh

Dengan menyetel indices=F (bawaan), kita akan mendapatkan jawaban atas pertanyaan "Berapa kali setiap elemen dari dataset asli muncul di setiap sampel bootstrap?". Misalnya, dalam sampel pertama: elemen ke-1 dari dataset muncul sekali, elemen ke-2 tidak muncul sama sekali, elemen ke-3 muncul sekali, elemen ke-4 dua kali, dan seterusnya.

 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 0Didukung oleh

Tabel seperti ini memungkinkan kita membuat ulang realisasi bootstrap di luar R. Atau di R sendiri ketika kita tidak ingin menggunakan set.seed dan melakukan semua perhitungan sekali lagi.

 onceAgain<-apply(tableOfIndices, 1, foo, data=iris, cor.type='s')Didukung oleh

Mari periksa apakah hasilnya sama:

 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] TRUEDidukung oleh

Ya, sama!

Jika Anda ingin mempelajari lebih lanjut tentang Machine Learning di R, ikuti kursus Machine Learning Toolbox dari DataCamp dan baca tutorial Machine Learning di R untuk pemula.

Topik

R Data Science

Teknologi

Topik

Tutorial Bootstrap di R

Jelajahi bersama AI

Apa itu Bootstrap?

Bootstrap dalam Praktik

Berbagai Jenis CI Bootstrap

Sedikit notasi (Maaf!)

CI Persentil

CI Normal

CI Basic

BC_α CI

Replikasi Hasil

Topik

Pelajari lebih lanjut tentang R dan Machine Learning

Machine Learning dengan caret di R

Machine Learning untuk Analitik Pemasaran dengan R

Unsupervised Learning di R

Terkait

Tutorial Korelasi di R

Spaghetti Plot dan Jalur Badai

40 Pertanyaan Wawancara DBMS Teratas di 2026

12 Alternatif ChatGPT Terbaik yang Bisa Anda Coba pada 2026

Teknologi

Topik

Pelajari lebih lanjut tentang R dan Machine Learning

Machine Learning dengan caret di R

Machine Learning untuk Analitik Pemasaran dengan R

Unsupervised Learning di R

Apa itu Bootstrap?

Bootstrap dalam Praktik

Berbagai Jenis CI Bootstrap

Sedikit notasi (Maaf!)

CI Persentil

CI Normal

CI Basic

BCα CI

Replikasi Hasil

Topik

Pelajari lebih lanjut tentang R dan Machine Learning

Terkait

Pelajari lebih lanjut tentang R dan Machine Learning

BC_α CI