Metode Regresi Linier Bayesian untuk Menentukan Kedaluarsa Obat
Dalam artikel sebelumnya saya telah menulis beberapa metode untuk menentukan waktu kedaluarsa obat berdasarkan regresi linier sederhana. Anda dapat melihat kembali artikel tersebut pada tautan ini. Metode yang telah saya paparkan sebelumnya adalah metode berdasarkan pendekatan frekuentis yang memerlukan data dalam jumlah cukup besar dan lengkap.
Kali ini saya akan menuliskan cara menentukan waktu kedaluarsa berdasarkan pendekatan Bayesian. Pendekatan Bayesian sendiri merupakan hal yang baru bagi saya akan tetapi ini sangat menarik karena memungkinkan pengambilan sampel yang lebih kecil.
Dalam artikel ini saya hanya akan memberikan panduan praktisnya dalam analisis dengan R. Mari kita mulai.
Kinetika degradasi dianggap mengikuti orde nol dengan persamaan:
$$y_n=\alpha - \beta x_n + \epsilon_n$$
di mana y adalah persen obat yang tersisa, alfa adalah kadar obat mula-mula (C0), beta adalah konstanta degradasi obat orde nol, x adalah waktu sejak pembuatan obat (dalam bulan), dan epsilon adalah residual error.
residual dianggap terdistribusi normal sehingga:
$$\epsilon_n \sim \operatorname{normal}(0,\sigma),$$
karena resiudal adalah selisih dari data eksperimental dengan data prediksi maka persamaan di atas ekivalen dengan:
$$y_n - (\alpha + \beta X_n) \sim \operatorname{normal}(0,\sigma),$$
ini masih bisa disederhanakan lagi menjadi:
$$y_n \sim \operatorname{normal}(\alpha + \beta X_n, \, \sigma).$$
persmaaan inilah yang menjadi dasar dalam pemodelan dengan statistika Bayesian.
Sekarang, mari kita coba implementasikan di R.
Pertama, implementasi metode Bayesian di R dapat dilakukan dengan memanggil library rstan.
> # memanggil library rstan
> library(rstan)
> rstan_options(auto_write = TRUE)
Kedua, menginput data. Dalam pemrograman ini, kita akan membuat data sebagai list bukan sebagai data.frame. Untuk sementara ini, kita abaikan dulu pengaruh bets (dianggap antar bets seragam), sehingga data nomor bets tidak dimasukkan.
> # memanggil library rstan
> library(rstan)
> rstan_options(auto_write = TRUE)
> data<- list(x= x, y = y, N = length(x))
> data
$x
[1] 0 0 0 3 3 3 6 6 6 9 9 9 12 12 12 18 18 18
$y
[1] 100.0 100.0 100.0 100.0 98.0 98.1 98.0 102.0 90.6 96.1 100.0 96.2
[13] 96.1 94.1 88.7 92.2 88.2 92.5
$N
[1] 18
Ketiga, membuat struktur model Bayesian dalam format stan. Kita dapat membuatnya dalam notepad kemudian disimpan dalam ekstensi .stan.
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
}
model {
y ~ normal(alpha + beta * x, sigma);
}
Selanjutnya, model ini dikompilasi ke dalam R.
> linear_regression <- stan_model("linear_regression.stan")
Keempat, melakukan fitting model terhadap data.
> fit1 <- sampling(linear_regression, data = data, chains = 4,
+ iter = 1000, refresh = 0)
> print(fit1)
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha 100.30 0.04 1.22 97.86 99.52 100.31 101.10 102.68 861 1
beta -0.52 0.00 0.12 -0.76 -0.60 -0.52 -0.44 -0.27 821 1
sigma 3.13 0.02 0.61 2.20 2.71 3.04 3.43 4.55 968 1
lp__ -27.64 0.04 1.28 -31.02 -28.25 -27.30 -26.72 -26.19 842 1
Samples were drawn using NUTS(diag_e) at Thu May 8 09:52:49 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
lalu memastikan apakah kalkulasi parameter sudah stasioner
> traceplot(fit1)
dari plot ini kita tahu bahwa perhitungan parameter sudah stasioner sehingga bisa lanjut ke tahap selanjutnya.
Kelima, menghitung t90 (waktu kedaluarsa)
> # mengekstrak distribusi posterior dari parameter
> posterior <- extract(fit1, pars = c('alpha', 'beta'))
> # menghitung t90
> t90 <- (90 - posterior$alpha)/posterior$beta
> summary_t90 <- quantile(t90, probs=c(0.025, 0.5, 0.975))
> summary_t90
2.5% 50% 97.5%
15.31435 19.91444 31.39603
berdasarkan perhitungan ini, kita tahu bahwa median waktu kedaluarsa adalah ~20 bulan dengan persentil 2,5% (dianggap batas bawah) adalah ~15 bulan. Oleh karena itu, waktu kedaluarsa ditetapkan 15 bulan.
keenam, visualisai data
> # Ambil 1000 nilai posterior dari alpha dan beta
> alpha <- posterior$alpha
> beta <- posterior$beta
> # Buat grid x (waktu) untuk prediksi, misal dari 0 sampai 20 bulan
> x_pred <- seq(0, 40, by = 0.1)
> # Buat matriks prediksi: untuk setiap x_pred, hitung y_pred dari setiap pasangan (alpha, beta)
> y_pred_matrix <- sapply(x_pred, function(x_val) alpha + beta * x_val)
> # Hitung kuantil 2.5%, 50%, dan 97.5% untuk setiap titik waktu prediksi
> y_lower <- apply(y_pred_matrix, 2, quantile, probs = 0.025)
> y_median <- apply(y_pred_matrix, 2, quantile, probs = 0.5)
> y_upper <- apply(y_pred_matrix, 2, quantile, probs = 0.975)
> # Plot hasil
> library(ggplot2)
> # Data asli
> data_plot <- data.frame(x = x, y = y)
> # Data prediksi
> pred_plot <- data.frame(
+ x = x_pred,
+ lower = y_lower,
+ median = y_median,
+ upper = y_upper
+ )
> ggplot() +
+ geom_point(data = data_plot, aes(x = x, y = y), color = "black", size = 2) +
+ geom_line(data = pred_plot, aes(x = x, y = median), color = "blue", size = 1) +
+ geom_ribbon(data = pred_plot, aes(x = x, ymin = lower, ymax = upper), fill = "skyblue", alpha = 0.4) +
+ geom_line(data = pred_plot, aes(x = x, y = lower), color = "blue", size = .5, lty=2) +
+ geom_line(data = pred_plot, aes(x = x, y = upper), color = "blue", size = .5, lty=2) +
+
+ labs(
+ title = "Fitting Data Stabilitas Obat dengan Regresi Linear Bayesian",
+ x = "Waktu (bulan)",
+ y = "Persen Obat Tersisa"
+ ) +
+ geom_hline(yintercept=90, col='red', lty=2)+
+ theme_minimal()
dari visualisasi data ini terlihat bahwa model dengan data sudah cukup fit. Dan waktu kedaluarsa hasil estimasi baik median maupun pada batas bawah sudah sesuai.
sekian dulu penjelasan tentang perhitungan waktu kedaluarsa dengan metode Bayesian ini. Mudah-mudahan bermanfaat.