afrezza_GIR.R

library(tidyverse)
library(tidymodels)
library(cowplot)

afrezza_data_GIR <- tibble(
  time = c(
    3.29696260543123, 10.7000146262981, 13.3128565160158, 15.851045780313,
    17.7920140412462, 18.5385402954512, 20.0315928038613, 21.5246453122715,
    22.4578031300278, 23.3909609477841, 24.8840134561942, 28.6166447272195,
    30.6695919262835, 32.9091706888986, 35.4371800497294, 34.9621178879625,
    43.9204329384233, 46.9065379552435, 51.2612744381064, 56.6113792599093,
    61.0905367851397, 64.0766418019599, 66.5295137800623, 70.7953780898055,
    74.5280093608308, 80.2762615182097, 86.6383241512684, 91.1008922041831,
    96.1772707327775, 100.283165130905, 103.642533274828, 109.614743308469,
    117.826532104724, 126.03832090098, 134.250109697236, 142.461898493491,
    150.673687289747, 158.885476086003, 167.097264882258, 175.309053678514,
    183.52084247477, 191.732631271025, 199.944420067281, 208.156208863537,
    216.367997659792, 224.579786456048, 232.791575252304, 241.003364048559,
    249.215152844815, 257.426941641071, 265.638730437326, 273.850519233582,
    282.062308029838, 290.274096826093, 298.485885622349, 306.697674418605,
    314.90946321486, 323.121252011116, 331.333040807372, 339.544829603627,
    347.756618399883, 356.279459802057),
  GIR = c(
    0.492673546520211, 0.747580609442821, 1.06606726311749, 1.44377712817787,
    1.79840995659492, 2.23993489237311, 2.62965423568666, 2.98169678448047,
    3.39496412436885, 3.80469926477102, 4.14732261493489, 4.38868957982697,
    4.75132872707946, 5.18932146337142, 5.51492603419258, 5.80392417397467,
    5.10925827501698, 4.89968110550093, 4.97503469454041, 4.70591473368513,
    4.32887766852688, 3.96718044113739, 3.56148210014806, 3.18727079457879,
    2.94708122951546, 2.58397112233147, 2.59982677335853, 2.54864912746922,
    2.31612948129026, 1.98208432987864, 1.66701213570733, 1.33793958720879,
    1.20821153335105, 1.12600761803526, 0.914075648861725, 0.775356541766321,
    0.693152626450527, 0.582691115244929, 0.56342457259279, 0.43626539108867,
    0.268518026397379, 0.189867717837424, 0.138276198068917, 0.129713290223521,
    0.107021584433225, 0.0988868219800994, 0.0963179496264797, 0.0997431127646378,
    0.0817610062893088, 0.0817610062893088, 0.0817610062893088, 0.0817610062893088,
    0.0817610062893088, 0.0817610062893088, 0.0779076977588806, 0.07919213393569,
    0.0766232615820712, 0.0817610062893088, 0.0817610062893088, 0.0817610062893088,
    0.0817610062893088, 0.0417294121120859))

get_insulin_activity_curve <- function(time, peakActivityTime, actionDuration, delay) {
  time_after_delay <- time - delay
  tau = peakActivityTime * (1 - peakActivityTime / actionDuration) / (1 - 2 * peakActivityTime / actionDuration)
  a = 2 * tau/ actionDuration
  S = 1 / (1 - a + (1 + a) * exp(-actionDuration / tau))
  activity_curve <- (S / tau^2) * time_after_delay * (1 - time_after_delay / actionDuration) * exp(-time_after_delay / tau)
  activity_curve[time_after_delay < 0] <- 0
  activity_curve[time_after_delay > actionDuration] <- 0
  activity_curve
}

function_to_optimize <- function(x, loss_type = "rmse") {
  peakActivityTime <- x[1]
  actionDuration <- x[2]
  delay <- x[3]
  raw_activity_curve <- get_insulin_activity_curve(afrezza_data_GIR[["time"]], peakActivityTime, actionDuration, delay)
  scaled_activity_curve <- raw_activity_curve * max(afrezza_data_GIR$GIR) / max(raw_activity_curve)
  case_when(
    loss_type == "rmse" ~ rmse_vec(afrezza_data_GIR[["GIR"]], scaled_activity_curve),
    loss_type == "msd" ~ msd_vec(afrezza_data_GIR[["GIR"]], scaled_activity_curve),
    loss_type == "mae" ~ mae_vec(afrezza_data_GIR[["GIR"]], scaled_activity_curve),
    loss_type == "mape" ~ mape_vec(afrezza_data_GIR[["GIR"]], scaled_activity_curve),
    loss_type == "smape" ~ smape_vec(afrezza_data_GIR[["GIR"]], scaled_activity_curve),
    loss_type == "mase" ~ mase_vec(afrezza_data_GIR[["GIR"]], scaled_activity_curve))
}


optim_sols <- 
  tibble(loss_type = c("rmse", "msd", "mae", "mape", "smape", "mase")) %>% 
  mutate(optim_sol = map(loss_type,
                         ~optim(par = c(10, 240, 5),
                                fn = function(x) function_to_optimize(x, loss_type = .),
                                lower = c(0, 30, 0),
                                upper = c(60, 360, 20),
                                method = "L-BFGS-B"))) %>% 
  unnest_wider(optim_sol)


# Plot Solutions ----------------------------------------------------------
optim_sols %>% 
  filter(convergence == 0) %>% 
  select(loss_type, par) %>% 
  mutate(activity_curve = 
           map(par, ~{
             raw_activity_curve <- get_insulin_activity_curve(afrezza_data_GIR[["time"]], .[1], .[2], .[3])
             scaled_activity_curve <- raw_activity_curve * max(afrezza_data_GIR$GIR) / max(raw_activity_curve)
             tibble(time = afrezza_data_GIR[["time"]], GIR = scaled_activity_curve)
  })) %>% 
  select(loss_type, activity_curve) %>% 
  unnest(activity_curve) %>% 
  ggplot(aes(time, GIR)) +
  geom_point(data = afrezza_data_GIR) +
  geom_line(mapping = aes(color = loss_type)) +
  scale_color_discrete(name = "Loss Type") +
  scale_x_continuous(name = "Time (Minutes)") +
  theme_minimal_grid()

afrezza_data_GIR %>% 
  mutate(optimized_curve = get_insulin_activity_curve)


optim_sols %>% 
  filter(convergence == 0) %>% 
  select(loss_type, par) %>% 
  mutate(activity_curve = 
           map(par, ~{
             raw_activity_curve <- get_insulin_activity_curve(afrezza_data_GIR[["time"]], .[1], .[2], .[3])
             scaled_activity_curve <- raw_activity_curve * max(afrezza_data_GIR$GIR) / max(raw_activity_curve)
             tibble(time = afrezza_data_GIR[["time"]], GIR = scaled_activity_curve)
           })) %>% 
  select(loss_type, activity_curve) %>% 
  unnest(activity_curve) %>% 
  ggplot(aes(time, GIR)) +
  geom_point(data = afrezza_data_GIR) +
  geom_line(mapping = aes(color = loss_type)) +
  scale_color_discrete(name = "Loss Type") +
  scale_x_continuous(name = "Time (Minutes)") +
  theme_minimal_grid() +
  ggtitle("Several Optimized Curves for Afrezza Insulin Model")


ggplot(mapping = aes(time, GIR)) +
  geom_point(data = afrezza_data_GIR) +
  geom_line(data = tibble(time = afrezza_data_GIR[["time"]],
                          GIR_raw = get_insulin_activity_curve(afrezza_data_GIR[["time"]], 29, 300, 10)) %>% 
              mutate(GIR = GIR_raw * max(afrezza_data_GIR$GIR) / max(GIR_raw))) +
  theme_minimal_grid() +
  scale_x_continuous(name = "Time (Minutes)") +
  ggtitle("Final Exponential Model Parameter Choices")



# Parameter Summary -------------------------------------------------------
optim_sols %>% 
  mutate(par = map(par, ~set_names(., c("peakActivityTime", "actionDuration", "delay")))) %>% 
  select(loss_type, par) %>%
  unnest_wider(par) %>%
  knitr::kable()