Skip to contents

Derivative Function Estimation (Step 2 of Kernel ODE)

Source: R/kernelODE_step2.R
kernelODE_step2.Rd

This function implements an iterative optimization algorithm as described in the Kernel ODE paper.

Usage

kernelODE_step2(
  Y,
  obs_time,
  yy_smth,
  tt,
  kernel,
  kernel_params,
  interaction_term = FALSE,
  theta_initial = NULL,
  adj_matrix = NULL,
  nzero_thres = NULL,
  eval_loss = FALSE,
  tol = 0.001,
  max_iter = 10,
  verbose = 0
)

Arguments

Y

A numeric matrix of dimension (n, p), where each column corresponds to the observed trajectory of a variable. Rows align with obs_time.

obs_time

A numeric vector of length n representing observation time points.

yy_smth

A numeric matrix of dimension (length(tt), p), where each column contains the smoothed trajectory of a variable evaluated on tt. This is typically the output of kernelODE_step1().

tt

A numeric vector representing a finer time grid used for evaluating the smoothed trajectories and their derivatives.

kernel

Kernel function to use.

kernel_params

A list of length p, where each element is a named list of parameters for a specific variable (e.g., list(bandwidth = 1) for Gaussian kernel). If the list has length 1, the same parameter set is used for all variables. This is typically the output of auto_select_kernel_params().

interaction_term

A logical value specifying whether to include interaction effects in the model.

theta_initial

A numeric matrix containing the initial \(\theta_j\) values for each variable at the start of optimization:

  • If interaction_term = FALSE, theta_initial must have dimensions (p, p).

  • If interaction_term = TRUE, theta_initial must have dimensions (p^2, p). If NULL, defaults to all ones.

adj_matrix

An adjacency matrix (p × p) representing a fixed regulatory network, where entry (k, j) = 1 indicates variable k regulates variable j. If supplied, no Lasso selection is performed, and non-penalized regression is used instead. If NULL (default), the set of regulators (edges) is selected via Lasso-type penalization.

nzero_thres

A number in [0, 1] specifying the maximum proportion of nonzero regulators (edges) allowed for each variable (i.e., at most p * nzero_thres regulators). Only used when adj_matrix is NULL. This provides a faster alternative to the computationally expensive pruning process.

eval_loss

A logical value specifying whether to evaluate the loss function at each iteration.

tol

Convergence tolerance for the relative improvement in the Frobenius norm of the \(\theta\) matrix.

max_iter

Maximum number of iterations for the optimization.

verbose

Integer; if greater than 0, prints progress messages during optimization.

Value

A list with components:

res_theta

A numeric matrix giving the estimated \(\theta_j\) values, with each column corresponding to one variable \(j\). It dimension matches that of theta_initial.

res_best_kappa

A numeric vector of length p with selected hyperparameter \(\kappa_j\) values for the \(\theta_j\) estimation step.

res_bj

A numeric vector of length p with estimated \(b_j\) values.

res_cj

A numeric matrix (n × p) with estimated \(c_j\) values in columns.

res_best_eta

A numeric vector of length p with selected hyperparameter \(\eta_j\) values from the \(F_j\) estimation step (i.e., estimating \(b_j\) and \(c_j\)), chosen by generalized cross-validation.

res_loss_path

Optional list of loss values over iterations (present only if eval_loss = TRUE).

network_est

Estimated regulatory network, in the form of adjacency matrix.

num_iter

Number of iterations performed.

config

List of input arguments for reproducibility.

References

Dai, X., & Li, L. (2022). Kernel ordinary differential equations. Journal of the American Statistical Association, 117(540), 1711-1725.

Examples

set.seed(1)
obs_time <- seq(0, 1, length.out = 10)
Y <- cbind(sin(2 * pi * obs_time), cos(4 * pi * obs_time)) + 0.1 * matrix(rnorm(20), 10, 2)  # each col is a variable
tt <- seq(0, 1, length.out = 100)
res_step1 <- kernelODE_step1(Y = Y, obs_time = obs_time, tt = tt)

kernel <- "gaussian"
kernel_params <- auto_select_kernel_params(kernel = kernel, Y = Y)
res_step2 <- kernelODE_step2(Y = Y, obs_time = obs_time, yy_smth = res_step1$yy_smth, tt = tt, kernel = kernel, kernel_params = kernel_params, interaction_term = FALSE)
network_est <- res_step2$network_est
network_est
#>      [,1] [,2]
#> [1,]    1    1
#> [2,]    1    1