Description Usage Arguments Details Value References Examples
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.
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(),
...
)
|
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 |
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 |
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 |
Named list. The observations. When |
obs_sets |
Number of observations sets. When |
nlin_pars |
Names of parameters or initial conditions
that will be estimated in stage 1 using nonlinear least squares optimization.
The parameter names in |
likelihood_pars |
Names of likelihood parameters not appearing in the ODE system,
which are needed for the the user-defined function |
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 |
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 |
... |
Additional arguments passed to |
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:
equations
the ODE equations
pars
the parameter values
x0
the initial conditions
time
the timing of the observations (vector or list)
obs
the observations
...
additional parameters passed from the call to simode
The function should return a list with two items:
time
the vector or list of time points of the observations
obs
the list of observations with the newly generated observations
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:
pars
the parameter values
time
the timing of the observations (vector or list)
obs
the observations
model_out
the model output returned from a call to solve_ode
.
If time
is a list with possibly different times for each variable
then model_out
will contain a union of all these times.
...
additional parameters passed from the call to simode
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:
pars
the parameter values
...
additional parameters passed from the call to simode
The function should return the rescaled parameter values.
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
.
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.
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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.