search

kingaa / pomp

R package for statistical inference using partially observed Markov processes
https://kingaa.github.io/pomp
GNU General Public License v3.0
111 stars 27 forks source link

Event-based annual aging in an age-structured ODE #89

Closed luccoffeng closed 5 years ago

luccoffeng commented 5 years ago

Hi, I'm using pomp2 and am interested in implementing an age-structured version of a stochastic SIR model. I have checked out an earlier question about this here and the pomp FAQ that it refers to (which covers state vectorisation and matricisation - if that's a word). Still figuring out that FAQ, but I'm pretty sure it doesn't tell me one essential thing that I need for my model: actual aging of the population.

The aging needs to happen as an annual event, N.B. not as a continuous flow between age bins (which would introduce exponential "creeping"). I tried to of implementing this as part of the covariate table, but that only allows defined parameters to vary over time in a step-wise or linear interpolation. And one needs to assume that all covariates change over time either step-wise or linearly, no mix allowed.

Even nicer would be to be able to define events, like the "events" functionality in the deSolve package. From what I can find, there is no such functionality in the pomp or pomp2 package. Is this on the list of things to develop? Or is it already there and have I missed it somehow?

Thanks for the awesome package!

kingaa commented 5 years ago

Have you tried accomplishing this via a time-dependent rprocess function? Since you're not in control of at what times such a function would be evaluated, you'd have to be a little clever in making sure that it is valid for all times, but I can think of several ways of accomplishing this.

With respect to the 'events' functionality, I'm curious as to the purpose you envision for it. Also, have you tried simply passing the argument 'events' to 'trajectory' or 'flow'? Additional arguments to those functions are simply passed down to 'deSolve::ode' when the deterministic skeleton is a vectorfield.

luccoffeng commented 5 years ago

Thanks for the prompt response!

Have you tried accomplishing this via a time-dependent rprocess function? Since you're not in control of at what times such a function would be evaluated, you'd have to be a little clever in making sure that it is valid for all times, but I can think of several ways of accomplishing this.

Nope, but I'm using rprocess = euler(step.fun = Csnippet("..."), delta.t = 1/365) to simulate fixed steps of one day, so I guess it should be possible to build in some sort of IF statement in the Csnippet to check whether current time t is an integer, and then if it is age everybody by one year by moving the bins up and adding newborns. Given my limited experience working in C++, I don't know how sensitive this IF statement will be to rounding errors. You mentioned there are other ways as well?

With respect to the 'events' functionality, I'm curious as to the purpose you envision for it.

I was thinking of planning an annual birthday event during which the whole population ages by one year by calling a function. Not uncommon to do in age-structured SIR modelling. The "forcing" functionality of deSolve::ode is essentially captured by the covariates argument in pomp. It's the "event" functionality that now needs to be hardcoded (as you suggest), which depends on the rprocess to be exactly evaluated at the planned time of event.

Also, have you tried simply passing the argument 'events' to 'trajectory' or 'flow'? Additional arguments to those functions are simply passed down to 'deSolve::ode' when the deterministic skeleton is a vectorfield.

No I haven't as the model I envisage is stochastic. Agreed, it should work for a deterministic version!

kingaa commented 5 years ago

Nope, but I'm using rprocess = euler(step.fun = Csnippet("..."), delta.t = 1/365) to simulate fixed steps of one day, so I guess it should be possible to build in some sort of IF statement in the Csnippet to check whether current time t is an integer, and then if it is age everybody by one year by moving the bins up and adding newborns. Given my limited experience working in C++, I don't know how sensitive this IF statement will be to rounding errors. You mentioned there are other ways as well?

That's right. And you're also right to be worrying about how to make sure that the 'if' statement triggers only at the right times. C arithmetic is just the native arithmetic on the computer, so roundoff error is a potential problem. One thing I've done in this situation is to create an additional state variable that tracks whether the aging has been done recently, to make sure that you don't accidentally age everyone twice. Keep in mind that 'euler' will guarantee that t is advanced by no more than delta.t, but it might be advanced less.

Interesting about the 'events' functionality. Let me think about it a little. Can you remind me about this in a few weeks' time?

luccoffeng commented 5 years ago

Interesting about the 'events' functionality. Let me think about it a little. Can you remind me about this in a few weeks' time?

Will do. And another motivation for having event functionality: when modelling real-life events (e.g. vaccination or a mass drug administration), the event should instantaneously change the state of the system by moving (parts of) the population to other states.

luccoffeng commented 5 years ago

Keep in mind that 'euler' will guarantee that t is advanced by no more than delta.t, but it might be advanced less.

In that case, if delta.t is available within the rprocess (don't know whether it is, would have to check), an IF statement of the sort "if(floor(t) - floor(t - delta.t) > 0)" should do the trick, as long as delta.t is set to be equal to or smaller than 1.

kingaa commented 5 years ago

Yes, delta.t will be available. In the C snippet, it's 'dt'.

To be extra safe, you might do if (floor(t) - floor(t-dt) > 0.5)

kingaa commented 5 years ago

Going to close this now. Feel free to reopen if more discussion is warranted.

luccoffeng commented 5 years ago

Hi Aaron, a little later than intended: end of April you asked me to remind you in a few weeks about the following question I had.

Could the "events" functionality that dsolve has be implemented in pomp2? It would allow the user to specificy times at which user-defined events abruptly change the states in the system. Examples of applications in epidemiological models are annual aging of the population (i.e. moving individuals/populations 1 position ahead in an age-stratified vector), vaccination (i.e. moving the proportion of a population to a different vaccinated state), screening of a population for disease (i.e. moving a proportion of the population to another disease/diagnostic/treatment stage). Events could occur at regular intervals and/or user-defined moments.

So I'm thinking of functionality similar to the covariate_table(), which allows the user to let model parameters vary over time. The 'only' difference is that events would require states within the model to abruptly change at particular points in time.

luccoffeng commented 4 years ago

Hi @kingaa Not sure if the above message with the reminder that you asked for came through. I'm still very much interested in functionality to program events at specific times into pomp models. Recap: by event I mean some event by which the state of the simulation changes instantaneously at a specified moment in time. Example applications would be:

I imagine that this functionality could also be of interest to other fields where events instantaneously disrupt or change the state of the simulation.

Best,

Luc