episode (version 1.0.0)

rodeo.aim: Regularised Ordinary Differential Equation Optimisation (RODEO) initialised via Adaptive Integral Matching

Description

Fit the parameters for an ODE model with data sampled across different contexts.

Usage

# S3 method for aim
rodeo(x, adjusts = c("lambda"), trace = FALSE, ...)

Arguments

x

aim-object created via aim-function.

adjusts

Character vector holding names of what quantities to adjust during algorithm. Possible quantities are: "lambda", "scales" and "weights".

trace

Logical indicating if status messages should be printed during rodeo.

...

Additional arguments passed to rodeo.

Value

An object with S3 class "rodeo":

o

Original ode-object.

op

Original opt-object with default values for lambda_min_ratio, lambda and (if needed) a inserted, if these were originally NULL.

params

Parameter estimates, stored as list of sparse column format matrices, "dgCMatrix" (or a list of those if multiple initialisations). Rows represent coordinates and columns represent the lambda value.

x0s

Initial state estimates stored in a matrix (or array). Rows represent coordinates, columns represent the lambda value and (if multiple initialisations) slices represent initialisations.

dfs

A matrix (or array, if multiple initialisations) of degrees of freedom. Row represents a parameter (the first is always the initial state parameter), columns represent lambda, slices represent initialisation, if multiple are provided.

codes

A matrix (or array) of convergence codes organised as dfs.

0:

The convergence criteria is met (see details in opt). Current estimate is probably a local minimum.

1:

Backtracking in the last iteration yields no numerical improvement, but no unsual behavior observed. Current estimate is probably a local minimum. However, if exact_gradient = FALSE in the reg-object in the ode-object, changing this may improve the code. Alternatively one can adjust backtracking via backtrack_max and tau_min in reg objects in ode object.

2:

The optimisation procedure exceeded maximal number of steps (step_max in reg objects).

3:

The last gradient was unsually large. Either the tolerances in reg objects are off or the ODE systems is very sensitive and runs over long time spans. In the latter case, initialisation(s) may have inappropriate zeros (change initialisation and/or make sure they start at smaller lambda value).

4:

The numeric ODE solver exceeded maximal number of steps. Check if supplied initial states were out of bounds, if not increase step_max (or tol) in reg-objects in ode-object.

steps

A matrix (or array) holding number of steps used in optimisation procedure. Organised as dfs.

losses

A vector (or matrix) of unpenalised losses at optimum for each lambda value (stored row-wise if multiple are provided).

penalties

A matrix (or array) of penalties for each parameter, organised as dfs.

jerr

A matrix (or array) of summary codes (for internal debugging), organised as dfs.

Details

The adapted quantities (scales, weights, penalty_factor) of x (returned by aim) are fed to the exact estimator rodeo. This estimator then traverses the lambda sequence in reverse order initialised in the last estimates from aim.

If desired, the quantities lambda, scales and weights are adjusted as in aim.

See Also

rodeo, aim, rodeo.ode

Examples

# NOT RUN {
set.seed(123)
# Example: Power Law Kinetics
A <- matrix(c(1, 0, 1,
  1, 1, 0), byrow = TRUE, nrow = 2)
p <- plk(A)
x0 <- c(10, 4, 1)
theta <- matrix(c(0, -0.25,
  0.75, 0,
  0, -0.1), byrow = TRUE, nrow = 3)
Time <- seq(0, 1, by = .025)

# Simulate data
y <- numsolve(p, Time, x0, theta)
y[, -1] <- y[, -1] + matrix(rnorm(prod(dim(y[, -1])), sd = .25), nrow = nrow(y))

# Estimation via aim
a <- aim(p, opt(y, nlambda = 10))
a$params$theta

# Supply to rodeo
rod <- rodeo(a)
rod$params$theta

# Compare with true parameter on column vector form
matrix(theta, ncol = 1)


# Example: include data from an intervened system
# where the first complex in A is inhibited
contexts <- cbind(1, c(0, 0, 0, 1, 1, 1))
y2 <- numsolve(p, Time, x0 + 1, theta * contexts[, 2])
y2[, -1] <- y2[, -1] + matrix(rnorm(prod(dim(y2[, -1])), sd = .25), nrow = nrow(y2))

# Estimation via aim
a <- aim(plk(A, r = reg(contexts = contexts)), opt(rbind(y, y2), nlambda = 10))
a$params$theta

# Supply to rodeo
rod <- rodeo(a)
rod$params$theta


# }