Inference of parameter and non-parametric function for ODE model

Machine Learning and Modeling
1

Hi all,

I'm working on an ODE model where I need to infer a parameter μ and a non-parametric function φ(t), and I'm looking for recommendations for tools or approaches in R to solve this problem. My dataset is small, involving just 7 data points in the form (t, n).

The model is given by:

dn/dt = φ(t) - μ * n(t), n(0) = 0

It has an analytical solution in integral form:

n(t) = ∫0t exp(-μ (t-s)) φ(s) ds

I've been exploring the fungp and magi R packages, as they both use Gaussian processes to model φ(t). However, I’ve encountered a couple of challenges:

  • fungp doesn’t seem to support parameter fitting, which is crucial for estimating μ.
  • magi doesn't allow for parameters that vary over time (which is the case for φ(t)).

Can anyone suggest other R packages or approaches that might be a better fit for this problem? I’m particularly interested in methods that can handle both the parameter fitting for μ and the time-varying nature of φ(t).

Thanks in advance for your help!

2

Hi,

I am not exactly sure what you mean by non-parametric function φ(t)... But assuming that you can numerically express it as a function of t, you can use maximum-likelihood-based approaches to estimate μ. You can use the nlmixr2 package to express your model as an ODE and optimize the value of μ given your time-varying data. This package is typically applied for non-linear mixed effects in a population setting in the field of pharmacometrics, but you should be able to adapt it to your case (eg, drop the variability component).

This might be a case of "when you have hammer, everything looks like a nail", but this should do the trick.

3

Thank you for your response! I appreciate the suggestion to use nlmixr2. I hadn’t considered adapting it for this problem, so I’ll definitely look into it. However, I’m still a bit unsure about how to incorporate the non-parametric modeling of φ(t) as a Gaussian process within that framework.

To clarify my goal: I’m trying to estimate μ as a fixed parameter while treating φ(t) as a time-varying function modeled non-parametrically (ideally as a Gaussian process). The challenge is that φ(t) is not directly observed, instead it’s part of the ODE structure, so I need a method that can handle this kind of latent function estimation alongside parameter inference.

Thanks again for your help!

4

You might be able to do this with Stan, which fits Bayesian models, including ODE models, Gaussian processes, and latent parameters. I haven't worked with ODE models in Stan but you might ask your question on the Stan discourse site, where folks are super helpful.

It's even possible that you could fit your model using the R brms package, which provides a high-level formula interface to Stan, similar to the model formulas used for fitting frequentist models (e.g., the glm and lme4::lmer functions), so you wouldn't need to write Stan code (here's an example of fitting an ODE model using brms).

2 Likes
5

Thank you for the suggestion! I think brms might indeed work for my case, and I’ve been exploring it further. I’ll definitely give it a try and see how it performs with my specific problem. Thanks again for pointing me in this direction, it’s been really helpful!

6

Have you looked into using the deSolve package for solving the ODE and then optimizing the parameters with something like optim or a Bayesian optimization package like rstan or greta? You could approximate φ(t) with splines or a similar basis function approach within the ODE.

7

Thanks for the suggestion! I’d prefer to keep φ(t) non-parametric rather than approximating it with splines or a basis function approach. Ideally, I’d like a method that allows for flexible inference of φ(t) while also estimating μ.

8

The suggestion by @joels seems like an obvious one to me, because a Bayesian approach would help in a small sample size setting: Using brms (or Stan) using the integration routines in Stan in the background. There's even an [old example around]( Non-linear model with integrate_1d in Bmrs - Interfaces / brms - The Stan Forums (I hope the issue the poster mentioned is fixed by now) that shows how to do integration in brms. brms should also let you use a Gaussian process and the times I've used them it's been pretty straighforward.

9

Thank you for your suggestion! I agree that brms seems like a good approach for this problem, and I’ve tried implementing it with the following code:

library(brms)

data <- data.frame(
  t = c(0, 1, 2, 3, 4, 5, 6, 7),
  n = c(0.00, 0.39, 0.57, 0.50, 0.12, 0.06, 0.02, 0.00),
  dndt = c(0.39, 0.28, 0.05, -0.23, -0.22, -0.05, -0.03, -0.02)
)

formula <- bf(
  dndt ~ phi - nu * n,
  phi ~ gp(t),
  nu ~ 1,
  nl = TRUE
)

priors <- set_prior("gamma(2, 0.1)", nlpar = "nu", lb = 0)

fit <- brm(
  formula = formula,
  data = data,
  prior = priors,
  family = gaussian(),
  control = list(adapt_delta = 0.95),
  chains = 4,
  iter = 2000,
  cores = 4,
  algorithm = "sampling",
  silent = 0
)

summary(fit)
plot(fit)

However, I’m running into an error, and I’m not sure how to fix it. I’ve also posted this issue on the Stan forum here, but I haven’t had much success yet.

Could you or anyone else in the community help me troubleshoot this? Specifically:

  1. Is there something wrong with how I’ve specified the formula?
  2. Are there any common pitfalls when using Gaussian processes in brms that I might be missing?

Any guidance would be greatly appreciated!

Thanks in advance!