Bootsrtap merupakan salah satu metode resampling dengan pengembalian dari data asli. Pada bootstrap dilakukan resampling data asli, kemudian statistik tiap sampel dihitung dan kemudian distribusi statitstik sampel tersebut dapat digunakan untuk memperkirakan ketidakpastian parameter.

# Membangun Data

set.seed(3)

# Jumlah observasi
n <- 100

# Generate variabel x dari distribusi normal (mean=20, sd=5)
x <- rnorm(n, mean = 20, sd = 5)

# Generate variabel y dengan pola hubungan linear terhadap x plus error
y <- 3 + 1.5 * x + rnorm(n, mean = 0, sd = 2)

# Gabungkan menjadi data frame
data <- data.frame(x, y)

# Menjadikan missing value secara acak pada 10 observasi x
data[sample(1:n, 10), "x"] <- NA

# Lihat 6 baris pertama
head(data)
##          x        y
## 1 15.19033 27.02020
## 2 18.53737 29.99590
## 3 21.29394 37.04712
## 4 14.23934 25.56358
## 5 20.97891 36.50329
## 6 20.15062 34.44226

Bootstrap Pada Regresi (tanpa missing value)

set.seed(3)

# Hapus baris yang mengandung NA
clean_data <- na.omit(data)

# Fungsi untuk bootstrap regresi
boot_regression <- function(data, indices) {
  # Ambil sampel bootstrap sesuai indices
  d <- data[indices, ]
  # Fit model regresi linear
  model <- lm(y ~ x, data = d)
  # Return koefisien model
  return(coef(model))
}

# Load library boot
library(boot)

# Lakukan bootstrap dengan 1000 replikasi
boot_result <- boot(data = clean_data, 
                    statistic = boot_regression,
                    R = 1000)

# Tampilkan hasil
boot_result
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = clean_data, statistic = boot_regression, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1* 3.743313 -0.04324745  1.03825132
## t2* 1.468492  0.00231232  0.05084442

Pada bootstrap ini: Ambil sampel dengan penggantian (sampling with repleacment) dari clean_data, buat model regresi (y ~ x), kemudian ambil keofisien regresinya. Ulangi sampai 1000 kali.

Parameter dari data di atas (original) adalah koefisien regresi dari data “clean_data”. Karena pada bootstrap ini kita mengambil statistik koefisien regresinya maka:

t1: Intercept

t2: slope

Rata-rata statistik koefisien intercept memiliki selisih dengan parameternya sebesar 0,04324745. Sedangkan pada rata-rata statistik koefisien slope memiliki selisih dengan parameternya sebesar 0.00231232.

# Plot hasil bootstrap
plot(boot_result)

Distribusi dari statistik koefisien regresi yang dihasilkan adalah distribusi normal yang ditandai dengan plot histogram yang berbentuk menyerupai lonceng serta plot quantile di mana titik-titik plotnya mengikuti garis.

# Hitung confidence interval 95% untuk koefisien x (index=2)
boot.ci(boot_result, type = "perc", index = 2)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_result, type = "perc", index = 2)
## 
## Intervals : 
## Level     Percentile     
## 95%   ( 1.373,  1.571 )  
## Calculations and Intervals on Original Scale

Pada index 2 (slope), parameter koefisien slope berada pada interval (1.373, 1.571) dengan tingkat kepercayaan 95%. Hal ini terbukti bahwa parameternya (1.468492) berada pada interval tersebut.

Bootstrap Pada Regresi (Dengan Missing Value), Input Missing Value Dengan Mean

# Hitung mean x (abaikan NA)
mean_x <- mean(data$x, na.rm = TRUE)

# Buat variabel baru dengan imputasi mean pada missing value x
data$ximp <- ifelse(is.na(data$x), mean_x, data$x)

# Fit model regresi linear setelah imputasi
model_imp <- lm(y ~ ximp, data = data)
summary(model_imp)
## 
## Call:
## lm(formula = y ~ ximp, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.0781  -1.6221   0.0452   1.4950   9.4395 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.59296    1.46420   2.454   0.0159 *  
## ximp         1.46849    0.07138  20.572   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.893 on 98 degrees of freedom
## Multiple R-squared:  0.812,  Adjusted R-squared:  0.8101 
## F-statistic: 423.2 on 1 and 98 DF,  p-value: < 2.2e-16
set.seed(3)

# Fungsi bootstrap setelah imputasi
boot_imp <- function(data, indices) {
  d <- data[indices, ]
  # Model regresi linear
  model <- lm(y ~ ximp, data = d)
  return(coef(model))
  }

# Jalankan bootstrap
boot_result_imp <- boot(data = data, statistic = boot_imp, R = 1000)

# Hasil
boot_result_imp
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = data, statistic = boot_imp, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original        bias    std. error
## t1* 3.592961 -0.0187105280  1.01503001
## t2* 1.468492  0.0008718996  0.05062401

Parameter dari data di atas (original) adalah koefisien regresi dari data yang missing valuenya diisi dengan rata-ratanya. Karena pada bootstrap ini kita juga mengambil statistik koefisien regresinya maka:

t1: Intercept

t2: slope

Rata-rata statistik koefisien intercept memiliki selisih dengan parameternya sebesar 0.0187105280. Sedangkan pada rata-rata statistik koefisien slope memiliki selisih dengan parameternya sebesar 0.0008718996.

# Plot hasil bootstrap
plot(boot_result_imp)

Distribusi statistik koefisien regresi yang dihasilkan juga berdistribusi normal.

# Convidence Interval bootstrap
boot.ci(boot_result_imp, type = "perc", index = 2)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_result_imp, type = "perc", index = 2)
## 
## Intervals : 
## Level     Percentile     
## 95%   ( 1.369,  1.567 )  
## Calculations and Intervals on Original Scale

Pada index 2 (slope), parameter koefisien slope berada pada interval (1.369, 1.567) dengan tingkat kepercayaan 95%. Hal ini terbukti bahwa parameternya (1.468492) berada pada interval tersebut.

Multiplication Impputation + Bootstrap

library(mice)
## Warning: package 'mice' was built under R version 4.4.3
## 
## Attaching package: 'mice'
## The following object is masked from 'package:stats':
## 
##     filter
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
# Lakukan multiple imputation (m=7) dengan Predictive Mean Matching
imp <- mice(data[, c("x", "y")],
            m = 7,
            method = 'pmm',
            seed = 123
            )
## 
##  iter imp variable
##   1   1  x
##   1   2  x
##   1   3  x
##   1   4  x
##   1   5  x
##   1   6  x
##   1   7  x
##   2   1  x
##   2   2  x
##   2   3  x
##   2   4  x
##   2   5  x
##   2   6  x
##   2   7  x
##   3   1  x
##   3   2  x
##   3   3  x
##   3   4  x
##   3   5  x
##   3   6  x
##   3   7  x
##   4   1  x
##   4   2  x
##   4   3  x
##   4   4  x
##   4   5  x
##   4   6  x
##   4   7  x
##   5   1  x
##   5   2  x
##   5   3  x
##   5   4  x
##   5   5  x
##   5   6  x
##   5   7  x
# Gabungkan dataset imputasi dalam long format
imp_data <- complete(imp, "long")

# Fit model di setiap dataset imputasi dan gabungkan hasilnya
model_mi <- with(imp, lm(y ~ x))
summary(pool(model_mi))
##          term estimate  std.error statistic       df      p.value
## 1 (Intercept) 3.834013 1.13187424  3.387314 76.43278 1.117917e-03
## 2           x 1.462267 0.05462417 26.769606 82.04022 2.661416e-42
# Perbandingan Koefisien Model Tanpa Bootstrap

# Library
library(mice)
library(broom)
## Warning: package 'broom' was built under R version 4.4.3
# 1. Model Data Lengkap
model_clean <- lm(y ~ x, data = clean_data)
clean_summary <- tidy(model_clean, conf.int = TRUE)

# 2. Model Mean Imputation
# Asumsi boot_result_imp yang sudah dibuat sebelumnya
boot_ci <- boot.ci(boot_result_imp, type = "perc", index = 2)
boot_summary <- tidy(model_imp, conf.int = TRUE)

# 3. Model MICE
model_mice <- with(imp, lm(y ~ x))
mice_summary <- summary(pool(model_mice), conf.int = TRUE)

# Membuat data frame yang lebih robust
results_table <- data.frame(Metode = c("Data Lengkap", "Mean Imputation", "MICE"),
                            Intercept = c(clean_summary$estimate[1],
                                          boot_summary$estimate[1],
                                          mice_summary$estimate[1]),
                            Slope = c(clean_summary$estimate[2],
                                      boot_summary$estimate[2],
                                      mice_summary$estimate[2]),
                            SE_Slope = c(clean_summary$std.error[2],
                                         boot_summary$std.error[2],
                                         mice_summary$std.error[2]),
                            CI_Slope = c(sprintf("(%.3f, %.3f)", clean_summary$conf.low[2],
                                                 clean_summary$conf.high[2]),
                                         sprintf("(%.3f, %.3f)", boot_ci$percent[4], boot_ci$percent[5]),
                                         sprintf("(%.3f, %.3f)", mice_summary$`2.5 %`[2], mice_summary$`97.5 %`[2])
                                         ),
                            stringsAsFactors = FALSE
                            )
# Tampilkan hasil
print(results_table)
##            Metode Intercept    Slope   SE_Slope       CI_Slope
## 1    Data Lengkap  3.743313 1.468492 0.05489248 (1.359, 1.578)
## 2 Mean Imputation  3.592961 1.468492 0.07138316 (1.369, 1.567)
## 3            MICE  3.834013 1.462267 0.05462417 (1.354, 1.571)

Dari output di atas terlihat nilai statitistik koefeisien regresi slope dari data tanpa missing value dan data dengan mengganti missing value dengan nilai mean adalah sama, yaitu 1.468492. Hal ini menandakan bahwa ketika kita mengganti missing value dengan rata-rata maka nilai koefisien slopenya tidak berubah, yang berubah hanya koefisien intersept-nya.