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