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)

• [x] #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)