User interaction with Stan code

[athowes]

The #70 PR raised some issues about how we want the user to be able to interact with the Stan code passed in for specific models (in this case ltcad). This will be an ongoing issue as more models are added.

Currently:

• epidist_stancode method available and exposed
• Takes as argument at least the data, perhaps also family and eventually priors

Questions:

• Extent to which we show and explain the custom Stan code
• Extent to which we expect users to make alterations to custom Stan code

Moving forward:

• I suggest a first thing will be to document via the returns of epidist_stancode.epidist_ltcad what the different objects are and how they change the model. This would then allow the user to (perhaps) alter the output of that Stan code
• Think / discuss how we want this to be handled in future (sorry for vague)

[athowes]

As an example of bad thing due to current implementation:

fit_gamma <- epidist(
  data = prep_obs,
  formula = formula,
  family = epidist_family(prep_obs, family = "gamma"),
  priors = priors,
  stancode = epidist_stancode(
    prep_obs,
    family = epidist_family(prep_obs, family = "gamma")
  )
)

I think if you just did:

stancode = epidist_stancode(prep_obs)

then it would put in the default as epidist_stancode(prep_obs, epidist_family(prep_obs) which would give you the lognormal family. Perhaps easy fix to this with default arguments another way.

[seabbs]

For this specific example can't you pass family in the same way you can pass data

[seabbs]

I.e epidist_standata(data, family)

[athowes]

Issue here about how close to put the family related code together. Prefer to have it all in function but then that means having a combination of the family object and the Stan object together. Might be preferable (if) when we also have the family as a class

[athowes]

I propose to:

• Move epidist_stancode methods into unexported functions
• Move calling of epidist_stancode inside of epidist model fitting
• Move any parts of the epidist_stancode (here so far the family related parts) into the relevant epidist_ functions
  • In future perhaps this could also be the prior part

Thoughts on this proposal @seabbs @kgostic @parksw3? (I will bring this up next meeting [today] and record here anyway.)

My reasons for preferring this:

• Bugs like Writing family insufficient to pass family into epidist #110
• I don't think users are going to be modifying bits of the Stan code in epidist_stancode as stands
• Extensible to modifying other parts of the Stan code
• I think a cleaner pattern for us to have for other models

Uncertainties:

• Which output format favoured for epidist_family now? A named list? family$family and family$stancode?

[athowes]

Third option:

• Keep all Stan code together but still put it into epidist
• Called internally
• Rough agreement to move forward with this for now

Obvious option to change:

• Prior on primary event changed based on data

[athowes]

I have started working on this in #115.

So far what I have done is move stancode <- epidist_stancode(data = data, family = family) inside of epidist.default rather than have it as an argument. This fixes the bug I had (#110) as follows:

> meanlog <- 1.8
> sdlog <- 0.5
> obs_time <- 25
> sample_size <- 200
> 
> sim_obs <- simulate_gillespie() |>
+   simulate_secondary(
+     meanlog = meanlog,
+     sdlog = sdlog
+   ) |>
+   observe_process() |>
+   filter_obs_by_obs_time(obs_time = obs_time) %>%
+   .[sample(seq_len(.N), sample_size, replace = FALSE)]
> prep_obs <- epidist_prepare(sim_obs, model = "latent_individual")
> fit_gamma <- epidist(data = prep_obs, family = epidist_family(prep_obs, family = "gamma"))
Compiling Stan program...
Start sampling
Running MCMC with 4 sequential chains...

Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 1 Iteration:  100 / 2000 [  5%]  (Warmup) 
Chain 1 Iteration:  200 / 2000 [ 10%]  (Warmup) 
Chain 1 Iteration:  300 / 2000 [ 15%]  (Warmup) 
Chain 1 Iteration:  400 / 2000 [ 20%]  (Warmup) 

This matches the specification under "third option" in comment #78 (comment) but I don't think is the complete solution we had in mind.

It's not clear to me the best next step now. Options as I see them:

1. Leave it as it is

• Note here: is it possible to put another argument stancode_overwrite = NULL such that if it's used then the user could put their own custom Stan code in?

2. Move the "family" related function into the epidist_family.latent_individual function which then outputs some other object (I suggested a list in User interaction with Stan code #78 (comment) and don't recall any alternative proposals).

• epidist.default would then contain code to reassemble the Stan code from the family-specific part and other part
• The family argument could be removed from epidist_stancode (this is good because currently it's bad to pass it)
• This raises point about how "default" the epidist method is now becoming

Tagging @seabbs to give input on this before I move forward with implementing.

[seabbs]

Nice issue.

but I don't think is the complete solution we had in mind.

Could you flesh this out a bit please?

[athowes]

Yep, as to why I had that impression:

1. I recall an interest in keeping the family function code with the family object (this solution does not do that)
2. Perhaps something I don't like about the current use of stancode <- epidist_stancode(data = data, family = family). Specifically I am unsure how this will work within epidist.default once we have other models. But perhaps it's fine.

If these do not seem important issues I would be happy to move forward with the solution as written in the PR currently, and then we can revisit it when we add more models.

I would also open an issue for "add custom Stan code prior" would would entail:

• Add ability to pass in custom Stan prior to epidist_stancode.latent_individual for additional Stan parameters
• Alter epidist_stancode within epidist to take prior as an argument
  • But then the blocker here is that epidist_prior should be edited to contain brmsprior and custom Stan code or something like that.. so it's unclear how well it works

[athowes]

@seabbs ping here -- could you give thoughts before I either look to get a "leave as is" solution merged, or try implementing the "move the family into epidist_family option"?

[seabbs]

I think leave as is (the 2. Point is handled by S3 dispatch)

And make issues for your point about custom priors and reorg the family function I think?

[athowes]

I have made an issue for the priors point (#123).

I don't fully understand what making an issue for the family function reorganisation would entail but I am going to try.