set.seed(123)
# Jumlah observasi
n <- 100
# Generate variabel x dari distribusi normal (mean=10, sd=2)
x <- rnorm(n, mean = 10, sd = 2)
# 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)
# Introduksi missing value secara acak pada 10 observasi x
data[sample(1:n, 10), "x"] <- NA
# Lihat 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
# 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)
## Warning: package 'boot' was built under R version 4.3.3
# 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.581084 0.06067069 1.1482885
## t2* 1.412127 -0.00547455 0.1074228
# Plot distribusi bootstrap
plot(boot_result)
# 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.176, 1.596 )
## Calculations and Intervals on Original Scale
#Praktikum 2
# Hitung mean x (abaikan NA)
mean_x <- mean(data$x, na.rm = TRUE)
# Buat variabel baru dengan imputasi mean
data$ximp <- ifelse(is.na(data$x), mean_x, data$x)
# Fit model 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 setelah imputasi
boot_imp <- function(data, indices) {
d <- data[indices, ]
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.653794 0.053055397 1.1350004
## t2* 1.412127 -0.005093136 0.1064137
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.188, 1.603 )
## Calculations and Intervals on Original Scale