pCODE (version 0.9.4)

pcode: Parameter Cascade Method for Ordinary Differential Equation Models

Description

Obtain estimates of both structural and nuisance parameters of an ODE model by parameter cascade method.

Usage

pcode(data, time, ode.model, par.names, state.names,
             likelihood.fun, par.initial, basis.list,lambda,controls)

Value

structural.par

The structural parameters of the ODE model.

nuisance.par

The nuisance parameters or the basis coefficients for interpolating observations.

Arguments

data

A data frame or a matrix contain observations from each dimension of the ODE model.

time

A vector contain observation times or a matrix if time points are different between dimensions.

ode.model

An R function that computes the time derivative of the ODE model given observations of states variable and structural parameters.

par.names

The names of structural parameters defined in the 'ode.model'.

state.names

The names of state variables defined in the 'ode.model'.

likelihood.fun

A likelihood function passed to PCODE in case of that the error terms do not have a Normal distribution.

par.initial

Initial value of structural parameters to be optimized.

basis.list

A list of basis objects for smoothing each dimension's observations. Can be the same or different across dimensions.

lambda

Penalty parameter for controling the fidelity of interpolation.

controls

A list of control parameters. See Details.

Details

The controls argument is a list providing addition inputs for the nonlinear least square optimizer or general optimizer optim:

nquadpts

Determine the number of quadrature points for approximating an integral. Default is 101.

smooth.lambda

Determine the smoothness penalty for obtaining initial value of nuisance parameters.

tau

Initial value of Marquardt parameter. Small values indicate good initial values for structural parameters.

tolx

Tolerance for parameters of objective functions. Default is set at 1e-6.

tolg

Tolerance for the gradient of parameters of objective functions. Default is set at 1e-6.

maxeval

The maximum number of evaluation of the outter optimizer. Default is set at 20.

Examples

library(fda)
library(deSolve)
library(MASS)
library(pracma)
#Simple ode model example
#define model parameters
model.par   <- c(theta = c(0.1))
#define state initial value
state       <- c(X     = 0.1)
#Define model for function 'ode' to numerically solve the system
ode.model <- function(t, state,parameters){
 with(as.list(c(state,parameters)),
      {
        dX <- theta*X*(1-X/10)
        return(list(dX))
      })
}
#Observation time points
times <- seq(0,100,length.out=101)
#Solve the ode model
desolve.mod <- ode(y=state,times=times,func=ode.model,parms = model.par)
#Prepare for doing parameter cascading method
#Generate basis object for interpolation and as argument of pcode
#21 konts equally spaced within [0,100]
knots <- seq(0,100,length.out=21)
#order of basis functions
norder <- 4
#number of basis funtions
nbasis <- length(knots) + norder - 2
#creating Bspline basis
basis  <- create.bspline.basis(c(0,100),nbasis,norder,breaks = knots)
#Add random noise to ode solution for simulating data
nobs  <- length(times)
scale <- 0.1
noise <- scale*rnorm(n = nobs, mean = 0 , sd = 1)
observation <- desolve.mod[,2] + noise
#parameter estimation
pcode(data = observation, time = times, ode.model = ode.model,
                     par.initial = 0.1, par.names = 'theta',state.names = 'X',
                     basis.list = basis, lambda = 1e2)