Simulasi Estimasi Parameter Menggunakan Regresi Logistik pada Kasus Paparan Merkuri di Sekitar Tambang
Rika Pramestia Putri
2026-04-16
Menentukan Parameter dan Skenario
set.seed(123)
# Parameter model
beta0 <- -1.5
beta1 <- -0.002 # jarak
beta2 <- 0.25 # konsumsi ikan
# Jumlah replikasi simulasi (Monte Carlo)
n_rep <- 50
# Skenario ukuran sampel
n_values <- c(30, 100, 300)
# Skenario outlier
outlier_status <- c("tanpa", "dengan")#Parameter ini menjadi dasar pembentukan data simulasi. Nilai beta1 negatif menunjukkan bahwa jarak menurunkan risiko paparan, sedangkan beta2 positif menunjukkan konsumsi ikan meningkatkan risiko.
Generate data awal
n_awal <- 100
jarak <- runif(n_awal, 50, 1000)
konsumsi <- sample(1:10, n_awal, replace = TRUE)
#Probabilitas dihitung dengan fungsi logistik
lin_pred <- beta0 + beta1*jarak + beta2*konsumsi
p <- 1 / (1 + exp(-lin_pred))
p <- pmin(pmax(p, 0.01), 0.99)
y <- rbinom(n_awal, 1, p)
data_awal <- data.frame(y, jarak, konsumsi)
head(data_awal)## y jarak konsumsi
## 1 0 323.19864 9
## 2 0 798.88988 4
## 3 0 438.52808 6
## 4 0 888.86653 9
## 5 0 943.44392 9
## 6 1 93.27867 7as.data.frame(data_awal)## y jarak konsumsi
## 1 0 323.19864 9
## 2 0 798.88988 4
## 3 0 438.52808 6
## 4 0 888.86653 9
## 5 0 943.44392 9
## 6 1 93.27867 7
## 7 0 551.70021 3
## 8 1 897.79809 8
## 9 0 573.86326 9
## 10 0 483.78400 3
## 11 0 958.99168 7
## 12 0 480.66745 3
## 13 0 693.69210 7
## 14 1 594.00173 6
## 15 1 147.77845 10
## 16 0 904.83372 5
## 17 0 283.78335 5
## 18 0 89.95656 8
## 19 0 361.52468 3
## 20 0 956.77847 10
## 21 0 895.06235 2
## 22 0 708.16324 10
## 23 0 658.48147 2
## 24 0 994.55629 10
## 25 0 672.92051 6
## 26 0 723.10394 4
## 27 0 566.86272 1
## 28 0 614.43492 6
## 29 0 324.70175 3
## 30 1 189.75796 8
## 31 0 964.87302 3
## 32 0 907.18409 8
## 33 0 706.17001 1
## 34 1 805.69405 7
## 35 1 73.38300 7
## 36 0 503.90617 7
## 37 1 770.53656 10
## 38 1 255.58754 6
## 39 0 352.27196 7
## 40 0 270.04450 10
## 41 1 185.66002 5
## 42 0 443.81902 6
## 43 0 443.03811 8
## 44 0 400.40318 5
## 45 0 194.82251 7
## 46 0 181.86576 4
## 47 0 271.38239 3
## 48 0 492.66433 9
## 49 0 302.67401 7
## 50 0 864.93633 6
## 51 1 93.53961 10
## 52 0 470.09007 9
## 53 1 808.97860 7
## 54 1 165.80430 2
## 55 0 582.90058 3
## 56 1 246.20482 8
## 57 1 171.15507 4
## 58 1 765.64247 7
## 59 0 900.29309 4
## 60 0 405.73964 1
## 61 1 681.85943 8
## 62 0 140.09863 4
## 63 0 414.77116 9
## 64 1 310.66446 8
## 65 0 823.90804 6
## 66 0 476.09052 4
## 67 1 819.56114 8
## 68 0 821.77003 3
## 69 0 804.62521 4
## 70 0 467.84010 4
## 71 1 766.75140 6
## 72 0 647.76007 1
## 73 0 724.67328 10
## 74 1 50.59353 4
## 75 1 501.55075 9
## 76 1 259.11294 7
## 77 1 410.82571 8
## 78 0 632.13245 5
## 79 0 384.20801 2
## 80 0 155.57865 6
## 81 1 281.43850 9
## 82 0 684.65281 8
## 83 1 446.76444 10
## 84 1 798.78604 4
## 85 0 147.72141 5
## 86 0 463.14810 7
## 87 0 985.70913 1
## 88 0 898.39856 8
## 89 1 892.14561 8
## 90 1 216.30002 10
## 91 1 174.16091 9
## 92 0 670.44683 8
## 93 1 376.34065 2
## 94 0 673.92022 5
## 95 1 354.35458 9
## 96 0 228.30656 7
## 97 1 793.17959 7
## 98 1 138.91524 10
## 99 1 493.44009 10
## 100 1 535.93019 8#Data yang dihasilkan merepresentasikan hubungan antara faktor risiko (jarak dan konsumsi) dengan status paparan. milai y=1 menunjukkan individu terpapar
Menyimpan Hasil Simulasi
results <- data.frame()#Digunakan untuk menyimpan hasil estimasi parameter dari setiap replikasi simulasi.
Simulasi Monte Carlo, Pembangkitan Data, Estimasi, dan Penyimpanan
for (n in n_values) {
for (out in outlier_status) {
# Monte Carlo simulation with repeated sampling
for (i in 1:n_rep) {
# Generate data
jarak <- runif(n, 50, 1000)
konsumsi <- sample(1:10, n, replace = TRUE)
# Tambahkan outlier
if (out == "dengan") {
jarak[1] <- 2000
}
# Probabilitas
lin_pred <- beta0 + beta1*jarak + beta2*konsumsi
p <- 1 / (1 + exp(-lin_pred))
p <- pmin(pmax(p, 0.01), 0.99)
# Generate Y
y <- rbinom(n, 1, p)
# Model logistik
model <- glm(y ~ jarak + konsumsi, family = binomial)
coef_est <- coef(model)
# Simpan hasil
results <- rbind(results, data.frame(
n = n,
outlier = out,
beta0_hat = coef_est[1],
beta1_hat = coef_est[2],
beta2_hat = coef_est[3]
))
}
}
}## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred#Kode di atas merupakan inti dari simulasi yang dilakukan menggunakan pendekatan Monte Carlo. Proses ini menggabungkan beberapa tahap utama, yaitu pembangkitan data, penambahan outlier, perhitungan probabilitas, pembentukan variabel respon, estimasi model regresi logistik, serta penyimpanan hasil estimasi.
Simulasi dilakukan dengan tiga skenario ukuran sampel (n = 30, 100, 300) dan dua kondisi data (dengan dan tanpa outlier). Untuk setiap kombinasi skenario, proses diulang sebanyak 50 kali guna memperoleh hasil estimasi yang lebih stabil.
Data dibangkitkan secara acak, di mana variabel jarak mengikuti distribusi uniform dan variabel konsumsi ikan bersifat diskrit. Outlier ditambahkan dengan memberikan nilai ekstrem pada variabel jarak. Probabilitas paparan dihitung menggunakan fungsi logistik, kemudian variabel respon dibentuk menggunakan distribusi binomial.
Selanjutnya, model regresi logistik diestimasi menggunakan fungsi glm() dengan distribusi binomial, yang secara implisit menggunakan metode Maximum Likelihood Estimation (MLE). Hasil estimasi parameter dari setiap iterasi kemudian disimpan untuk dianalisis lebih lanjut.
Proses ini bertujuan untuk melihat bagaimana kondisi data yang berbeda mempengaruhi hasil estimasi parameter. Dengan melakukan pengulangan (Monte Carlo), dapat diperoleh gambaran yang lebih akurat mengenai bias dan kestabilan model.
Hasil dari tahap ini berupa kumpulan estimasi parameter untuk setiap kombinasi skenario, yang kemudian digunakan untuk menghitung rata-rata dan bias pada tahap selanjutnya.
Rata-rata Estimasi Parameter
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionsummary_results <- results %>%
group_by(n, outlier) %>%
summarise(
mean_b0 = mean(beta0_hat),
mean_b1 = mean(beta1_hat),
mean_b2 = mean(beta2_hat),
.groups = "drop"
)#Rata-rata ini digunakan untuk melihat seberapa dekat estimasi dengan nilai sebenarnya.
Menghitung Bias
summary_results <- summary_results %>%
mutate(
bias_b0 = mean_b0 - beta0,
bias_b1 = mean_b1 - beta1,
bias_b2 = mean_b2 - beta2
)
as.data.frame(summary_results)## n outlier mean_b0 mean_b1 mean_b2 bias_b0 bias_b1
## 1 30 dengan -4.941376 -0.002518146 0.6514522 -3.44137617 -5.181458e-04
## 2 30 tanpa -1.929550 -0.002762736 0.3354421 -0.42955016 -7.627362e-04
## 3 100 dengan -1.772251 -0.002069195 0.2930153 -0.27225098 -6.919542e-05
## 4 100 tanpa -1.742629 -0.001857062 0.2735748 -0.24262913 1.429380e-04
## 5 300 dengan -1.555965 -0.001895168 0.2504480 -0.05596538 1.048316e-04
## 6 300 tanpa -1.516832 -0.002060239 0.2532775 -0.01683204 -6.023890e-05
## bias_b2
## 1 0.4014521945
## 2 0.0854421183
## 3 0.0430153426
## 4 0.0235748116
## 5 0.0004480072
## 6 0.0032775483#Bias dihitung sebagai selisih antara estimasi dan nilai parameter sebenarnya. Bias kecil → estimasi akurat Bias besar → estimasi kurang baik
Visualisasi Bias
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.5.3# Bias beta1 (jarak)
ggplot(summary_results, aes(x = factor(n), y = bias_b1, fill = outlier)) +
geom_bar(stat = "identity", position = "dodge") +
labs(title = "Bias Parameter Jarak (β1)",
x = "Ukuran Sampel",
y = "Bias") +
theme_minimal()
# Bias beta2 (konsumsi ikan)
ggplot(summary_results, aes(x = factor(n), y = bias_b2, fill = outlier)) +
geom_bar(stat = "identity", position = "dodge") +
labs(title = "Bias Parameter Konsumsi Ikan (β2)",
x = "Ukuran Sampel",
y = "Bias") +
theme_minimal()
#Grafik menunjukkan bahwa:
bias menurun saat sampel meningkat outlier meningkatkan bias
Melihat Data Hasil
head(results)## n outlier beta0_hat beta1_hat beta2_hat
## (Intercept) 30 tanpa -1.285733 -1.358879e-03 0.1553141
## (Intercept)1 30 tanpa -2.168954 -4.506916e-04 0.2508174
## (Intercept)2 30 tanpa -1.257952 -7.519158e-03 0.5085000
## (Intercept)3 30 tanpa -1.472757 -2.182430e-03 0.1996182
## (Intercept)4 30 tanpa -2.504990 -9.483405e-04 0.4064265
## (Intercept)5 30 tanpa -4.702662 -9.808505e-05 0.4015644#Menampilkan sebagian hasil estimasi dari seluruh simulasi.