Pengantar Bootstrapping dalam R

Dalam statistik, kita sering ingin melakukan resampling untuk menguji keakuratan estimasi sampel kita. Ini disebut bootstrapping—sebuah tes yang didasarkan pada pengambilan sampel acak berulang denganpengembalian (replacement). Dalam pembelajaran 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

Kita akan membuat sebuah vektor berisi angka yang terdistribusi secara normal dengan nama myData.

#Membuat Distribusi Normal Acak Menggunakan rnorm()

Fungsi rnorm() memungkinkan kita untuk menghasilkan urutan angka acak dari distribusi normal hipotetisdengan rata-rata (mean) 20 dan standar deviasi 4.5. Di sini, kita memilih untuk menghasilkan 2000 observasi.

set.seed(300) # Setting the seed for replication purposes
myData <- rnorm(2000,20,4.5) # Creating a random normal distribution (n=300, mean=20, sd=4.5)

Penjelasan:

  • set.seed(300) 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(2000, 20, 4.5) menghasilkan 2000 observasi dari distribusi normal dengan mean 20 dan standardeviasi 4.5. Hasilnya disimpan dalam vektor myData.

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

Untuk memastikan bahwa vektor myData telah dibuat dengan benar, kita memverifikasi bahwa: 1. Terdapat 2000 observasi (menggunakan length()), 2. Nilai rata-rata (mean) mendekati 20 (menggunakan mean()), dan 3. Standar deviasi mendekati 4.5 (menggunakan sd()).

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

length(myData) # How many observations?
## [1] 2000

length(myData) menghitung jumlah observasi dalam vektor myData. Harusnya menghasilkan 2000 karena kitamembuat 2000 observasi.

mean(myData) # What is the mean?
## [1] 20.25773

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

sd(myData) # What is the standard deviation?
## [1] 4.590852

sd(myData) menghitung standar deviasi dari vektor myData. Harusnya mendekati 4.5 karena kita menetapkanstandar deviasi = 4.5 saat membuat data.

Membuat Grafik untuk myData

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

# Membuat histogram dari myData
hist(myData, breaks = 30, col = "lightblue", main = "Distribusi myData", xlab = "Nilai", ylab = "Frekuensi")
# Menambahkan garis vertikal untuk mean
abline(v = mean(myData), col = "red", lwd = 2, lty = 2)
# Menambahkan legenda
legend("topright", legend = paste("Mean =", round(mean(myData), 5)), col = "red", lty = 2, lwd = 2)

Penjelasan: 1. hist() digunakan untuk membuat histogram dari myData. breaks = 30 menentukan jumlah bin (kelompok) pada histogram. col = “lightblue” memberikan warna biru muda pada histogram. main, xlab, dan ylab masing-masing digunakan untuk judul grafik, label sumbu x, dan label sumbu y.2. abline() menambahkan garis vertikal pada nilai mean. v = mean(myData) 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. 3. legend() menambahkan legenda ke grafik untuk menjelaskan garis merah. Output dari kode ini adalah histogram yang menunjukkan distribusi myData dengan garis merah putusputus yang menandakan nilai mean.

2. Bootstrapping Menggunakan Fungsi Dasar R

Apa itu bootstrapping?

Bootstrapping adalah metode untuk memperhitungkan ketidakpastian saat kitamengukur estimasi sampel. Dengan menggunakan pengambilan sampel acak dengan pengembalian (randomsampling 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 akanmengamati estimasi sampel tertentu (misalnya, mean, standar deviasi). Jika kita mengambil sampel ulangsecara acak dari distribusi dasar ini dengan penggantian sejumlah besar kali. Alasan kita mengambil sampelulang dengan penggantian adalah untuk memungkinkan angka-angka tertentu dipilih beberapa kali, yangmenciptakan varians yang lebih besar dalam distribusi yang dihasilkan.

Mengapa Menggunakan Pengembalian (Replacement)?

Alasan kita menggunakan pengembalian adalahuntuk memungkinkan angka tertentu dipilih lebih dari satu kali, yang menciptakan variasi yang lebih besardalam 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 myData Sebanyak 1000 Kali Menggunakan for(iin x)

Secara sederhana, yang kita lakukan di sini adalah menghitung mean dari 1000 sampel yang masing-masingterdiri dari 2000 observasi dari myData 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). Untukmelakukannya, kita menghitung standar deviasi dari setiap mean yang dihitung menggunakan pengambilansampel acak dengan pengembalian (random sampling with replacement).

set.seed(200) # Setting the seed for replication purposes
sample.size <- 2000 # Sample size
n.samples <- 1000 # Number of bootstrap samples
bootstrap.results <- c() # Creating an empty vector to hold the results
for (i in 1:n.samples)
{
 obs <- sample(1:sample.size, replace=TRUE)
 bootstrap.results[i] <- mean(myData[obs]) # Mean of the bootstrap sample}
}

Penjelasan:

  • set.seed(200) menetapkan seed untuk replikasi.
  • sample.size <- 2000 menentukan ukuran sampel (sama dengan jumlah observasi dalam myData).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 2000. bootstrap.results[i] <- mean(myData[obs]) menghitung mean dari sampel yang diambil danmenyimpannya dalam vektor bootstrap.results.
length(bootstrap.results) # Sanity check: this should contain the mean of 1000 different samples
## [1] 1000

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

summary(bootstrap.results) # Sanity check
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   19.92   20.19   20.26   20.26   20.33   20.57

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

sd(bootstrap.results) # Checking the standard deviation of the distribution of means (this is what we are interested in!)
## [1] 0.1021229

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

par(mfrow=c(2,1), pin=c(5.8,0.98)) # Combining plots (2 rows, 1 column) and setting the plots size
hist(bootstrap.results, # Creating an histogram
 col="#d83737", # Changing the color
 xlab="Mean", # Giving a label to the x axis
 main=paste("Means of 1000 bootstrap samples from myData")) # Giving a title to the graphhist(myData, # Creating an histogram
hist(myData,  # Creating an histogram
     col = "#37aad8",  # Changing the color
     xlab = "Value",  # Giving a label to the x axis
     main = "Distribution of myData")  # Giving a title to the graph

### Penjelasan:

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

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

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

set.seed(200) # Setting the seed for replication purposes
sample.size <- 2000 # Sample size
n.samples <- 1000 # Number of bootstrap samples
bootstrap.results <- c() # Creating an empty vector to hold the results

for (i in 1:n.samples) {
  bootstrap.results[i] <- mean(rnorm(sample.size, 20, 4.5)) # Mean of the bootstrap sample
}

Penjelasan:

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

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 kitamemahami distribusi dari mean yang dihasilkan oleh bootstrapping.

sd(bootstrap.results) # Checking the standard deviation of the distribution of means (this is what we are interested in!)
## [1] 0.1041927

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)) # Combining plots (2 rows, 1 column) and setting the plots size
hist(bootstrap.results, # Creating an histogram
 col="#d83737", # Changing the color
 xlab="Mean", # Giving a label to the x axis
 main=paste("Means of 1000 bootstrap samples from the DGP")) # Giving a title to the graph
hist(myData, # Creating an histogram
 col="#37aad8", # Changing the color
 xlab="Value", # Giving a label to the x axis
 main=paste("Distribution of myData")) # Giving a title to the graph

Penjelasan:

Exercises

Latihan 1

Tetapkan benih Anda pada angka 150. Hasilkan distribusi normal acak dari 1000 observasi, dengan rata-rata30 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().

Latihan 1: Simulasi Distribusi Normal & Rata-rata Sampel Pada latihan ini, kita akan:

Menetapkan seed ke 150 agar hasil dapat direplikasi. Menghasilkan 1000 observasi dari distribusi normal dengan rata-rata 30 dan simpangan baku 2,5. Mengambil 50 sampel acak dari kumpulan data ini, masing-masing berukuran 50 observasi. Menghitung rata-rata dari setiap sampel dan menyimpannya dalam vektor sample_means. Kode R untuk Latihan 1

# 1. Menetapkan seed untuk memastikan hasil yang replikatif
set.seed(150)  

# 2. Membuat data dengan distribusi normal (1000 observasi, mean = 30, sd = 2.5)
data <- rnorm(1000, mean = 30, sd = 2.5)  

# 3. Inisialisasi vektor kosong untuk menyimpan rata-rata dari setiap sampel
sample_means <- numeric(50)  

# 4. Mengambil 50 sampel acak dari data (masing-masing berukuran 50) dan menghitung rata-ratanya
for (i in 1:50) {
  sample_means[i] <- mean(sample(data, size = 50, replace = TRUE))  
}

# 5. Menampilkan hasil rata-rata dari 50 sampel pertama
sample_means
##  [1] 30.14309 29.02990 29.76690 30.48360 29.73881 29.71241 29.86385 29.50678
##  [9] 30.56253 29.79078 30.12034 29.86928 29.54955 29.55482 30.04611 30.07912
## [17] 30.16797 29.64843 30.13613 30.47296 30.24592 29.79687 29.70930 29.79885
## [25] 30.11966 30.25528 29.48354 29.84970 30.02096 30.38388 29.65810 30.26251
## [33] 30.29870 29.64380 30.04316 30.29065 30.60772 29.66573 29.80451 29.08342
## [41] 29.63084 30.34124 29.88440 29.82377 29.02525 29.77596 29.64970 29.86317
## [49] 29.88076 29.31903

Penjelasan kode:

set.seed(150) → Agar hasil dapat direplikasi. rnorm(1000, mean = 30, sd = 2.5) → Membuat 1000 angka acak dari distribusi normal dengan rata-rata 30 dan simpangan baku 2,5. numeric(50) → Membuat vektor kosong berisi 50 elemen untuk menyimpan rata-rata sampel. for (i in 1:50) → Looping sebanyak 50 kali untuk mengambil sampel acak dari data. sample(data, size = 50, replace = TRUE) → Mengambil 50 angka acak dari data dengan pengembalian. mean() → Menghitung rata-rata dari setiap sampel

Latihan 2

Hasilkan dua histogram untuk menampilkan secara grafis distribusi rata-rata yang diperoleh dalam Latihan 1serta nilai dari 1000 observasi dalam kumpulan data asli Anda. Gabungkan histogram ini menjadi satu grafikkeseluruhan. Fungsi yang relevan: par(), hist().

Membuat histogram untuk distribusi rata-rata dari 50 sampel. Membuat histogram untuk distribusi asli dari 1000 observasi. Menggabungkan kedua histogram dalam satu grafik menggunakan fungsi par(). Kode R untuk Latihan 2

# Mengatur tata letak grafik menjadi 2 baris dan 1 kolom
par(mfrow = c(2, 1))  

# 1. Membuat histogram untuk distribusi rata-rata sampel
hist(sample_means, 
     col = "red", 
     main = "Distribusi Rata-rata dari 50 Sampel",
     xlab = "Rata-rata Sampel",
     ylab = "Frekuensi",
     breaks = 15, 
     border = "black")

# 2. Membuat histogram untuk data asli
hist(data, 
     col = "blue", 
     main = "Distribusi Data Asli (1000 Observasi)",
     xlab = "Nilai Data",
     ylab = "Frekuensi",
     breaks = 30, 
     border = "black")

# Mengembalikan pengaturan grafik ke default (1 grafik per tampilan)
par(mfrow = c(1, 1))

Penjelasan kode:

  • par(mfrow = c(2, 1)) → Mengatur layout grafik agar ada 2 histogram dalam 1 tampilan.
  • hist(sample_means, col = “red”, breaks = 15) → Membuat histogram untuk distribusi rata-rata sampel dengan warna merah.
  • hist(data, col = “blue”, breaks = 30) → Membuat histogram untuk distribusi data asli dengan warna biru.
  • border = “black” → Memberikan garis hitam di setiap batang histogram agar lebih jelas.
  • par(mfrow = c(1, 1)) → Mengembalikan layout ke default agar grafik berikutnya tidak terpengaruh.

Interpretasi Hasil

  • Histogram pertama (merah) menunjukkan distribusi rata-rata dari 50 sampel. Distribusi ini lebih - sempit dibandingkan data asli, karena rata-rata sampel memiliki variabilitas yang lebih kecil.
  • Histogram kedua (biru) menunjukkan distribusi asli dari 1000 observasi, yang lebih lebar dan mencerminkan seluruh data.
  • Kesimpulan: Bootstrapping membantu kita memahami variabilitas estimasi rata-rata dan bagaimana distribusi sampel dibandingkan dengan data asli.