simode: Statistical inference of ordinary differential equations...

Description Usage Arguments Details Value References Examples

View source: R/simode.R

Description

Estimating the parameters of an ODE system in two stages: 1) Estimate the parameters using separable integral-matching, 2) Estimate the parameters using nonlinear least squares starting from the values obtained in stage 1.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
simode(
  equations,
  pars,
  time,
  obs,
  obs_sets = 1,
  nlin_pars = NULL,
  likelihood_pars = NULL,
  fixed = NULL,
  start = NULL,
  lower = NULL,
  upper = NULL,
  im_method = c("separable", "non-separable"),
  decouple_equations = FALSE,
  gen_obs = NULL,
  calc_nll = NULL,
  scale_pars = NULL,
  simode_ctrl = simode.control(),
  ...
)

Arguments

equations

Named vector. The equations describing the ODE system. Each element of the vector should contain a character representation of the right-hand side of an equation, and should be named according to the left-hand side of the equation (i.e., the variable name). An equation can contain parameters appearing in pars, variables appearing in the equations names, observed non-modeled variables appearing in obs, and/or any function of 't', which is a reserved symbol for the time domain.

pars

The names of the parameters and initial conditions to be estimated. An initial condition name for a certain variable is the name given to the relevant equation in equations (e.g., if an equation is named 'x' than its initial condition should be named 'x' as well). Note: The symbol 't' is reserved for the time domain and cannot be used as a parameter name.

time

Time points of the observations. Either a vector, if the same time points were used for observing all variables, or a list of vectors the length of obs, of which each element is the length of the relevant element in obs.

obs

Named list. The observations. When obs_sets=1, obs should contain a list of vectors with observations of either a variable described by one of the equations (named according to the relevant equation name) or a non-modeled variable appearing in one of the equations. Each observations vector should be the length of the relevant time vector. When obs_sets>1, obs should contain a list, where each list member is a list that fits the description in the case of obs_sets=1.

obs_sets

Number of observations sets. When obs_sets>1, the function will fit the set of observations according to the value of obs_sets_fit in simode_ctrl.

nlin_pars

Names of parameters or initial conditions that will be estimated in stage 1 using nonlinear least squares optimization. The parameter names in nlin_pars must appear in pars.

likelihood_pars

Names of likelihood parameters not appearing in the ODE system, which are needed for the the user-defined function calc_nll. The parameter names in likelihood_pars must appear in pars.

fixed

Named vector. Fixed values for one or more of the ODE system parameters or initial conditions. Parameters in this list will not be estimated.

start

Named vector. Starting values for optimization of parameters/initial conditions. Must contain starting values for all the parameters in nlin_pars and likelihood_pars. If im_method="non-separable", can optionally contain starting values for any other parameter/initial condition.

lower

Named vector. Lower bounds for any parameter/initial condition.

upper

Named vector. Upper bounds for any parameter/initial condition.

im_method

The method to use for integral-matching. Default "separable" means that linear parameters are estimated directly while "non-separable" means that linear parameters are estimated using nonlinear least squares optimization. If none of the parameters are linear then the default can be used.

decouple_equations

Whether to fit each equation separately in the integral-matching stage.

gen_obs

An optional user-defined function for completing missing observations (see Details).

calc_nll

An optional user-defined function for calculating negative log-likelihood for the model (see Details).

scale_pars

An optional user-defined function for scaling transformations of parameters estimated using non-linear optimization (see Details).

simode_ctrl

Various control parameters. See simode.control.

...

Additional arguments passed to optim, gen_obs and calc_nll

Details

gen_obs can be used in cases of a partially observed system, for which observations of the missing variables can be generated given values for the system parameters. The function will be called during the optimization using integral-matching.
It must be defined as gen_obs <- function(equations, pars, x0, time, obs, ...), where:

The function should return a list with two items:

calc_nll allows the user to pass his own likelihood function to be used in the optimization in the second stage (if not defined, the default nonlinear least squares optimization will be used). The likelihood function will also be used in a following call to profile, for the calculation of likelihood profiles. It must be defined as calc_nll <- function(pars, time, obs, model_out, ...), where:

The function should return the negative log-likelihood.

scale_pars allows the user to pass a function for rescaling transformations of some/all of the parameters estimated using non-linear optimization. The function will be called at each step of the optimization with the current values of the parameters, as the optimization algorithm sees them, and will return rescaled parameters values to be used in estimation of the direct parameters (stage 1) and in solving the ODE equations (stage 2). It must be defined as scale_pars <- function(pars, ...), where:

The function should return the rescaled parameter values.

Value

If obs_sets=1, the function returns a simode object containing the parameter estimates and solutions after integral-matching (stage 1) and after nonlinear least squares optimization (stage 2). If obs_sets>1 and obs_sets_fit!="together" in simode_ctrl, the function returns a list.simode object which is a list of simode objects the length of obs_sets.

References

Dattner & Klaassen (2015). Optimal Rate of Direct Estimators in Systems of Ordinary Differential Equations Linear in Functions of the Parameters, Electronic Journal of Statistics, Vol. 9, No. 2, 1939-1973.

Dattner, Miller, Petrenko, Kadouriz, Jurkevitch & Huppert (2017). Modelling and Parameter Inference of Predator-prey Dynamics in Heterogeneous Environments Using The Direct Integral Approach, Journal of The Royal Society Interface 14.126: 20160525.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
## =================================================
## Predator-Prey Lotka-Volterra model
## =================================================

## generate model equations and parameters (X=Prey,Y=Predator)
pars <- c('alpha','beta','gamma','delta')
vars <- c('X','Y')
eq_X <- 'alpha*X-beta*X*Y'
eq_Y <- 'delta*X*Y-gamma*Y'
equations <- c(eq_X,eq_Y)
names(equations) <- vars
x0 <- c(0.9,0.9)
names(x0) <- vars
theta <- c(2/3,4/3,1,1)
names(theta) <- pars

## generate observations
n <- 50
time <- seq(0,25,length.out=n)
model_out <- solve_ode(equations,theta,x0,time)
x_det <- model_out[,vars]
set.seed(1000)
sigma <- 0.05
obs <- list()
for(i in 1:length(vars)) {
  obs[[i]] <- pmax(0, rnorm(n,x_det[,i],sigma))
}
names(obs) <- vars

## estimate model parameters with known initial conditions
simode_fit1 <- simode(equations=equations, pars=pars, fixed=x0, time=time, obs=obs)
plot(simode_fit1, type='fit', time=seq(0,25,length.out=100), pars_true=theta, mfrow=c(2,1))
plot(simode_fit1, type='est', pars_true=theta)


## estimate model parameters and initial conditions
simode_fit2 <- simode(equations=equations, pars=c(pars,vars), time=time, obs=obs)
plot(simode_fit2, type='fit', time=seq(0,25,length.out=100), pars_true=c(theta,x0), mfrow=c(2,1))
plot(simode_fit2, type='est', pars_true=c(theta,x0))

profiles_fit2 <- profile(simode_fit2,step_size=0.01,max_steps=50)
plot(profiles_fit2,mfrow=c(2,3))
ci_fit2 <- confint(profiles_fit2)
ci_fit2
plot(ci_fit2,pars_true=c(theta,x0),legend=TRUE)

simode documentation built on July 1, 2020, 10:30 p.m.