Pendahuluan

Dalam dunia farmasi, waktu kedaluarsa (shelf-life) adalah janji mutu dari produsen kepada pasien: selama periode ini, kualitas, keamanan, dan efektivitas obat tetap terjamin. Setelah melewati batas ini, produsen tidak lagi bertanggung jawab atas mutu produk tersebut. Itulah sebabnya, sangat penting bagi apoteker, tenaga kesehatan, maupun pasien untuk selalu memeriksa tanggal kedaluarsa sebelum menggunakan obat.

Penetapan waktu kedaluarsa tidak dilakukan sembarangan. Industri farmasi wajib mengikuti standar yang ketat, salah satunya yang diatur dalam International Conference on Harmonization (ICH) Q1A (R2). Pedoman ini menguraikan prosedur lengkap, mulai dari desain uji stabilitas, kondisi penyimpanan, hingga analisis data untuk memastikan keakuratan estimasi waktu kedaluarsa.

Bagaimana Waktu Kedaluarsa Ditetapkan?

Dalam analisis data stabilitas, salah satu pendekatan yang banyak digunakan adalah model efek campuran linier (linear mixed effect model). Model ini sangat berguna ketika kita memiliki data dari lima batch atau lebih, dan memperlakukan batch sebagai faktor acak (random effect). Pendekatan ini mencerminkan variasi alami antar batch, membuat hasil estimasi menjadi lebih realistis dan andal.

Masa simpan produk obat ditentukan berdasarkan konsep statistik yang cukup sederhana namun kuat:

Kita ingin yakin, dengan tingkat kepercayaan 90%, bahwa paling tidak 95% dari semua batch masih memenuhi standar kualitas (Critical Quality Attributes, CQA).
Secara visual, ini digambarkan sebagai batas bawah dari grafik prediksi. Selama batas bawah ini masih berada di atas ambang batas mutu (acceptance limit), maka obat masih dianggap stabil.

Analogi Sederhana

Bayangkan Anda memiliki 100 botol obat dari berbagai batch. Tugas Anda adalah memastikan, sampai bulan ke berapa obat-obat ini masih bisa disimpan, sambil tetap yakin bahwa minimal 95 botol tetap memenuhi standar.

Namun, Anda juga menginginkan keputusan ini cukup pasti — setidaknya 90% yakin.

Dengan model efek campuran, Anda dapat memperkirakan kapan kualitas "terburuk" dari obat mulai menyentuh batas mutu. Titik itulah yang kita sebut sebagai akhir masa simpan.

Contoh Kasus

Mari kita lihat sebuah contoh kasus sederhana untuk memahami penerapan konsep ini.

Misalnya, sebuah studi stabilitas dilakukan dengan menyimpan produk obat selama 6 bulan dalam kondisi real-time, dan dilakukan pengukuran kadar zat aktif pada beberapa waktu tertentu. Data dikumpulkan dari lima batch berbeda.

> runOrder <- 1:14
> month <- c(rep(c(0,3), each=5), rep(6,4))
> batch <- factor(c(4, 3, 5, 2, 1, 2, 4, 1, 5, 3, 2, 5, 1, 3))
> concDrug <- c(99.27, 100.026, 100.403, 100.108, 99.903, 
+               99.478, 99.759, 99.432, 100.596, 99.619, 
+               98.001, 98.957, 99.048, 99.808)
> dataStab <- as.data.frame(cbind(runOrder, month, batch, concDrug))
> dataStab
   runOrder month batch concDrug
1         1     0     4   99.270
2         2     0     3  100.026
3         3     0     5  100.403
4         4     0     2  100.108
5         5     0     1   99.903
6         6     3     2   99.478
7         7     3     4   99.759
8         8     3     1   99.432
9         9     3     5  100.596
10       10     3     3   99.619
11       11     6     2   98.001
12       12     6     5   98.957
13       13     6     1   99.048
14       14     6     3   99.808

# runOrder : urutan pengujian; month : waktu pengujian; 
# batch : nomor bets; concDrug : konsentrasi obat yang tersisa

Menggunakan model efek campuran linier, kita dapat mengestimasi penurunan kadar zat aktif dari waktu ke waktu, serta menentukan kapan batas kualitas minimum tercapai.

Estimasi Waktu Kedaluarsa Berdasarkan Kinetika Orde Nol

Berdasarkan data yang tersedia, kita melakukan estimasi waktu kedaluarsa menggunakan model kinetika orde nol, yang dirumuskan sebagai:

$$C = C_0 - k \times t$$

Di mana:

C = konsentrasi obat pada waktu t (persen)
C0 = konsentrasi awal (persen)
k = laju penurunan konsentrasi (persen per bulan)
t = waktu yang dihitung sejak obat selesai dibuat (bulan)

Waktu kedaluarsa, yaitu saat konsentrasi turun hingga 90% dari kadar awal, dapat dihitung dengan rumus:

$$t_{90} = \frac{90 - C_0}{k}$$

Namun, dalam studi ini ada data yang tidak lengkap — tepatnya, data pengujian bulan ke-6 pada batch ke-4 tidak tersedia. Kondisi ini menyebabkan metode konvensional (menghitung regresi per batch lalu dirata-rata) tidak bisa digunakan.

Sebagai gantinya, pendekatan yang lebih tepat adalah menggunakan model efek campuran linier yang bisa mengatasi ketidaklengkapan data.

Alur Perhitungan Waktu Kedaluarsa

Langkah-langkah analisisnya adalah sebagai berikut:

Panggil semua paket yang diperlukan untuk analisis data

 > # load packages
 > library(nlme)
 > library(lme4)
 > library(dplyr)
 > library(ggplot2)
 > library(lattice)
 > library(outliers)
 > library(tidyverse)
 > library(modelr)
 > library(arm)

Menyiapkan data uji stabilitas.

 # membuat data terkelompok
 dataStabGrouped <- groupedData(concDrug ~ month | batch,
                         data = dataStab,
                         labels=list(x="time (month)",
                                     y = "remaining concentration (%)"))
 # plotting data
 xyplot(concDrug ~ month | batch, data = dataStabGrouped,
        type = c("p", "r"))

Membuat model efek campuran:

Menyertakan slope (k) dan intercept (C0) sebagai efek acak, dan menganggap keduanya saling berkorelasi.

  > # fitting with linear mixed effect model  with random slopes and intercepts 
  > model_lme_01 <- lmer(concDrug ~ month + (month | batch),
  +                  data = dataStabGrouped)
  > summary(model_lme_01)
  Linear mixed model fit by REML ['lmerMod']
  Formula: concDrug ~ month + (month | batch)
     Data: dataStabGrouped

  REML criterion at convergence: 25.6

  Scaled residuals: 
      Min      1Q  Median      3Q     Max 
  -1.1938 -0.3105 -0.1690  0.2685  2.1239 

  Random effects:
   Groups   Name        Variance Std.Dev. Corr 
   batch    (Intercept) 0.08874  0.2979        
            month       0.01719  0.1311   -0.84
   Residual             0.18090  0.4253        
  Number of obs: 14, groups:  batch, 5

  Fixed effects:
               Estimate Std. Error t value
  (Intercept) 100.03268    0.22004 454.620
  month        -0.14574    0.07724  -1.887

  Correlation of Fixed Effects:
        (Intr)
  month -0.772

Membuat model alternatif:

Menyertakan slope dan intercept sebagai efek acak, namun menganggap keduanya tidak berkorelasi.

  > model_lme_02 <- lmer(concDrug ~ month + (month || batch),
  +                      data = dataStabGrouped)
  > summary(model_lme_02)
  Linear mixed model fit by REML ['lmerMod']
  Formula: concDrug ~ month + ((1 | batch) + (0 + month | batch))
     Data: dataStabGrouped

  REML criterion at convergence: 25.9

  Scaled residuals: 
       Min       1Q   Median       3Q      Max 
  -1.61444 -0.28110 -0.01999  0.31604  1.97477 

  Random effects:
   Groups   Name        Variance  Std.Dev. 
   batch    (Intercept) 4.583e-08 0.0002141
   batch.1  month       5.515e-03 0.0742602
   Residual             2.302e-01 0.4798174
  Number of obs: 14, groups:  batch, 5

  Fixed effects:
               Estimate Std. Error t value
  (Intercept) 100.04464    0.19718 507.370
  month        -0.15770    0.06458  -2.442

  Correlation of Fixed Effects:
        (Intr)
  month -0.643

Membuat model ketiga:

Hanya menyertakan intercept sebagai efek acak.

  > # fitting with linear mixed effect model  with random intercepts only and uncorrelated
  > model_lme_03 <- lmer(concDrug ~ month + (1 | batch),
  +                      data = dataStabGrouped)
  > summary(model_lme_03)
  Linear mixed model fit by REML ['lmerMod']
  Formula: concDrug ~ month + (1 | batch)
     Data: dataStabGrouped

  REML criterion at convergence: 26.6

  Scaled residuals: 
       Min       1Q   Median       3Q      Max 
  -1.90012 -0.23099 -0.01921  0.40027  1.79832 

  Random effects:
   Groups   Name        Variance Std.Dev.
   batch    (Intercept) 0.02309  0.1519  
   Residual             0.27846  0.5277  
  Number of obs: 14, groups:  batch, 5

  Fixed effects:
               Estimate Std. Error t value
  (Intercept) 100.04895    0.22677 441.193
  month        -0.16202    0.05908  -2.742

  Correlation of Fixed Effects:
        (Intr)
  month -0.722

Membandingkan ketiga model:

Memilih model terbaik berdasarkan kriteria statistik (seperti AIC atau BIC).

  > anova(model_lme_01, model_lme_02, model_lme_03)
  refitting model(s) with ML (instead of REML)
  Data: dataStabGrouped
  Models:
  model_lme_03: concDrug ~ month + (1 | batch)
  model_lme_02: concDrug ~ month + ((1 | batch) + (0 + month | batch))
  model_lme_01: concDrug ~ month + (month | batch)
               npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)
  model_lme_03    4 28.693 31.249 -10.3463   20.692                     
  model_lme_02    5 30.106 33.301 -10.0528   20.106 0.5869  1     0.4436
  model_lme_01    6 31.926 35.760  -9.9628   19.926 0.1801  1     0.6713

berdasarkan nilai probabilitas (p-value) dari uji Chi-square, tidak ditemukan perbedaan yang signifikan antara ketiga model tersebut (p>0,05). Pada kasus seperti ini, dipilih model yang paling sederhana yaitu model_lme_03.

Memeriksa plot diagnostik:
- Memastikan model sesuai dengan asumsi statistik (misal: normalitas residu, homoskedastisitas).
- plot residual terstandard vs hasil prediksi (fitted value)
```
  > # standardized residuals versus fitted values by batch
  > plot(model_lme_03, resid(., scaled=TRUE) ~ fitted(.) | batch, abline = 0)
```
  pada plot ini terlihat bahwa residual terstandard tidak berada di luar rentang -1,96 sampai +1,96 yang berarti bahwa residual terdistribusi normal. selain itu residual berada di sekitar nol dan tersebar secara acak yang menunjukkan tidak adanya masalah heteroskedastisitas.
- plot hasil prediksi vs eksperimental
```
  > #observed versus fitted values 
  > plot(model_lme_03, concDrug ~ fitted(.) , abline = c(0,1))
  > plot(model_lme_03, concDrug ~ fitted(.) | batch, abline = c(0,1))
```
  Dari plot antara konsentrasi terprediksi dan data eksperimental, diketahui bahwa keduanya memiliki korelasi linier yang baik. Dengan demikian, model_lme_03 dapat memodelkan hubungan antara konsentrasi vs waktu penyimpanan pada lima bets berdasarkan data percobaan.
- hubungan antara konsentrasi tersisa vs waktu penyimpanan pada lima bets obat dapat divisualisasikan sebagai berikut:
```
  > predicted_values<- data_grid(dataStabGrouped, month, batch) %>% 
  +                    add_predictions(model_lme_03)
  > predicted_values %>% 
  +   ggplot(aes(month, pred, color = batch)) +
  +   geom_line() +
  +   geom_point(data = dataStabGrouped, aes(month, concDrug, color = batch)) +
  +   facet_wrap(~ batch) +  # Ini memisahkan per batch
  +   labs(title = "Prediction per batch",
  +        x = "month",
  +        y = "concentration (%)") +
  +   theme_minimal()
```

Estimasi waktu kedaluarsa berdasarkan efek tetap:

Menggunakan efek tetap dari parameter model (fixed-effect)

  > # Ambil koefisien tetap (fixed effects)
  > coef_fix <- fixed.effects(model_lme_03)
  > intercept <- coef_fix[[1]]
  > slope <- coef_fix[[2]]
  > # Tentukan batas minimal kadar obat (misalnya 90%)
  > threshold <- 90
  > # Hitung waktu kedaluwarsa (shelf life)
  > shelf_life <- (threshold - intercept) / slope
  > shelf_life
  [1] 62.02263

berdasarkan fixed-effect maka waktu kedaluarsa adalah sekitar 62 bulan.

Estimasi waktu kedaluarsa berdasarkan simulasi Monte Carlo pada parameter efek tetap (fixed-effect)
- Dalam tulisan ini, pembahasan dibatasi pada estimasi shelf-life berdasarkan confidence interval yang diperoleh dari simulasi Monte Carlo terhadap fixed-effect. Untuk estimasi dengan random-effect (berarti tiap bets akan memiliki hasil estimasi yang berbeda), nanti bisa dibuatkan jika banyak yang request.
- pertama, dibuat dulu simulasi 1000 data virtual
```
  ## membuat simulasi 1000 data virtual
  set.seed(123)
  sim_model <- sim(model_lme_03, n.sims = 1000)
```
- kedua, ekstrak data slope dan intercept hasil simulasi dari fixed-effect
```
  intercepts <- sim_model@fixef[,1]
  slopes <- sim_model@fixef[,2]
```
- ketiga, hitung t90, dan gunakan fungsi pmax( , 0) untuk mengubah hasil estimasi negatif menjadi nol, kemudian periksa histogram distribusi frekuensi datanya
```
  > shelf_life_sim <- (threshold - intercepts) / slopes
  > shelf_life_sim <- shelf_life_sim[shelf_life_sim > 0]
  > hist(shelf_life_sim)
```
```
  > hist(log(shelf_life_sim))
```
- berdasarkan histogram ini, kita tahu bahwa data t90 terdistribusi log-normal dan memiliki outlier.
- keempat, hilangkan outlier dari data logaritma natural t90
```
  > Q1 <- quantile(log(shelf_life_sim), .25)
  > Q3 <- quantile(log(shelf_life_sim), .75)
  > IQR <- IQR(log(shelf_life_sim))
  > no_outliers <- subset(log(shelf_life_sim), 
  +                       log(shelf_life_sim) > (Q1 - 1.5*IQR) & 
  +                         log(shelf_life_sim) < (Q3 + 1.5*IQR))
  > hist(no_outliers)
```
- kelima, hitung nilai rata-rata dan confidence interval dari data yang ditransformasikan, kemudian kembalikan ke skala linier dengan fungsi eksponensial atau exp()
```
  > # Hitung mean dan CI 95%
  > log_ci <- quantile(no_outliers, probs = c(0.025, 0.975))
  > log_ci
      2.5%    97.5% 
  3.550765 4.914701 
  > exp(log_ci)
       2.5%     97.5% 
   34.83996 136.27861 
  > mean.t90 <- mean(no_outliers)
  > exp(mean.t90)
  [1] 63.29042
```
- Rata-rata waktu kedaluarsa yang diperoleh adalah sekitar 63 bulan dengan batas bawah ~34 bulan dan batas atas ~136 bulan. Dalam konteks penjaminan mutu untuk keamanan pasien, dipilih batas bawah sebagai waktu kedaluarsa yaitu 34 bulan.
- Dengan pendekatan ini, hasil estimasi waktu kedaluarsa akan lebih robust, realistis, dan sesuai dengan standar industri farmasi.

Cara Menentukan Kedaluarsa Obat dengan Model Efek Campuran Linier