# Visualisasi data stabilitas obat

Hai semua, dalam tulisan ini saya akan melanjutkan analisis data stabilitas obat dari artikel [Cara Menentukan Tanggal Kedaluwarsa Obat dengan Regresi Linier dan Monte Carlo](https://hashnode.com/edit/cm9zuapyk000709l89qfd09ba).

Data stabilitas dalam artikel tersebut akan lebih mudah dipahami bila divisualisasikan dalam bentuk rata-rata +- standar deviasi. XY-plot yang dihasilkan lebih baik seperti ini:

![](https://cdn.hashnode.com/res/hashnode/image/upload/v1746241897307/16490551-824f-4aa6-b76b-e79247ad2ceb.png align="center")

Untuk memperoleh plot ini, jalankan perintah dengan pemrograman R berikut ini secara berurutan.

```plaintext
> 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
> 
> # memanggil paket yang diperlukan
> library(dplyr)
> library(MASS)
> 
> 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
> 
> # Menghitung rata-rata dan standar deviasi untuk setiap nilai x
> summary_data <- aggregate(conc ~ time, data = data2, FUN = function(z) c(mean = mean(z), sd = sd(z),  n = length(z)))
> 
> # Memisahkan kolom rata-rata dan standar deviasi
> summary_data$mean_conc <- summary_data$conc[, "mean"]
> summary_data$sd_conc <- summary_data$conc[, "sd"]
> summary_data$n_conc <- summary_data$conc[, "n"]
> summary_data$se_conc <- summary_data$sd_conc / sqrt(summary_data$n_conc)
> 
> confidence_level <- 0.95
> alpha <- 1 - confidence_level
> t_critical <- qt(1 - alpha/2, df = summary_data$n_conc - 1)
> summary_data$ci_lower <- summary_data$mean_conc - t_critical * summary_data$se_conc
> summary_data$ci_upper <- summary_data$mean_conc + t_critical * summary_data$se_conc
> # Melakukan regresi linier pada data individual (bukan rata-rata)
> lm1 <- lm(conc ~ time, data = data2)
> newx <- data.frame(time=seq(0,30, length.out=18))
> newx
        time
1   0.000000
2   1.764706
3   3.529412
4   5.294118
5   7.058824
6   8.823529
7  10.588235
8  12.352941
9  14.117647
10 15.882353
11 17.647059
12 19.411765
13 21.176471
14 22.941176
15 24.705882
16 26.470588
17 28.235294
18 30.000000
> 
> newy <- predict(lm1, interval = "confidence", level=.95, newdata = newx)
> data3 <- as.data.frame(cbind(newx, newy[,-1]))
> y.lwr <- aggregate(lwr ~ time, data = data3, FUN = function(z) c(mean = mean(z), sd = sd(z),  n = length(z)))
> y.upr <- aggregate(upr ~ time, data = data3, FUN = function(z) c(mean = mean(z), sd = sd(z),  n = length(z)))
> 
> # Membuat plot menggunakan ggplot2
> library(ggplot2)
> 
> plot <- ggplot(summary_data, aes(x = time, y = mean_conc)) +
+   geom_point(size = 3) +  # Menampilkan titik data rata-rata
+   geom_errorbar(aes(ymin = mean_conc - sd_conc, ymax = mean_conc + sd_conc)) + # Menambahkan error bar (standar deviasi)
+   geom_abline(intercept = coef(lm1)[1], slope = coef(lm1)[2], linetype = "dashed", color = "blue") + # Menambahkan garis regresi
+   geom_line(data = y.lwr, aes(x = time, y = lwr[,1]), linetype = "dotted", color = "red")+
+   geom_line(data = y.upr, aes(x = time, y = upr[,1]), linetype = "dotted", color = "red")+
+   geom_hline(yintercept=90, linetype='dotted', color='gray')+
+   labs(
+     title = "",
+     x = "Time (month)",
+     y = "%-Remaining"
+   ) +
+   theme_bw()+
+   coord_cartesian(xlim = c(0, 30), ylim = c(80, 105))
> plot
```