Bootstrap dalam Regresi & Estimasi Missing Value

#Simulasi Dataset Awal

set.seed(123)  # Set angka acak agar hasil replikasi tetap sama

n <- 100  # Jumlah total data (observasi)

# Membuat variabel x dari distribusi normal dengan rata-rata 10 dan deviasi standar 2
x <- rnorm(n, mean = 10, sd = 2)

# Membuat variabel y berdasarkan persamaan linear y = 3 + 1.5x ditambah noise acak
y <- 3 + 1.5 * x + rnorm(n, mean = 0, sd = 2)

# Gabungkan x dan y dalam sebuah data frame
data <- data.frame(x, y)

# Secara acak hilangkan nilai pada kolom x sebanyak 10 observasi
data[sample(1:n, 10), "x"] <- NA

# Tampilkan sebagian kecil data (6 baris pertama)
head(data)
##           x        y
## 1  8.879049 14.89776
## 2  9.539645 17.82323
## 3 13.117417 22.18274
## 4 10.141017 17.51644
## 5 10.258575 16.48463
## 6 13.430130 23.05514

Penjelasan Ringkas:

set.seed() berfungsi untuk memastikan proses random menghasilkan output identik jika dijalankan ulang. rnorm() digunakan untuk menghasilkan nilai dari distribusi normal. Nilai y dibuat agar mengikuti hubungan linier dari x dengan tambahan variasi acak. sample() membantu memilih indeks secara acak untuk membuat missing value pada variabel x.

Praktikum 1: Bootstrap Regresi pada Data Tanpa Missing

# Hilangkan baris dengan NA
clean_data <- na.omit(data)

# Definisikan fungsi regresi yang akan dibootstrap
boot_regression <- function(data, indices) {
  sampel <- data[indices, ]  # Ambil subset data berdasarkan indeks bootstrap
  model <- lm(y ~ x, data = sampel)  # Regresi linear
  return(coef(model))  # Kembalikan koefisien regresi
}

# Load pustaka yang dibutuhkan
library(boot)

# Eksekusi bootstrap sebanyak 1000 kali
boot_result <- boot(
  data = clean_data,
  statistic = boot_regression,
  R = 1000
)

# Lihat hasil koefisien dan distribusinya
boot_result
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = clean_data, statistic = boot_regression, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1* 3.581084  0.06067069   1.1482885
## t2* 1.412127 -0.00547455   0.1074228
plot(boot_result)

# Estimasi interval kepercayaan 95% untuk slope (x)
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.176,  1.596 )  
## Calculations and Intervals on Original Scale

Penjelasan Ulang:

na.omit() membuang semua baris yang mengandung NA sehingga analisis dilakukan hanya pada data lengkap. Fungsi boot_regression() digunakan untuk membuat model regresi pada sampel acak (dengan pengembalian) yang diambil dari data asli. Fungsi boot() digunakan untuk menjalankan proses bootstrap non-parametrik. boot.ci() dipakai untuk menghasilkan batas bawah dan atas dari interval kepercayaan berdasarkan hasil bootstrap.

Praktikum 2: Mengatasi Missing dengan Mean Imputation + Bootstrap

# Rata-rata dari x (tidak memperhitungkan NA)
mean_x <- mean(data$x, na.rm = TRUE)

# Gantikan nilai NA dengan rata-rata
data$ximp <- ifelse(is.na(data$x), mean_x, data$x)

# 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 
## -4.1153 -1.4394 -0.0902  1.2053  6.5280 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.6538     1.2332   2.963  0.00383 ** 
## ximp          1.4121     0.1191  11.854  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.109 on 98 degrees of freedom
## Multiple R-squared:  0.5891, Adjusted R-squared:  0.5849 
## F-statistic: 140.5 on 1 and 98 DF,  p-value: < 2.2e-16

Catatan Penjelasan:

Rata-rata x digunakan untuk menggantikan data yang hilang, metode ini disebut mean imputation. Dengan menggunakan ifelse(), NA diganti menjadi rata-rata agar semua baris dapat digunakan. Bootstrap dilakukan seperti sebelumnya namun menggunakan variabel ximp yang telah diimputasi.

Praktikum 3: Multiple Imputation (MICE) + Bootstrap

# Load package mice
library(mice)
## 
## Attaching package: 'mice'
## The following object is masked from 'package:stats':
## 
##     filter
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
# Multiple imputation menggunakan metode Predictive Mean Matching
imp <- mice(data[, c("x", "y")], m = 5, method = "pmm", seed = 123)
## 
##  iter imp variable
##   1   1  x
##   1   2  x
##   1   3  x
##   1   4  x
##   1   5  x
##   2   1  x
##   2   2  x
##   2   3  x
##   2   4  x
##   2   5  x
##   3   1  x
##   3   2  x
##   3   3  x
##   3   4  x
##   3   5  x
##   4   1  x
##   4   2  x
##   4   3  x
##   4   4  x
##   4   5  x
##   5   1  x
##   5   2  x
##   5   3  x
##   5   4  x
##   5   5  x
# Transform hasil imputasi ke format long
imp_data <- complete(imp, "long")

# Jalankan regresi pada masing-masing dataset hasil imputasi
model_mi <- with(imp, lm(y ~ x))

# Gabungkan hasil model dengan teknik pooling
summary(pool(model_mi))
##          term estimate std.error statistic       df      p.value
## 1 (Intercept) 3.619991 1.1112706  3.257524 78.99385 1.657655e-03
## 2           x 1.408248 0.1068028 13.185496 78.10532 1.472407e-21

Penjelasan Singkat:

mice() digunakan untuk melakukan multiple imputation—berbeda dari mean imputation, pendekatan ini mempertahankan variasi alami antar observasi. method = “pmm” menghasilkan prediksi nilai yang realistis berdasarkan nilai-nilai sebenarnya yang mirip. with() digunakan untuk menjalankan model regresi pada masing-masing hasil imputasi, lalu pool() menggabungkan hasil regresi menjadi satu ringkasan.

Ringkasan Gabungan Hasil Analisis

# Library yang dibutuhkan
library(mice)
library(boot)
library(broom)

# Contoh data dengan missing value
set.seed(123)
x <- rnorm(100)
y <- 2 + 3 * x + rnorm(100)
y[sample(1:100, 20)] <- NA  # tambahkan NA
data <- data.frame(x, y)

# Data lengkap (hapus NA)
clean_data <- na.omit(data)

model_clean <- lm(y ~ x, data = clean_data)
clean_summary <- tidy(model_clean, conf.int = TRUE)

# Mean imputation
mean_y <- mean(data$y, na.rm = TRUE)
data_imp <- data
data_imp$y[is.na(data_imp$y)] <- mean_y

# Fungsi bootstrap
boot_fun <- function(data, indices) {
  d <- data[indices, ]
  model <- lm(y ~ x, data = d)
  return(coef(model))
}

# Bootstrap 1000 kali
boot_result_imp <- boot(data = data_imp, statistic = boot_fun, R = 1000)
# Ekstrak hasil koefisien rata-rata
boot_coef <- colMeans(boot_result_imp$t)
# Buat objek tidy buatan agar sesuai format
boot_summary <- data.frame(term = c("(Intercept)", "x"),
                           estimate = boot_coef)

# Imputasi menggunakan MICE
imp <- mice(data, m = 5, method = "pmm", seed = 500)
## 
##  iter imp variable
##   1   1  y
##   1   2  y
##   1   3  y
##   1   4  y
##   1   5  y
##   2   1  y
##   2   2  y
##   2   3  y
##   2   4  y
##   2   5  y
##   3   1  y
##   3   2  y
##   3   3  y
##   3   4  y
##   3   5  y
##   4   1  y
##   4   2  y
##   4   3  y
##   4   4  y
##   4   5  y
##   5   1  y
##   5   2  y
##   5   3  y
##   5   4  y
##   5   5  y
model_mice <- with(imp, lm(y ~ x))
mice_summary <- tidy(pool(model_mice), conf.int = TRUE)

# Buat ringkasan hasil akhir
results_table <- data.frame(
  Metode = c("Data Lengkap", "Mean Imputation + Bootstrap", "Multiple 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])
)

print(results_table)
##                        Metode Intercept    Slope
## 1                Data Lengkap  1.859555 2.949724
## 2 Mean Imputation + Bootstrap  1.963003 2.449028
## 3  Multiple Imputation (MICE)  1.868170 2.913091