Cara Menentukan Tanggal Kedaluwarsa Obat dengan Regresi Linier dan Monte Carlo
Halo para pejuang Industri Farmasi,
Setelah sebelumnya menulis artikel tentang Cara Menentukan Kedaluarsa Obat dengan Model Efek Campuran Linier, kini saya akan menulis cara menentukan tanggal kedaluarsa yang lebih sederhana. Perhitungan dengan cara ini dapat dilakukan bila variasi antar bets dapat dianggap homogen dan dengan data yang lengkap.
Sebagai contoh kasus, mari kita cermati data uji stabilitas obat berikut ini:
> # menyiapkan data
> batch <- factor(rep(c(1:3), 3))
> time <- rep(c(0, 3, 6, 9, 12, 18), each = 3)
> assay <- c(51, 51, 53,
+ 51, 50, 52,
+ 50, 52, 48,
+ 49, 51, 51,
+ 49, 48, 47,
+ 47, 45, 49)
> data1 <- as.data.frame(cbind(time, assay, batch))
> data1
time assay batch
1 0 51 1
2 0 51 2
3 0 53 3
4 3 51 1
5 3 50 2
6 3 52 3
7 6 50 1
8 6 52 2
9 6 48 3
10 9 49 1
11 9 51 2
12 9 51 3
13 12 49 1
14 12 48 2
15 12 47 3
16 18 47 1
17 18 45 2
18 18 49 3
Dataset ini terdiri dari tiga kolom (tiga variabel) yaitu time atau waktu pengukuran kadar obat (bulan), assay yaitu hasil penetapan kadar obat (mg/tablet), dan nomor batch. Pada data ini, tidak ada data yang hilang dan diasumsikan bahwa variasi antar bets adalah homogen (bisa diuji dengan uji homogenitas varian dengan metode Bartlett atau Levene, dalam tulisan ini di skip dulu ya).
Dalam metode estimasi waktu kedaluarsa ini, kinetika degradasi obat dianggap mengikuti kinetika orde nol. Langkah-langkah estimasinya adalah seperti ini:
pertama, memanggil paket-paket yang dibutuhkan
> # memanggil paket yang diperlukan
> library(dplyr)
> library(MASS)
kedua, menghitung konsentrasi obat yang tersisa dalam satuan persen yaitu dengan membagi nilai assay pada waktu tertentu dengan assay pada waktu ke-nol sesuai dengan masing-masing bets.
> data2 <- data1 %>%
+ group_by(batch) %>% # Kelompokkan berdasarkan bets
+ mutate(cp0 = assay[time == 0][1], # Ambil kadar saat waktu = 0
+ conc = (assay / cp0) * 100) # Hitung persen
> data2
# A tibble: 18 × 5
# Groups: batch [3]
time assay batch cp0 conc
<dbl> <dbl> <dbl> <dbl> <dbl>
1 0 51 1 51 100
2 0 51 2 51 100
3 0 53 3 53 100
4 3 51 1 51 100
5 3 50 2 51 98.0
6 3 52 3 53 98.1
7 6 50 1 51 98.0
8 6 52 2 51 102.
9 6 48 3 53 90.6
10 9 49 1 51 96.1
11 9 51 2 51 100
12 9 51 3 53 96.2
13 12 49 1 51 96.1
14 12 48 2 51 94.1
15 12 47 3 53 88.7
16 18 47 1 51 92.2
17 18 45 2 51 88.2
18 18 49 3 53 92.5
Pada data tersebut yaitu data2, nilai konsentrasi obat yang tersisa (dalam persen) telah dihitung dan dituliskan pada kolom kelima yaitu conc.
ketiga, melakukan analisis regresi linier dan memvisualisasikan hasil fitting.
> # regresi linier
> lm1 <- lm(conc ~ time, data=data2)
> summary(lm1)
Call:
lm(formula = conc ~ time, data = data2)
Residuals:
Min 1Q Median 3Q Max
-6.6202 -0.5396 0.2378 1.2407 4.7746
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 100.2875 1.1470 87.436 < 2e-16 ***
time -0.5169 0.1153 -4.484 0.000376 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.893 on 16 degrees of freedom
Multiple R-squared: 0.5569, Adjusted R-squared: 0.5292
F-statistic: 20.11 on 1 and 16 DF, p-value: 0.0003758
> # visualisasi data
> newx <- seq(min(time), max(25), by = .05)
> conf_interval <- predict(lm1, newdata= data.frame(time=newx),
+ interval = "confidence", level=.95)
> plot(time, data2$conc, ylim=c(80, 120), xlim=c(0, 25), xlab='time (month)',
+ ylab = 'percent remaining')
> lines(newx, conf_interval[,2], col='blue', lty=2)
> lines(newx, conf_interval[,3], col='blue', lty=2)
> abline(lm1)
> abline(h=c(90,110), col=2)
berdasarkan hasil analisis regresi dan visualisasi data, kita dapat mengetahui bahwa terdapat hubungan linier yang signifikan antara konsentrasi obat tersisa dengan waktu penyimpanan.
Sebenarnya, dari plot ini kita bisa mengetahui bahwa secara rata-rata (garis kontinyu warna hitam), waktu kedaluarsa ada di sekitar 20 bulan. Akan tetapi, jika mengacu pada lower 95%-confidence interval, waktu kedaluarsa di sekitar 15 bulan lebih sedikit.
Selanjutnya, kita periksa lebih lanjut dengan uji formal.
keempat, mensimulasikan slope dan intercept (dalam tulisan ini digunakan 10 ribu kali simulasi)
> ## Ambil rata-rata dan variansi kovarians dari koefisien
> coef_mean <- coef(lm1)
> coef_var <- vcov(lm1)
> ## Simulasi 10.000 kombinasi slope dan intercept
> set.seed(123)
> sim_coef <- mvrnorm(n = 10000, mu = coef_mean, Sigma = coef_var)
kelima, menghitung t90 dan mengeluarkan hasil perhitungan negatif
> ## Hitung t90 untuk masing-masing simulasi
> t90_sim <- (90 - sim_coef[,1]) / sim_coef[,2]
> ## Filter hanya t90 positif
> t90_sim <- t90_sim[t90_sim > 0]
keenam, memeriksa ringkasan dan distribusi hasil perhitungan t90
> ## Lihat ringkasan hasil
> summary(t90_sim)
Min. 1st Qu. Median Mean 3rd Qu. Max.
13.25 18.16 19.96 20.63 22.24 64.97
dari ringkasan ini, diketahui bahwa estimasi waktu kedaluarsa sekitar 20 bulan sama seperti pengamatan visual di atas.
> hist(t90_sim)
> hist(log(t90_sim))
berdasarkan pengamatan histogram distribusi frekuensi, diketahui bahwa t90 mengikuti distribusi log-normal. Selain itu, terdapat outlier dalam data ini.
ketujuh, menghilangkan outlier
> ## menghilangkan outlier
> Q1 <- quantile(log(t90_sim), .25)
> Q3 <- quantile(log(t90_sim), .75)
> IQR <- IQR(log(t90_sim))
> no_outliers <- subset(log(t90_sim),
+ log(t90_sim) > (Q1 - 1.5*IQR) &
+ log(t90_sim) < (Q3 + 1.5*IQR))
> hist(no_outliers)
kedelapan, menghitung nilai rata-rata dan confidence interval (95%)
> ## estimasi t90 berdasarkan ci 95%
> log_ci <- quantile(no_outliers, probs = c(0.025, 0.975))
> exp(log_ci)
2.5% 97.5%
15.63450 27.52062
> mean <- mean(no_outliers)
> exp(mean)
[1] 20.10318
Berdasarkan hasil perhitungan ini, nilai rata-rata waktu kedaluarsa (t90) adalah ~20 bulan. Batas bawah interval kepercayaan adalah ~15. Dengan demikian, waktu kedaluarsa dapat dianggap 15 bulan.