set.seed(1)
n <- 100
x <- rnorm(n, mean = 10, sd = 2)
y <- 3 + 1.5 * x +rnorm(n, mean = 0, sd = 2)
data <- data.frame(x, y)
data[sample(1:n, 10), "x"] <- NA
# 1. Bootstrap untuk Regresi
clean_data <- na.omit(data)
boot_regression <- function(data, indices) {
d <- data [indices, ]
model <- lm(y ~ x, data = d)
return(coef(model))
}
library(boot)
## Warning: package 'boot' was built under R version 4.4.3
# Bootstrap dengan 1000 replikasi
boot_result <- boot(
data = clean_data,
statistic = boot_regression,
R = 1000
)
boot_result
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = clean_data, statistic = boot_regression, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 2.330541 0.021682695 0.98665816
## t2* 1.557043 -0.001401521 0.09427229
plot(boot_result)
# Interval kepercayaan 95% untuk koefisien 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.372, 1.737 )
## Calculations and Intervals on Original Scale
# 2. Estimasi pada Missing Value dengan Bootstrap
mean_x <- mean(data$x, na.rm = TRUE)
data$ximp <- ifelse(is.na(data$x), mean_x, data$x) # variabel baru dengan imputasi mean
model_imp <- lm(y ~ ximp, data = data) # fit model setelah imputasi
summary(model_imp)
##
## Call:
## lm(formula = y ~ ximp, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5492 -1.2132 -0.2255 1.0754 4.7252
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.2661 1.1981 1.891 0.0615 .
## ximp 1.5570 0.1151 13.524 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.962 on 98 degrees of freedom
## Multiple R-squared: 0.6511, Adjusted R-squared: 0.6476
## F-statistic: 182.9 on 1 and 98 DF, p-value: < 2.2e-16
boot_imp <- function(data, indices) {
d <- data[indices, ]
model <- lm(y ~ ximp, data = d)
return(coef(model))
}
boot_result_imp <- boot(data = data, statistic = boot_imp, R = 1000)
boot_result_imp
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = data, statistic = boot_imp, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 2.266098 -0.036717339 0.90732541
## t2* 1.557043 0.003549375 0.08795762
plot(boot_result_imp)
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.394, 1.734 )
## Calculations and Intervals on Original Scale
# 3. Multiple Imputation + 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
library(broom)
## Warning: package 'broom' was built under R version 4.4.3
# Multiple imputation sebanyak 5 dengan 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
imp_data <- complete(imp, "long")
# fit model di setiap dataset imputasi dan gabungan hasilnya
model_mi <- with(imp, lm(y ~ x))
summary(pool(model_mi))
## term estimate std.error statistic df p.value
## 1 (Intercept) 2.531091 1.169222 2.164765 92.82244 3.297083e-02
## 2 x 1.530233 0.112619 13.587703 91.37253 1.179867e-23
model_clean <- lm(y ~ x, data = clean_data)
clean_summary <- tidy(model_clean, conf.int = TRUE)
# Model Mean Imputation + Bootstrap
boot_ci <- boot.ci(boot_result_imp, type = "perc", index = 2)
boot_summary <- tidy(model_imp, conf.int = TRUE)
# Model MICE
model_mice <- with(imp, lm(y ~ x))
mice_summary <- summary(pool(model_mice), conf.int = TRUE)
# Data frame yang lebih robust
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]
),
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
)
results_table
## Metode Intercept Slope SE_Slope CI_Slope
## 1 Data Lengkap 2.330541 1.557043 0.1131522 (1.332, 1.782)
## 2 Mean Imputation + Bootstrap 2.266098 1.557043 0.1151300 (1.394, 1.734)
## 3 MICE 2.531091 1.530233 0.1126190 (1.307, 1.754)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
# Plot perbandingan estimasi slope dengan berbagai metode
results <- data.frame(
Method = c("Data Lengkap", "Mean Imp + Bootstrap", "MICE"),
Slope = c(1.412127, 1.412127, 1.408248),
SE = c(0.1079083, 0.1191314, 0.1068028 ),
CI_lower = c(1.198, 1.188, 1.196),
CI_upper = c(1.627, 1.603, 1.621)
)
ggplot(results, aes(x = Method, y = Slope, color = Method)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2) +
labs(title = "Perbandingan Estimasi Slope dengan Berbagai Metode",
y = "Estimasi Slope (y ~ x)") +
theme_minimal()
