Priors functions

[zmcucunuba]

No description provided.

[zmcucunuba]

For priors, we are working on branch priors from dev
https://github.com/epiverse-trace/serofoi/tree/priors

[sumalibajaj]

library(dplyr)
sf_normal <- function(mean=0, sd=1) {
  # Restricting normal inputs to be non-negative
  if(mean < 0 | sd <= 0){
    print("Normal distribution here only accepts non-negative values for mean and standard deviation")
    stop()}
  
  return (list(mean=mean, sd=sd, name="normal"))
}

sf_uniform <- function(min=0, max=1) {
  # Restricting uniform inputs to be non-negative
  if(min < 0 | max < 0){
    print("Uniform distribution here only accepts non-negative values for min and max")
    stop()}
  if (min >= max) {
    print("Uniform distribution only accepts min < max")
  }
  
    return (list(min=min, max=max, name="uniform"))
}

sf_student_t <- function(mean=0, sd=1, nu=3) {
  # Restricting student_t inputs to be non-negative
  if(mean < 0 | sd < 0 | nu < 0){
    print("Student_t distribution here only accepts non-negative values for mean, sd, and nu")
    stop()}
  return (list(mean=mean, sd=sd, nu = nu, name="student_t"))
}

sf_none <- function() {
  return (list(name="none"))
}


dist_list_to_vector <- function(dist_list) {
  dist_codes <- list(none=0, normal=1, uniform=2, student_t=3)
  dist_list$name <- dist_codes[dist_list$name]
  return (unlist(dist_list, use.names=FALSE))
}

set_priors_rw <- function(my_model,
                          foi_1 = sf_normal(), 
                          foi_i = sf_normal(), 
                          sigma_cauchy_rw = 1,
                          seroreversion_rate=sf_none()) {
  # Restricting foi_1 to be only normal or uniform
  if(!(foi_1$name == "normal"| foi_1$name == "uniform")){
    print("foi_1 only accepts normal or uniform distributions")
    stop()}

  # Restricting FOI_i to be only normal or student_t
  if(!(foi_i$name=="normal"| foi_i$name=="student_t")){
    print("foi_i only accepts normal or student_t distributions for the random walk")
    stop()}
  
  # Restricting sigma_cauchy_rw to be > 0 
  if(sigma_cauchy_rw <= 0){
    print("The standard devation of the cauchy distrubtion can only be non-negative")
    stop()}
  
  # If seroreversion in model is TRUE:
  # and no input from user: set seroreversion_rate prior as normal
  # if there is an input from user: ERROR (if not normal or uniform)
  
  # If seroreversion in model is FALSE:
  # and no input from user: No problem- using none (default)
  # if there is an input from user: ERROR (don't supply prior)
  if(my_model$is_seroreversion == TRUE){
    if(seroreversion_rate$name == "none"){
      print("seroreversion rate can't be empty because the model has seroreversion, setting to default normal(0, 1)")
      seroreversion_rate = sf_normal() # default if is_seroreversion == TRUE
    } else if(!(seroreversion_rate$name %in% c("none", "normal", "uniform"))){
        print("seroreversion rate can only be normal or uniform")
        stop()
      }
  } else{
      if(seroreversion_rate$name != "none"){
        print("model does not have seroreversion, don't add any prior")
        stop()
        }
    }
    
  return (list(foi_1_vector=dist_list_to_vector(foi_1), 
              foi_i_vector=dist_list_to_vector(foi_i), 
              sigma_cauchy_rw=sigma_cauchy_rw, 
              seroreversion_rate_vector=dist_list_to_vector(seroreversion_rate)))
  
}




set_priors_iid <- function(my_model,
                           foi = sf_normal(), 
                           seroreversion_rate = sf_none()) {
  # Restricting foi to be only normal or uniform
  if(!(foi$name=="normal"| foi$name=="uniform")){
    print("foi only accepts normal or uniform distributions")
    stop()
  }
  
  # If seroreversion in model is TRUE:
  # and no input from user: set seroreversion_rate prior as normal
  # if there is an input from user: ERROR (if not normal or uniform)
  
  # If seroreversion in model is FALSE:
  # and no input from user: No problem- using none (default)
  # if there is an input from user: ERROR (don't supply prior)
  if(my_model$is_seroreversion == TRUE){
    if(seroreversion_rate$name == "none"){
      print("seroreversion rate can't be empty because the model has seroreversion, setting to default normal(0, 1)")
      seroreversion_rate = sf_normal() # default if is_seroreversion == TRUE
    } else if(!(seroreversion_rate$name %in% c("none", "normal", "uniform"))){
      print("seroreversion rate can only be normal or uniform")
      stop()
    }
  } else{
    if(seroreversion_rate$name != "none"){
      print("model does not have seroreversion, don't add any prior")
      stop()
    }
  }
  
  
  return (list(foi_vector=dist_list_to_vector(foi), 
               seroreversion_rate_vector = dist_list_to_vector(seroreversion_rate)))
  
}

[zmcucunuba]

Testing the wrapper set_priors function

set_sero_priors <- function(my_model,
                            foi,
                            foi_1,
                            foi_i,
                            seroreversion_rate,
                            sigma_cauchy_rw,
                            type = "RW_forward")
{

  # Restricting the parameters according to IID vs RW option
  # if (type == "IID") {
  #   if(exists("foi_1")| exists("foi_i")| exists("sigma_cauchy_rw"))
  #     print("this is IDD so there is only a foi value expected")
  #   stop()
  # Jaime is gonna solve this
  # }

  return(list_of_priors_to_go_to_stan)
}


set_priors(type="IID",
           foi = sf_normal(),
           foi_1=sf_normal(),
           foi_i=sf_normal(),
           sigma_cauchy_rw= 1,
           seroreversion_rate=sf_none())

set_priors(type="IID",
           foi = sf_normal())

[zmcucunuba]

List of priors in Excel

https://docs.google.com/spreadsheets/d/1tGm_XBoUUwqcrpNt7MyQeUm7y7gg9meKNEPAwb-CLeg/edit#gid=1102747719