Pegantar Bootstrapping dalam R

Dalam statistik, sering ingin dilakukan resampling untuk menguji keakuratan estimasi sampel. Ini disebut bootstrapping—sebuah tes yang didasarkan pada pengambilan sampel acak berulang denganpengembalian (replacement). Dalam pembahasan kali ini, kita akan:

  1. Menghasilkan 1000 sampel acak baru dari sebuah vektor (dengan pengembalian),
  2. Menghitung estimasi tertentu (misalnya, mean, median, dll.) untuk setiap sampel tersebut, dan
  3. Menghitung standar deviasi dari distribusi estimasi tersebut.

Dataset

Membuat sebuah vektor berisi angka yang terdistribusi secara normal dengan nama ‘Dataku’

Membuat Distribusi Normal Acak Menggunakan rnorm()

Fungsi rnorm() digunakan untuk menghasilkan urutan angka acak dari distribusi normal hipotesis dengan rata-rata (mean) 20 dan standar deviasi 4, sebanyak 4000 observasi.

set.seed(99)
Dataku <- rnorm(4000, 20, 4)

Penjelasan:

  • set.seed(99) digunakan untuk menetapkan seed (angka acak) agar hasil yang dihasilkan oleh fungsi acak (seperti rnorm) dapat direplikasi. Ini memastikan bahwa setiap kali kode dijalankan, hasilnya akan sama.
  • rnorm(4000, 20, 4) menghasilkan 4000 observasi dari distribusi normal dengan mean 20 dan standar deviasi 4. Hasilnya disimpan dalam vektor Dataku.

Melakukan Pemeriksaan Awal (Sanity Checks) pada Dataku Menggunakan length(), mean(), dan sd()

Untuk memastikan bahwa vektor Dataku telah dibuat dengan benar, kita memverifikasi bahwa:

  1. Terdapat 4000 observasi (menggunakan length()),
  2. Nilai rata-rata (mean) mendekati 20 (menggunakan mean()), dan
  3. Standar deviasi mendekati 4 (menggunakan sd()).

Catatan: Kita tidak seharusnya mengharapkan mean dan standar deviasi memiliki nilai yang persis sama dengan yang kita masukkan ke dalam rnorm(), karena ada unsur acak yang terlibat saat menggunakan fungsi ini (rnorm menghasilkan distribusi normal acak di sekitar mean dan standar deviasi yang diberikan).

length(Dataku)
## [1] 4000

length(Dataku) menghitung jumlah observasi dalam vektor Dataku. Harusnya menghasilkan 4000 karena kita membuat 4000 observasi.

mean(Dataku)
## [1] 20.02919

mean(Dataku) menghitung rata-rata dari vektor Dataku. Harusnya mendekati 20 karena kita menetapkan mean = 20 saat membuat data.

sd(Dataku)
## [1] 4.000716

sd(Dataku) menghitung standar deviasi dari vektor Dataku. Harusnya mendekati 4 karena kita menetapkan standar deviasi = 4 saat membuat data.

Membuat Grafik untuk Dataku

Untuk tujuan visual, berikut adalah grafik dari vektor Dataku. Garis putus-putus merah menunjukkan nilai mean (20.02919) dari distribusi ini.

# Membuat histogram dari Dataku
hist(Dataku, breaks = 30, col = "darkolivegreen1", main = "Distribusi Dataku", xlab = "Nilai", ylab = "Frekuensi")

# Menambahkan garis vertikal untuk mean
abline(v = mean(Dataku), col = "red", lwd = 2, lty = 2)

# Menambahkan legenda
legend("topright", legend = paste("Mean =", round(mean(Dataku), 5)), col = "red", lty = 2, lwd = 2)

Penjelasan:

  1. hist() digunakan untuk membuat histogram dari Dataku.
  • breaks = 30 menentukan jumlah bin (kelompok) pada histogram.
  • col = “darkolivegreen1” memberikan warna hijau zaitun muda pada histogram.
  • main, xlab, dan ylab masing-masing digunakan untuk judul grafik, label sumbu x, dan label sumbu y.
  1. abline() menambahkan garis vertikal pada nilai mean.
  • v = mean(Dataku) menentukan posisi garis pada nilai mean.
  • col = “red” memberikan warna merah pada garis.
  • lwd = 2 menentukan ketebalan garis.
  • lty = 2 membuat garis menjadi putus-putus.
  1. legend() menambahkan legenda ke grafik untuk menjelaskan garis merah. Output dari kode ini adalah histogram yang menunjukkan distribusi Dataku dengan garis merah putus-putus yang menandakan nilai mean.

Bootstrapping Menggunakan Fungsi Dasar R

Apa itu bootstrapping? Bootstrapping adalah metode untuk memperhitungkan ketidakpastian saat kita mengukur estimasi sampel. Dengan menggunakan pengambilan sampel acak dengan pengembalian (random sampling with replacement), kita menghitung estimasi kita untuk berbagai sampel acak yang diambil dari distribusi asli kita.

Analogi dengan “Mean of Means” pada Dadu Perhatikan kejadian saat menghitung “rata-rata dari rata-rata” (mean of means) untuk lemparan dadu, kita memperkirakan seberapa besar kemungkinan kita akan mengamati estimasi sampel tertentu (misalnya, mean, standar deviasi). Jika kita mengambil sampel ulang secara acak dari distribusi dasar ini dengan penggantian sejumlah besar kali. Alasan kita mengambil sampel ulang dengan penggantian adalah untuk memungkinkan angka-angka tertentu dipilih beberapa kali, yang menciptakan varians yang lebih besar dalam distribusi yang dihasilkan.

Mengapa Menggunakan Pengembalian (Replacement)? Alasan kita menggunakan pengembalian adalah untuk memungkinkan angka tertentu dipilih lebih dari satu kali, yang menciptakan variasi yang lebih besar dalam distribusi yang dihasilkan.

Langkah-langkah Bootstrapping:

  1. Ambil sampel acak dengan pengembalian dari distribusi asli.
  2. Hitung estimasi (misalnya, mean, median, standar deviasi) untuk sampel tersebut.
  3. Ulangi proses ini berkali-kali (misalnya, 1000 kali) untuk menghasilkan distribusi estimasi.
  4. Analisis distribusi estimasi ini untuk memahami variabilitas dan ketidakpastian dari estimasi sampel.

Resampling dari Dataku Sebanyak 1000 Kali Menggunakan for(i in x)

Secara sederhana, yang kita lakukan di sini adalah menghitung mean dari 1000 sampel yang masing-masing terdiri dari 4000 observasi dari Dataku menggunakan perulangan for(i in x). Mean dari setiap sampel ini disimpan dalam sebuah vektor (bootstrap.results).

Tujuan: Kita tertarik untuk mengukur ketidakpastian dalam distribusi estimasi kita (dalam hal ini, mean). Untuk melakukannya, kita menghitung standar deviasi dari setiap mean yang dihitung menggunakan pengambilan sampel acak dengan pengembalian (random sampling with replacement).

set.seed(99) 
sample.size <- 4000 
n.samples <- 1000 
bootstrap.results <- c() 

for (i in 1:n.samples)
{
  obs <- sample(1:sample.size, replace=TRUE)
  bootstrap.results[i] <- mean(Dataku[obs])
}

Penjelasan:

  • set.seed(99) menetapkan seed untuk replikasi.
  • sample.size <- 4000 menentukan ukuran sampel (sama dengan jumlah observasi dalam Dataku).
  • n.samples <- 1000 menentukan jumlah sampel bootstrap yang akan diambil.
  • bootstrap.results <- c() membuat vektor kosong untuk menyimpan hasil bootstrap.
  • for (i in 1:n.samples) melakukan perulangan sebanyak 1000 kali.
  • obs <- sample(1:sample.size, replace=TRUE) mengambil sampel acak dengan pengembalian dari indeks 1 hingga 4000.
  • bootstrap.results[i] <- mean(Dataku[obs]) menghitung mean dari sampel yang diambil dan menyimpannya dalam vektor bootstrap.results.
length(bootstrap.results)
## [1] 1000

length(bootstrap.results) memeriksa panjang vektor bootstrap.results. Harusnya 1000 karena kita mengambil 1000 sampel.

summary(bootstrap.results)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   19.86   19.98   20.03   20.03   20.07   20.28

summary(bootstrap.results) memberikan ringkasan statistik dari vektor bootstrap.results, termasuk min, max, mean, dan kuartil.

sd(bootstrap.results)
## [1] 0.06486814

sd(bootstrap.results) menghitung standar deviasi dari distribusi mean yang dihasilkan oleh bootstrapping.

par(mfrow=c(2,1), pin=c(5.8,0.98)) 

hist(bootstrap.results, 
     col="#d83737", 
     xlab="Mean", 
     main=paste("Means of 1000 bootstrap samples from Dataku")) 
hist(Dataku,
     col="#37aad8",
     xlab="Value",
     main=paste("Distribution of Dataku"))

Penjelasan:

  • par(mfrow=c(2,1), pin=c(5.8,0.98)) mengatur tata letak plot menjadi 2 baris dan 1 kolom, serta menentukan ukuran plot.
  • hist(bootstrap.results, col=“#d83737”, xlab=“Mean”, main=paste(“Means of 1000 bootstrap samples from Dataku”)) membuat histogram dari distribusi mean hasil bootstrapping.
  • hist(Dataku, col=“#37aad8”, xlab=“Value”, main=paste(“Distribution of Dataku”)) membuat histogram dari data asli Dataku.

Pengambilan sampel ulang sebanyak 1000 kali dari proses pembuatan data aktual menggunakan for(i in x)

Demikian pula, kita juga dapat mengambil 1000 sampel acak dari proses pembuatan data asli.

set.seed(99)
sample.size <- 4000
n.samples <- 1000 
bootstrap.results <- c()

for (i in 1:n.samples)
{
  bootstrap.results[i] <- mean(rnorm(4000,20,4.5)) 
}

Penjelasan:

  • set.seed(99) digunakan untuk menetapkan seed agar hasil yang dihasilkan dapat direplikasi. Ini memastikan bahwa setiap kali kode dijalankan, hasilnya akan sama.
  • sample.size <- 4000 menentukan ukuran sampel yang akan dihasilkan setiap kali (dalam hal ini, 4000 observasi).
  • n.samples <- 1000 menentukan jumlah sampel bootstrap yang akan diambil (dalam hal ini, 1000 sampel).
  • bootstrap.results <- c() membuat vektor kosong untuk menyimpan hasil bootstrap (yaitu mean dari setiap sampel).
  • for (i in 1:n.samples) melakukan perulangan sebanyak 1000 kali.
    • Di dalam perulangan, rnorm(4000, 20, 4) menghasilkan 4000 observasi dari distribusi normal dengan mean 20 dan standar deviasi 4.
    • mean(rnorm(4000, 20, 4)) menghitung mean dari sampel yang baru saja dihasilkan.
    • bootstrap.results[i] menyimpan mean tersebut ke dalam vektor bootstrap.results
length(bootstrap.results)
## [1] 1000

length(bootstrap.results) memeriksa panjang vektor bootstrap.results. Harusnya menghasilkan 1000 karena kita mengambil 1000 sampel bootstrap.

summary(bootstrap.results)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   19.75   19.95   20.00   20.00   20.05   20.24

summary(bootstrap.results) memberikan ringkasan statistik dari vektor bootstrap.results, termasuk nilai minimum, kuartil pertama (Q1), median (Q2), mean, kuartil ketiga (Q3), dan nilai maksimum. Ini membantu kita memahami distribusi dari mean yang dihasilkan oleh bootstrapping.

sd(bootstrap.results)
## [1] 0.06908104

sd(bootstrap.results) menghitung standar deviasi dari distribusi mean yang dihasilkan oleh bootstrapping. Standar deviasi ini mengukur variabilitas atau ketidakpastian dari estimasi mean yang kita peroleh melalui bootstrapping.

par(mfrow=c(2,1), pin=c(5.8,0.98)) 

hist(bootstrap.results, 
     col="#d83737", 
     xlab="Mean",
     main=paste("Means of 1000 bootstrap samples from the DGP")) 

hist(Dataku, 
     col="#37aad8", 
     xlab="Value", 
     main=paste("Distribution of Dataku"))

Penjelasan:

  • par(mfrow=c(2,1), pin=c(5.8,0.98)) mengatur tata letak plot menjadi 2 baris dan 1 kolom, serta menentukan ukuran plot.
  • hist(bootstrap.results, col=“#d83737”, xlab=“Mean”, main=paste(“Means of 1000 bootstrap samples from the DGP”)) membuat histogram dari distribusi mean hasil bootstrapping.
    • col=“#d83737” memberikan warna merah pada histogram.
    • xlab=“Mean” memberikan label “Mean” pada sumbu x.
    • main=paste(“Means of 1000 bootstrap samples from the DGP”) memberikan judul grafik.
  • hist(Dataku, col=“#37aad8”, xlab=“Value”, main=paste(“Distribution of Dataku”)) membuat data asli Dataku.
    • col=“#37aad8” memberikan warna biru pada histogram.
    • xlab=“Value” memberikan label “Value” pada sumbu x.
    • main=paste(“Distribution of Dataku”) memberikan judul grafik

Exercise

Latihan 1

Tetapkan benih Anda pada angka 150. Hasilkan distribusi normal acak dari 1000 observasi, dengan rata-rata 30 dan simpangan baku 2,5. Hitung rata-rata dari 50 sampel dari 1000 observasi dari kumpulan data tersebut. Simpan hasil Anda dalam vektor.

Fungsi yang relevan: set.seed(), rnorm(), for(i in x), sample().

set.seed(150)
data_awal <- rnorm(1000, mean = 30, sd = 2.5)
hasil_mean <- c()
for (i in 1:50) 
{
  sampel <- sample(data_awal, size = 30, replace = TRUE)
  hasil_mean[i] <- mean(sampel)
}
hasil_mean
##  [1] 30.62018 29.40175 29.14913 29.44911 29.61297 30.61114 30.31577 29.62514
##  [9] 29.30315 30.03616 29.51691 30.07556 29.68163 29.45476 31.15973 29.92226
## [17] 29.72954 29.99679 30.03054 29.95485 29.86738 29.12227 29.87347 29.30146
## [25] 30.41955 30.50686 29.34598 30.44362 29.65193 29.87748 30.19688 30.08418
## [33] 30.41547 30.22426 30.50423 29.51167 30.17520 29.73383 29.52756 29.89345
## [41] 30.22087 29.93959 30.41592 29.80145 29.38629 29.81230 30.09678 30.09313
## [49] 29.67071 30.75131

Penjelasan:

hist(hasil_mean, 
     col = "cornsilk2", 
     main = "Distribusi Mean dari 50 Sampel", 
     xlab = "Mean Sampel", 
     ylab = "Frekuensi")

abline(v = mean(hasil_mean), col = "red", lwd = 2, lty = 2)

legend("topright", legend = paste("Rata-rata =", round(mean(hasil_mean), 3)),
       col = "red", lty = 2, lwd = 2)

Latihan 2

Hasilkan dua histogram untuk menampilkan secara grafis distribusi rata-rata yang diperoleh dalam Latihan 1 serta nilai dari 1000 observasi dalam kumpulan data asli Anda. Gabungkan histogram ini menjadi satu grafik keseluruhan.

Fungsi yang relevan: par(), hist().

set.seed(150)
data_awal <- rnorm(1000, mean = 30, sd = 2.5)
hasil_mean <- c()

for (i in 1:50)
{
  sampel <- sample(data_awal, size = 30, replace = TRUE)
  hasil_mean[i] <- mean(sampel)
}

par(mfrow = c(1, 2)) 

# Histogram 1: Data asli
hist(data_awal, 
     col = "darkseagreen3", 
     main = "Distribusi Data Asli", 
     xlab = "Nilai Observasi", 
     ylab = "Frekuensi")

# Histogram 2: Distribusi Mean Sampel
hist(hasil_mean, 
     col = "darkslategray", 
     main = "Distribusi Mean dari Sampel", 
     xlab = "Mean Sampel", 
     ylab = "Frekuensi")

Penjelasan: