episode (version 1.0.0)

aim: Adaptive Integral Matching (AIM)

Description

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

Usage

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

Arguments

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.

Value

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.

Details

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.

See Also

rodeo, rodeo.ode, rodeo.aim, imd

Examples

# NOT RUN {
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)


# }