aim: Adaptive Integral Matching (AIM)
In episode: Estimation with Penalisation in Systems of Ordinary Differential Equations

Description Usage Arguments Details Value See Also Examples

Gives approximate parameter estimates using integral matching and optionally adapts weights and scales to these. Feed this to rodeo for initialising exact estimation.

1 2	aim(o, op, x = NULL, adapts = c("scales", "penalty_factor"), xout = FALSE, params = NULL, ...)

`o`	An object of class `ode` (created via `mak`, `plk`, etc.).
`op`	An object of class `opt` (created via `opt`). Compatibility check with `o` is conducted.
`x`	Matrix of dimension mx(d+1) containing custom time and smoothed values of the processes used for the integral matching, see details. If `NULL` (default) a linear-interpolation smoother is employed.
`adapts`	Character vector holding names of what quantities to adapt during algorithm. Possible quantities are: `"scales"`, `"weights"` and `"penalty_factors"`, see details. Default is `"scales"` and `"penalty_factor"`.
`xout`	Logical indicating if a matrix containing the process values at the time points in `y` in `op`, linearly interpolated from `x` should be returned. Default is `FALSE`.
`params`	List of vectors of initial parameter values. Default (`NULL`) translates to the zero-vector. If the ODE system is linear in the parameter vector and the boundary constraints on are not unusual, then `params` is ignored.
`...`	Additional arguments passed to `aim`.

Loss function

Integral matching requires a smoothed process in order to approximate the parameter estimates. More precisely, the loss function

RSS / (2 * (n - s)) + lambda*penalty

is minimised, where RSS is the sum of the squared 2-norms of

x_i - x_{i-1} - \int_{t_{i-1}}^{t_{i}}{f(x(s), context_l * param) ds}

Here f is the vector field of the ODE-system, x is the smoothed process and param is (internally) scaled with scales in reg.

Custom smoother

The supplied x is a way of customising how x in the loss function is made. Firstly, x must have similiar layout as y in op, i.e., first column is time and the remaining columns contain smoothed values of the process at those time points, see opt-documentation for details.

The number of decreases in time in x must match the number of decreases in time in y in op. The process x in the loss function is simply a linear interpolation of x, hence finer discretisations give more refined integral matching. Each context is handled seperately. Moreover, the compatibility checks in rodeo are also conducted here.

Adaptations

The adapts are ways to adapt quantities in the opt-object and reg-objects in o to the data:

"scales": A standardisation of the columns in the linearisation takes place and carry over to scales in reg-objects in o (generally recommended and default).
"weights": The observational weights (weights in op) are adjusted coordinate-for-coordinate (column-wise) by reciprocal average residual sum of squares across penalty parameters.
"penalty_factor": The penalty factors in reg are adjusted by the reciprocal average magnitude of the estimates (parameters whose average magnitude is 0 join exclude). For extremely large systems, this option can heavily reduce further computations if the returned aim-object is passed to rodeo.

If either "penalty_factor" or "weights" are in adapts a refitting takes place.

Finally, note that lambda and nlambda in returned opt-object may have changed.

An object with S3 class "aim":

`o`	Original `ode`-object with adapted quantities.
`op`	Original `opt`-object with adapted quantities and default values for `lambda_min_ratio`, `lambda` and (if needed) `a` inserted, if these were originally `NULL`.
`params`	A list of matrices (of type dgCMatrix) with the parameter estimates (scaled by `scales` in `reg`), one column for each lambda value. One element in the list for each parameter (other than initial state parameter).
`x0s`	A matrix with the initial state estimates.
`rss`	A vector of residual sum of squares.
`x`	Original `x`, or `y` in `op` if `x` was `NULL`.
`xout`	(If `xout` = `TRUE`) A matrix containing the values of the process at the time points in `y` in `op`, interpolated from `x`. Use it to check that `x` works as intended.

rodeo, rodeo.ode, rodeo.aim, imd

set.seed(123)
# Michaelis-Menten system with two 0-rate reactions
A <- matrix(c(1, 1, 0, 0,
  0, 0, 1, 0,
  0, 0, 1, 0,
  0, 0, 1, 1,
  0, 0, 1, 1), ncol = 4, byrow = TRUE)
B <- matrix(c(0, 0, 1, 0,
  1, 1, 0, 0,
  1, 0, 0, 1,
  0, 1, 0, 0,
  1, 0, 0, 1), ncol = 4, byrow = TRUE)
k <- c(0.1, 1.25, 0.5, 0, 0); x0 <- c(E = 5, S = 5, ES = 1.5, P = 1.5)
Time <- seq(0, 10, by = 1)

# Simulate data, in second context the catalytic rate has been inhibited
contexts <- cbind(1, c(1, 1, 0, 1, 1))
m <- mak(A, B, r = reg(contexts = contexts))
y <- numsolve(m, c(Time, Time), cbind(x0, x0 + c(2, -2, 0, 0)), contexts * k)
y[, -1] <- y[, -1] + matrix(rnorm(prod(dim(y[, -1])), sd = .25), nrow = nrow(y))

# Create optimisation object via data
o <- opt(y, nlambda = 10)

# Fit data using Adaptive Integral Matching on mak-object
a <- aim(m, o)
a$params$rate

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


# 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
a <- aim(p, opt(y, nlambda = 10))
a$params$theta

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


# Example: use custom lowess smoother on data
# Smooth each coordinate of data to get curve
# on extended time grid
ext_time <- seq(0, 1, by = 0.001)
x_smooth <- apply(y[, -1], 2, function(z) {
  # Create loess object
  data <- data.frame(Time = y[, -1], obs = z)
  lo <- loess(obs ~ Time, data)

  # Get smoothed curve on extended time vector
  predict(lo, newdata = data.frame(Time = ext_time))
})

# Bind the extended time
x_smooth <- cbind(ext_time, x_smooth)

# Run aim on the custom smoothed curve
a_custom <- aim(p, opt(y), x = x_smooth)