Restructuring inference methods

[ben18785]

First, we set up the data:

data(chagas2012)
serodata <- prepare_serodata(chagas2012)

Then we currently fit the model to the data:

seromodel_tv_normal <- fit_seromodel(
  serodata = serodata,
  foi_model = "tv_normal",
  foi_parameters = list(
    foi_location = 0.5,
    foi_scale = 0.1
  ),
  chunks = rep(c(1, 2, 3), c(23, 10, 15)),
  iter = 1500
)

Our proposed changes are:

Time-varying:

# create foi indices
years <- seq(1990, 2000, 1)
chunks <- c(rep(1, 5), rep(2, 6))
chunks_df <- tibble(year=years, foi_index=chunks)

• Priors functions #193

# specify prior
prior_example <- set_prior(type=”RW”, foi_0 = “normal(0,1)”, sigma=”cauchy(0, 2)”)

# specify likelihood
model_example <- set_model(type=”time”, seroreversion=FALSE, foi_indices=chunks_df)

# internal function
example_data_for_stan <- prepare_stan_data_time_model(
  data=serodata,
  model=model_example,
  prior=prior_example
)

# fit model
fit <- fit_seromodel(data=serodata, model=model_example, prior=prior_example)

Age-varying:

# create foi indices
ages <- seq(1, 11, 1)
chunks <- c(rep(1, 5), rep(2, 6))
chunks_df <- tibble(age=ages, foi_index=chunks)

# specify prior
prior_example <- set_prior(type=”IID”, foi = “uniform(0,100)”)

# specify likelihood
model_example <- set_model(type=”age”, seroreversion=TRUE, foi_indices=chunks_df)

# fit model
fit <- fit_seromodel(data=serodata, model=model_example, prior=prior_example)

Time- and age-varying:

# age
ages <- seq(1, 11, 1)
chunks <- c(rep(1, 5), rep(2, 6))
age_df <- tibble(age=ages, foi_index=chunks)

# time
years <- seq(1990, 2000, 1)
chunks <- c(rep(1, 5), rep(2, 6))
time_df <- tibble(year=years, foi_index=chunks)

# create list of foi index dfs
chunks_list <- list(
  age=age_df,
  time=time_df
)

# specify prior
prior_example <- set_prior(type=”IID”, foi = “uniform(0,100)”)

# specify likelihood
model_example <- set_model(type=”age-time”, seroreversion=TRUE, foi_indices=chunks_list)

# fit model
fit <- fit_seromodel(data=serodata, model=model_example, prior=prior_example)

[zmcucunuba]

to develop this idea

[zmcucunuba]

"Final" structure to the wrapper functions

#----- PREPARE SERODATA

prepare_sero_data <- function(raw_data) {

stan_sero_data <- data.frame(n_tested, n_seropositive, year_survey, age_min, age_max)
return(stan_sero_data)

}

my_data <- prepare_sero_data(my_raw_data)

#----- SET SEROMODEL

set_sero_model <- function(model = "time",
seroreversion = TRUE,
sf_data = my_data,
foi_windows = my_foi_windows) {

check_set_foi_windows <- function (sf_data, foi_windows){
#This fucntion will
# - check the foi_windows and corrects it when possible
# otherwise returns an error
return(final_foi_windows)
}

if (!exists(foi_windows)) {

# use sf_data and build a vector with consecutive birth_cohort
final_foi_windows <- data.frame (year = birth_cohort, window: 1:length(birth_cohort))
# we need to add also if model == "age" and if model == "age-time"

} else (check_set_foi_windows(foi_windows))

return(list(model = model,
is_foi_log = is_foi_log,
seroreversion = seroreversion,
sf_data = sf_data,
foi_windows = final_foi_windows))

}

my_model <- set_sero_model(model = "time",
is_foi_log = TRUE,
seroreversion = TRUE,
sf_data = my_data,
foi_windows = my_foi_windows)

set_sero_priors <- function(my_model,
is_foi_log = TRUE,
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)
}

my_priors <- set_sero_priors(my_model,
is_foi_log = TRUE,
foi,
foi_1,
foi_i,
seroreversion_rate,
sigma_cauchy_rw,
type = "RW_forward")

#----- FIT SEROMODEL

fit_sero_model <- function(sf_model = my_model,
sf_data = my_data,
sf_priors = my_priors,
...) {

This function does:

1. Select the stanfile

2. Compile the stanfile

3. Run the model

return (sero_stan_file)

}

my_fit <- fit_sero_model(sf_model = my_model,
sf_data = my_data,
sf_priors = my_priors)