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

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