nmode's Git Repositories - Rnaught/blob - inst/app/scripts/estimators.R

   1 # Main logic block for estimator-related interactions.
   2 estimators_logic <- function(input, output, react_values) {
   3   # Initialize a data frame to hold estimates.
   4   react_values$estimates_table <- data.frame(
   5     Estimator = character(0),
   6     `Serial interval` = character(0),
   7     check.names = FALSE
   8
   9   )
  10   # Initialize a list to hold added estimators.
  11   react_values$estimators <- list()
  12
  13   add_id(input, output, react_values)
  14   add_idea(input, output, react_values)
  15   add_seq_bayes(input, output, react_values)
  16   add_wp(input, output, react_values)
  17
  18   render_estimates(output, react_values)
  19   delete_estimators(input, react_values)
  20   export_estimates(output, react_values)
  21 }
  22
  23 # If an estimator is added, ensure it is not a duplicate and add it to the list
  24 # of estimators. This function should be called at the end of each
  25 # estimator-specific 'add' function, after validating their parameters.
  26 add_estimator <- function(method, new_estimator, output, react_values) {
  27   num_estimators <- length(react_values$estimators)
  28
  29   # Check whether the new estimator is a duplicate, and warn if so.
  30   for (i in seq_len(num_estimators)) {
  31     if (identical(new_estimator, react_values$estimators[[i]])) {
  32       showNotification("Error: This estimator has already been added.",
  33         duration = 3, id = "notify-error"
  34       )
  35       return()
  36     }
  37   }
  38
  39   # Add the new estimator to the list of estimators.
  40   react_values$estimators[[num_estimators + 1]] <- new_estimator
  41
  42   showNotification("Estimator added successfully.",
  43     duration = 3, id = "notify-success"
  44   )
  45
  46   # Evaluate the new estimator on all existing datasets and create a new row in
  47   # the estimates table.
  48   update_estimates_row(new_estimator, react_values)
  49 }
  50
  51 # Ensure serial intervals are specified as positive numbers.
  52 validate_mu <- function(method, input, output) {
  53   mu <- suppressWarnings(as.numeric(trimws(input[[paste0("mu_", method)]])))
  54   if (is.na(mu) || mu <= 0) {
  55     output[[paste0("mu_", method, "_warn")]] <- renderText(
  56       "The serial interval must be a positive number."
  57     )
  58     return(NULL)
  59   }
  60   output[[paste0("mu_", method, "_warn")]] <- renderText("")
  61   mu
  62 }
  63
  64 # Incidence Decay (ID).
  65 add_id <- function(input, output, react_values) {
  66   observeEvent(input$add_id, {
  67     mu <- validate_mu("id", input, output)
  68     if (!is.null(mu)) {
  69       new_estimator <- list(
  70         method = "id", mu = mu, mu_units = input$mu_id_units
  71       )
  72       add_estimator("id", new_estimator, output, react_values)
  73     }
  74   })
  75 }
  76
  77 # Incidence Decay and Exponential Adjustment (IDEA).
  78 add_idea <- function(input, output, react_values) {
  79   observeEvent(input$add_idea, {
  80     mu <- validate_mu("idea", input, output)
  81     if (!is.null(mu)) {
  82       new_estimator <- list(
  83         method = "idea", mu = mu, mu_units = input$mu_idea_units
  84       )
  85       add_estimator("idea", new_estimator, output, react_values)
  86     }
  87   })
  88 }
  89
  90 # Sequential Bayes (seqB).
  91 add_seq_bayes <- function(input, output, react_values) {
  92   observeEvent(input$add_seq_bayes, {
  93     mu <- validate_mu("seq_bayes", input, output)
  94
  95     kappa <- trimws(input$kappa)
  96     kappa <- if (kappa == "") 20 else suppressWarnings(as.numeric(kappa))
  97
  98     if (is.na(kappa) || kappa <= 0) {
  99       output$kappa_warn <- renderText(
 100         "The maximum prior must be a positive number."
 101       )
 102     } else if (!is.null(mu)) {
 103       output$kappa_warn <- renderText("")
 104       new_estimator <- list(
 105         method = "seq_bayes", mu = mu,
 106         mu_units = input$mu_seq_bayes_units, kappa = kappa
 107       )
 108       add_estimator("seq_bayes", new_estimator, output, react_values)
 109     }
 110   })
 111 }
 112
 113 # White and Pagano (WP).
 114 add_wp <- function(input, output, react_values) {
 115   observeEvent(input$add_wp, {
 116     if (input$wp_mu_known == "Yes") {
 117       mu <- validate_mu("wp", input, output)
 118       if (!is.null(mu)) {
 119         new_estimator <- list(method = "wp",
 120           mu = mu, mu_units = input$mu_wp_units
 121         )
 122         add_estimator("wp", new_estimator, output, react_values)
 123       }
 124     } else {
 125       grid_length <- trimws(input$grid_length)
 126       max_shape <- trimws(input$max_shape)
 127       max_scale <- trimws(input$max_scale)
 128
 129       suppressWarnings({
 130         grid_length <- if (grid_length == "") 100 else as.numeric(grid_length)
 131         max_shape <- if (max_shape == "") 10 else as.numeric(max_shape)
 132         max_scale <- if (max_scale == "") 10 else as.numeric(max_scale)
 133       })
 134
 135       valid <- TRUE
 136
 137       if (is.na(grid_length) || grid_length <= 0) {
 138         output$grid_length_warn <- renderText(
 139           "The grid length must be a positive integer."
 140         )
 141         valid <- FALSE
 142       } else {
 143         output$grid_length_warn <- renderText("")
 144       }
 145
 146       if (is.na(max_shape) || max_shape <= 0) {
 147         output$max_shape_warn <- renderText(
 148           "The maximum shape must be a positive number."
 149         )
 150         valid <- FALSE
 151       } else {
 152         output$max_shape_warn <- renderText("")
 153       }
 154
 155       if (is.na(max_scale) || max_scale <= 0) {
 156         output$max_scale_warn <- renderText(
 157           "The maximum scale must be a positive number."
 158         )
 159         valid <- FALSE
 160       } else {
 161         output$max_scale_warn <- renderText("")
 162       }
 163
 164       if (valid) {
 165         new_estimator <- list(method = "wp", mu = NA, grid_length = grid_length,
 166           max_shape = max_shape, max_scale = max_scale
 167         )
 168         add_estimator("wp", new_estimator, output, react_values)
 169       }
 170     }
 171   })
 172 }
 173
 174 # Convert an estimator's specified serial interval to match the time units of
 175 # the given dataset.
 176 convert_mu_units <- function(data_units, estimator_units, mu) {
 177   if (data_units == "Days" && estimator_units == "Weeks") {
 178     return(mu * 7)
 179   } else if (data_units == "Weeks" && estimator_units == "Days") {
 180     return(mu / 7)
 181   }
 182   mu
 183 }
 184
 185 # Add a row to the estimates table when a new estimator is added.
 186 update_estimates_row <- function(estimator, react_values) {
 187   dataset_rows <- seq_len(nrow(react_values$data_table))
 188   estimates <- c()
 189
 190   if (nrow(react_values$data_table) > 0) {
 191     estimates <- dataset_rows
 192     for (row in dataset_rows) {
 193       estimate <- eval_estimator(estimator, react_values$data_table[row, ])
 194       estimates[row] <- estimate
 195     }
 196   }
 197
 198   new_row <- data.frame(
 199     t(c(estimator_name(estimator), estimator_mu_text(estimator), estimates))
 200   )
 201   colnames(new_row) <- colnames(react_values$estimates_table)
 202
 203   react_values$estimates_table <- rbind(
 204     react_values$estimates_table, new_row
 205   )
 206 }
 207
 208 # Evaluate the specified estimator on the given dataset.
 209 eval_estimator <- function(estimator, dataset) {
 210   cases <- as.integer(unlist(strsplit(dataset[, 3], ",")))
 211
 212   if (estimator$method == "id") {
 213     mu <- convert_mu_units(dataset[, 2], estimator$mu_units, estimator$mu)
 214     estimate <- round(Rnaught::id(cases, mu), 2)
 215   } else if (estimator$method == "idea") {
 216     mu <- convert_mu_units(dataset[, 2], estimator$mu_units, estimator$mu)
 217     estimate <- round(Rnaught::idea(cases, mu), 2)
 218   } else if (estimator$method == "seq_bayes") {
 219     mu <- convert_mu_units(dataset[, 2], estimator$mu_units, estimator$mu)
 220     estimate <- round(Rnaught::seq_bayes(cases, mu, estimator$kappa), 2)
 221   } else if (estimator$method == "wp") {
 222     if (is.na(estimator$mu)) {
 223       estimate <- Rnaught::wp(cases, serial = TRUE,
 224         grid_length = estimator$grid_length,
 225         max_shape = estimator$max_shape, max_scale = estimator$max_scale
 226       )
 227       estimated_mu <- round(sum(estimate$supp * estimate$pmf), 2)
 228       estimate <- paste0(round(estimate$r0, 2), " (SI = ", estimated_mu,
 229         " ", tolower(dataset[, 2]), ")"
 230       )
 231     } else {
 232       mu <- convert_mu_units(dataset[, 2], estimator$mu_units, estimator$mu)
 233       estimate <- round(Rnaught::wp(cases, mu), 2)
 234     }
 235   }
 236
 237   return(estimate)
 238 }
 239
 240 # Create the name of an estimator to be added to the first column of the
 241 # estimates table.
 242 estimator_name <- function(estimator) {
 243   if (estimator$method == "id") {
 244     return("ID")
 245   } else if (estimator$method == "idea") {
 246     return("IDEA")
 247   } else if (estimator$method == "seq_bayes") {
 248     return(paste0("seqB", " (&#954; = ", estimator$kappa, ")"))
 249   } else if (estimator$method == "wp") {
 250     if (is.na(estimator$mu)) {
 251       return(paste0("WP (", estimator$grid_length, ", ",
 252         round(estimator$max_shape, 3), ", ", round(estimator$max_scale, 3), ")"
 253       ))
 254     } else {
 255       return("WP")
 256     }
 257   }
 258 }
 259
 260 # Create the text to be displayed for the serial interval in the second column
 261 # of the estimates table.
 262 estimator_mu_text <- function(estimator) {
 263   if (is.na(estimator$mu)) {
 264     return("—")
 265   }
 266   paste(estimator$mu, tolower(estimator$mu_units))
 267 }
 268
 269 # Render the estimates table whenever it is updated.
 270 render_estimates <- function(output, react_values) {
 271   observe({
 272     output$estimates_table <- DT::renderDataTable(react_values$estimates_table,
 273       escape = FALSE, rownames = FALSE,
 274       options = list(
 275         columnDefs = list(list(className = "dt-left", targets = "_all"))
 276       ),
 277     )
 278   })
 279 }
 280
 281 # Delete rows from the estimates table and the corresponding estimators.
 282 delete_estimators <- function(input, react_values) {
 283   observeEvent(input$estimators_delete, {
 284     rows_selected <- input$estimates_table_rows_selected
 285     react_values$estimators <- react_values$estimators[-rows_selected]
 286     react_values$estimates_table <-
 287       react_values$estimates_table[-rows_selected, ]
 288   })
 289 }
 290
 291 # Export estimates table as a CSV file.
 292 export_estimates <- function(output, react_values) {
 293   output$estimates_export <- downloadHandler(
 294     filename = function() {
 295       paste0(
 296         "Rnaught_estimates_", format(Sys.time(), "%y-%m-%d_%H-%M-%S"), ".csv"
 297       )
 298     },
 299     content = function(file) {
 300       output_table <- data.frame(
 301         lapply(react_values$estimates_table, sub_entity)
 302       )
 303       colnames(output_table) <- sub_entity(
 304         colnames(react_values$estimates_table)
 305       )
 306       write.csv(output_table, file, row.names = FALSE)
 307     }
 308   )
 309 }
 310
 311 # Substitute HTML entity codes with natural names.
 312 sub_entity <- function(obj) {
 313   obj <- gsub("&#954;", "kappa", obj)
 314   obj
 315 }