library(boot)
## Warning: package 'boot' was built under R version 4.4.3
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
library(broom)
## Warning: package 'broom' was built under R version 4.4.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
1. Membuat Data Simulasi
set.seed(123)
# Jumlah observasi
n <- 100
# Variabel x
x <- rnorm(n, mean = 10, sd = 2)
# Variabel y
y <- 3 + 1.5 * x + rnorm(n, mean = 0, sd = 2)
# Gabungkan menjadi data frame
data <- data.frame(x, y)
# Menambahkan missing value pada x
data[sample(1:n, 10), "x"] <- NA
# Melihat data awal
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
2. Bootstrap Regresi Tanpa Missing Value
# Menghapus data yang mengandung NA
clean_data <- na.omit(data)
# Fungsi bootstrap regresi
boot_regression <- function(data, indices) {
d <- data[indices, ]
model <- lm(y ~ x, data = d)
return(coef(model))
}
# Bootstrap 1000 kali
boot_result <- boot(
data = clean_data,
statistic = boot_regression,
R = 1000
)
# Hasil bootstrap
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 bootstrap
plot(boot_result)
# Confidence Interval 95%
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
3. Mean Imputation + Bootstrap
# Menghitung mean x tanpa NA
mean_x <- mean(data$x, na.rm = TRUE)
# Mengganti NA dengan mean
data$ximp <- ifelse(
is.na(data$x),
mean_x,
data$x
)
# Model regresi 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
# Fungsi bootstrap imputasi
boot_imp <- function(data, indices) {
d <- data[indices, ]
model <- lm(y ~ ximp, data = d)
return(coef(model))
}
# Bootstrap setelah imputasi
boot_result_imp <- boot(
data = data,
statistic = boot_imp,
R = 1000
)
# Hasil bootstrap
boot_result_imp
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = data, statistic = boot_imp, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 3.653794 0.053055397 1.1350004
## t2* 1.412127 -0.005093136 0.1064137
# Plot bootstrap
plot(boot_result_imp)
# Confidence Interval
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.188, 1.603 )
## Calculations and Intervals on Original Scale
4. Multiple Imputation (MICE)
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
# Model regresi pada hasil imputasi
model_mice <- with(
imp,
lm(y ~ x)
)
# Ringkasan hasil
summary(
pool(model_mice)
)
## 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
5. Perbandingan Hasil
# Model data lengkap
model_clean <- lm(
y ~ x,
data = clean_data
)
clean_summary <- tidy(
model_clean,
conf.int = TRUE
)
# Model mean imputation
boot_summary <- tidy(
model_imp,
conf.int = TRUE
)
# Confidence interval bootstrap
boot_ci <- boot.ci(
boot_result_imp,
type = "perc",
index = 2
)
# Model MICE
mice_summary <- summary(
pool(model_mice),
conf.int = TRUE
)
# Tabel hasil
results_table <- data.frame(
Metode = c(
"Data Lengkap",
"Mean Imputation + Bootstrap",
"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]
)
)
# Menampilkan hasil
print(results_table)
## Metode Intercept Slope
## 1 Data Lengkap 3.581084 1.412127
## 2 Mean Imputation + Bootstrap 3.653794 1.412127
## 3 MICE 3.619991 1.408248
6. Visualisasi Perbandingan
results <- data.frame(
Method = c(
"Data Lengkap",
"Mean Imp + Bootstrap",
"MICE"
),
Slope = c(
clean_summary$estimate[2],
boot_summary$estimate[2],
mice_summary$estimate[2]
)
)
ggplot(
results,
aes(x = Method, y = Slope)
) +
geom_point(size = 3) +
labs(
title = "Perbandingan Estimasi Slope",
y = "Koefisien Slope"
) +
theme_minimal()
