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.

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:

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

> 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