# Model: P(Diabetes = 1) = logit(b0 + b1*GulaDarah + b2*BMI)
# Parameter sebenarnya (true):
# b0 = -5.0 (intercept)
# b1 = 0.03 (efek kadar gula darah, mg/dL)
# b2 = 0.10 (efek BMI, kg/m²)
# Skenario:
# 1. Variasi ukuran sampel (n = 50, 200, 1000)
# 2. Adanya data salah catat / outlier (10% nilai dilipatgandakan)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
# PARAMETER SEBENARNYA
b0 <- -5.0 # intercept
b1 <- 0.03 # koefisien Gula Darah (mg/dL)
b2 <- 0.10 # koefisien BMI (kg/m²)
cat("STUDI KASUS: PREDIKSI RISIKO DIABETES\n")
## STUDI KASUS: PREDIKSI RISIKO DIABETES
cat("Parameter sebenarnya:\n")
## Parameter sebenarnya:
cat(" b0 (intercept) :", b0, "\n")
## b0 (intercept) : -5
cat(" b1 (Gula Darah) :", b1, "\n")
## b1 (Gula Darah) : 0.03
cat(" b2 (BMI) :", b2, "\n\n")
## b2 (BMI) : 0.1
# FUNGSI BANGKITKAN DATA
# Gula Darah ~ N(120, 25) mg/dL rentang realistis pasien
# BMI ~ N(27, 5) kg/m² rentang realistis populasi
buat_data <- function(n, seed, outlier = FALSE) {
set.seed(seed)
GulaDarah <- rnorm(n, mean = 120, sd = 25)
BMI <- rnorm(n, mean = 27, sd = 5)
prob <- 1 / (1 + exp(-(b0 + b1*GulaDarah + b2*BMI)))
Y <- rbinom(n, 1, prob) # 1 = diabetes, 0 = tidak
# Outlier: 10% nilai gula darah salah catat (dilipatgandakan)
if (outlier) {
idx <- sample(1:n, round(n * 0.1))
GulaDarah[idx] <- GulaDarah[idx] * 2
}
data.frame(Y, GulaDarah, BMI)
}
# FUNGSI ESTIMASI
estimasi <- function(df) {
tryCatch(coef(glm(Y ~ GulaDarah + BMI, data = df, family = binomial)),
error = function(e) c(NA, NA, NA))
}
# REPLIKASI 100x
replikasi <- function(n, B = 100, outlier = FALSE) {
t(sapply(1:B, function(b) estimasi(buat_data(n, seed = b*100, outlier))))
}
# RUN SEMUA SKENARIO
cat("Menjalankan simulasi 100 replikasi\n")
## Menjalankan simulasi 100 replikasi
r_n50 <- replikasi(n = 50) # S1a: pasien sedikit
r_n200 <- replikasi(n = 200) # S1b: pasien sedang
r_n1000 <- replikasi(n = 1000) # S1c: pasien banyak
r_clean <- replikasi(n = 200) # S2a: data bersih
r_out <- replikasi(n = 200, outlier = TRUE) # S2b: ada salah catat
cat("Selesai!\n\n")
## Selesai!
# HITUNG BIAS, VARIANSI, MSE
true_val <- c(b0, b1, b2)
ringkasan <- function(mat, label) {
mean_est <- colMeans(mat, na.rm = TRUE)
bias <- mean_est - true_val
variansi <- apply(mat, 2, var, na.rm = TRUE)
data.frame(
Skenario = label,
Parameter = c("b0 (Intercept)", "b1 (Gula Darah)", "b2 (BMI)"),
True = true_val,
Mean_Est = round(mean_est, 4),
Bias = round(bias, 4),
Variansi = round(variansi, 4),
MSE = round(bias^2 + variansi, 4)
)
}
hasil <- rbind(
ringkasan(r_n50, "S1: n=50 (Pasien Sedikit)"),
ringkasan(r_n200, "S1: n=200 (Pasien Sedang)"),
ringkasan(r_n1000, "S1: n=1000 (Pasien Banyak)"),
ringkasan(r_clean, "S2: Tanpa Salah Catat"),
ringkasan(r_out, "S2: Salah Catat 10%")
)
cat("=== TABEL RINGKASAN HASIL SIMULASI ===\n")
## === TABEL RINGKASAN HASIL SIMULASI ===
print(hasil, row.names = FALSE)
## Skenario Parameter True Mean_Est Bias Variansi
## S1: n=50 (Pasien Sedikit) b0 (Intercept) -5.00 -5.6563 -0.6563 10.2971
## S1: n=50 (Pasien Sedikit) b1 (Gula Darah) 0.03 0.0338 0.0038 0.0003
## S1: n=50 (Pasien Sedikit) b2 (BMI) 0.10 0.1138 0.0138 0.0065
## S1: n=200 (Pasien Sedang) b0 (Intercept) -5.00 -5.1731 -0.1731 2.6641
## S1: n=200 (Pasien Sedang) b1 (Gula Darah) 0.03 0.0315 0.0015 0.0001
## S1: n=200 (Pasien Sedang) b2 (BMI) 0.10 0.1009 0.0009 0.0016
## S1: n=1000 (Pasien Banyak) b0 (Intercept) -5.00 -5.1003 -0.1003 0.2928
## S1: n=1000 (Pasien Banyak) b1 (Gula Darah) 0.03 0.0299 -0.0001 0.0000
## S1: n=1000 (Pasien Banyak) b2 (BMI) 0.10 0.1041 0.0041 0.0002
## S2: Tanpa Salah Catat b0 (Intercept) -5.00 -5.1731 -0.1731 2.6641
## S2: Tanpa Salah Catat b1 (Gula Darah) 0.03 0.0315 0.0015 0.0001
## S2: Tanpa Salah Catat b2 (BMI) 0.10 0.1009 0.0009 0.0016
## S2: Salah Catat 10% b0 (Intercept) -5.00 -2.9527 2.0473 1.9940
## S2: Salah Catat 10% b1 (Gula Darah) 0.03 0.0128 -0.0172 0.0000
## S2: Salah Catat 10% b2 (BMI) 0.10 0.0942 -0.0058 0.0014
## MSE
## 10.7278
## 0.0003
## 0.0067
## 2.6941
## 0.0001
## 0.0016
## 0.3029
## 0.0000
## 0.0002
## 2.6941
## 0.0001
## 0.0016
## 6.1856
## 0.0003
## 0.0014
# Skenario terbaik & terburuk
bias_rank <- aggregate(abs(Bias) ~ Skenario, data = hasil, FUN = mean)
names(bias_rank)[2] <- "Mean_AbsBias"
bias_rank <- bias_rank[order(-bias_rank$Mean_AbsBias), ]
cat("\n=== RANKING SKENARIO (dari terburuk ke terbaik) ===\n")
##
## === RANKING SKENARIO (dari terburuk ke terbaik) ===
print(bias_rank, row.names = FALSE)
## Skenario Mean_AbsBias
## S2: Salah Catat 10% 0.69010000
## S1: n=50 (Pasien Sedikit) 0.22463333
## S1: n=200 (Pasien Sedang) 0.05850000
## S2: Tanpa Salah Catat 0.05850000
## S1: n=1000 (Pasien Banyak) 0.03483333
cat("\nSkenario bias TERBESAR :", bias_rank$Skenario[1], "\n")
##
## Skenario bias TERBESAR : S2: Salah Catat 10%
cat("Skenario estimasi TERBAIK:", bias_rank$Skenario[nrow(bias_rank)], "\n\n")
## Skenario estimasi TERBAIK: S1: n=1000 (Pasien Banyak)
# Grafik 1 — Distribusi estimasi b1 (Gula Darah), Skenario 1
df1 <- rbind(
data.frame(b1 = r_n50[,2], n = "n=50"),
data.frame(b1 = r_n200[,2], n = "n=200"),
data.frame(b1 = r_n1000[,2], n = "n=1000")
)
ggplot(df1, aes(x = b1, fill = n)) +
geom_density(alpha = 0.5) +
geom_vline(xintercept = b1, color = "red", linetype = "dashed", linewidth = 1) +
labs(title = "Distribusi Estimasi b1 (Gula Darah) — Variasi Ukuran Sampel",
subtitle = "Garis merah = nilai sebenarnya b1 = 0.03",
x = "Estimasi b1", y = "Densitas", fill = "Ukuran Sampel") +
theme_bw()
# Grafik 2 — Distribusi estimasi b1 (Gula Darah), Skenario 2
df2 <- rbind(
data.frame(b1 = r_clean[,2], kondisi = "Tanpa Salah Catat"),
data.frame(b1 = r_out[,2], kondisi = "Salah Catat 10%")
)
ggplot(df2, aes(x = b1, fill = kondisi)) +
geom_density(alpha = 0.5) +
geom_vline(xintercept = b1, color = "red", linetype = "dashed", linewidth = 1) +
labs(title = "Distribusi Estimasi b1 (Gula Darah) — Pengaruh Salah Catat",
subtitle = "Garis merah = nilai sebenarnya b1 = 0.03",
x = "Estimasi b1", y = "Densitas", fill = "Kondisi Data") +
theme_bw()
# Grafik 3 — Bias per skenario
ggplot(hasil, aes(x = Skenario, y = Bias, fill = Parameter)) +
geom_col(position = "dodge") +
geom_hline(yintercept = 0) +
labs(title = "Bias Estimasi Parameter per Skenario",
subtitle = "Studi Kasus: Prediksi Risiko Diabetes", x = NULL) +
theme_bw() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))
# Grafik 4 — MSE per skenario
ggplot(hasil, aes(x = Skenario, y = MSE, fill = Parameter)) +
geom_col(position = "dodge") +
labs(title = "MSE Estimasi Parameter per Skenario",
subtitle = "Studi Kasus: Prediksi Risiko Diabetes", x = NULL) +
theme_bw() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))
