vector components in fitode?

[bbolker]

Jonathan Dushoff has a current problem where he is experimenting with models using a version of the linear chain trick. He therefore wants to be able to specify models of the form (roughly)

$$ \begin{split} S' & = -\beta S I_{\rm tot} \\ I_1 & = \beta S I_{\rm tot} - n \gamma I_1 \\ I_i & = n \gamma (I_{i-1} - I_{i}), \quad i = 2 \ldots n \\ R & = n \gamma I_n \end{split} $$

(In fact in his case the $\gamma$ values form a geometric progression rather than being constant across boxes)

Do you think there might be a way to accomplish this reasonably straightforwardly?

[parksw3]

I talked to Jonathan about it, but I don't think there's a straightforward way to do this... there are two possibilities.

1. Write a separate function that will turn this kind of model to a fitode model object appropriately. We might need to come up with a reasonable syntax and a few other examples to test out. Doing so for this specific model won't be too difficult, I imagine (I think Jonathan said he'll try this for his model).

2. It's already possible to pass on custom gradient functions though that loses the automatic differentiation capability. Jonathan doesn't seem to like this option too much.

[bbolker]

Here's a simple code generator that should generate something like the model JD wants: I don't know if it's worth the effort in this case or not. (I also don't know if the geometric series of rates is parameterized correctly - it's possible that it would be easier to fit with a simpler parameterization and then transform back to whatever is easier for humans to understand)

n <- 4
Ivec <- paste0("I",1:n)
incidence <- sprintf("beta*S*(%s)", paste(Ivec, collapse = "+"))
Itrans <- function(i) sprintf("n*gamma^%d*I%d", i, i)
vars <- c("S", Ivec, "R")
resp <- c(sprintf("-%s", incidence),
          sprintf("%s - n*gamma*I1", incidence),
          sprintf("%s - %s", Itrans((2:n)-1), Itrans(2:n)),
          Itrans(n))
                  

## could also use base-R Map(), or a for loop ...
purrr::map2(vars, resp, ~reformulate(.x, response = .y))

[parksw3]

This is very neat... you should share with JD if you haven't already.

At the same time, I wonder if it's not worth the trouble writing a super generic function to do this given so many potential model structures... it seems like it would be so much easier to work on special cases like this than to try to figure out one function that covers all possible model structures... something to think about for a longer term

[bbolker]

Agreed (I did share it with him). I wouldn't try to write something general until we had collected a lot of specific examples.

[dushoff]

This is super-cool; I'll talk to Ningrui about trying it out.