Priors functions

[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

[ntorresd]

The work done by @zmcucunuba and @sumalibajaj on this issue has been essential to solve the lack of flexibility to specify the the stan models in serofoi. The way we are handling this (as of #200) is by means of functions like sf_normal() and sf_uniform(), that allows to specify the prior distributions to be used for inference (either for the FOI or the seroreversion rate). Setting the necessary stan data parameters is done internally by means of build_stan_data() and set_stan_data_defaults() using the priors and model specifications as input. For instance, to specify an age-varying model with a normal distribution centered around 0 and sd of 1 for the first FOI value (forward random walk), and a seroreversion uniformly distributed between 0 and 2:

seromodel <- fit_seromodel(
  serosurvey = serosurvey,
  model_type = "age",
  foi_prior = sf_normal(0, 1),
  is_seroreversion = TRUE,
  seroreversion_prior = sf_uniform(0, 2)
)

So far, we have only implemented forward random walk methods for age/time varying models, which is why we haven't had the need to use the functions set_priors_rw() and set_priors_iid() to set the models; these functions may come in handy once we decide to implement IID models. We may reopen this issue once we reach that point.