Package 'multivar'

Description

Simulate, estimate, and forecast vector autoregressive (VAR) models for multiple-subject data using structured penalization. Decomposes dynamics into shared (common) and subject-specific (unique) components via adaptive LASSO with FISTA optimization. Supports cross-validation and extended BIC model selection and subgroup detection, and time-varying parameters.

multivar is an R package for simulating, estimating and forecasting stationary Vector Autoregressive (VAR) models for multiple subject data using the penalized multi-VAR framework.

Author(s)

Maintainer: Zachary Fisher [email protected]

Authors:

• Christopher Crawford

Other contributors:

• Younghoon Kim [contributor]

• Vladas Pipiras [contributor]

Description

Creates a specification object that defines the column and row structure of the design matrix A. This is the single source of truth for all downstream functions that need to know column indices.

Usage

build_matrix_spec(
  k,
  d,
  n,
  tvp = FALSE,
  common_effects = TRUE,
  subgroup = FALSE,
  common_tvp_effects = FALSE,
  breaks = NULL,
  subgroup_membership = NULL
)

Arguments

Details

The design matrix A has the following column structure depending on model type:

k=1, non-TVP:

[d columns]

k=1, TVP, common_effects=TRUE:

[d common | d*P unique_tvp]

k=1, TVP, common_effects=FALSE:

[d*P unique_tvp]

k>1, non-TVP, no subgroup:

[d common | d unique_1 | ... | d unique_k]

k>1, non-TVP, with subgroup:

[d common | d subgrp_1 | ... | d subgrp_S | d unique_1 | ... | d unique_k]

k>1, TVP, common_effects=TRUE, common_tvp_effects=TRUE:

[d common | d unique_1 | ... | d unique_k | d*P common_tvp | d*P unique_tvp_1 | ... | d*P unique_tvp_k]

Value

A list with class "matrix_spec" containing:

• params: Model parameters (k, d, n, flags)

• cols: Column index ranges for each effect type

• rows: Row index ranges for subjects and periods

Description

This matches the logic in est_base_weight_mat(): we expand ratios_unique, and (optionally) ratios_subgroup and ratios_unique_tvp into a list of S scenarios, where S = dim(W)[3].

Usage

build_scenario_table(object)

Arguments

Description

Canonical VAR Fitting Function for multivar

Usage

canonical.multivar(object)

Arguments

Details

A function to fit a canonical VAR model to each individual dataset.

Value

A list of results.

See Also

constructModel,

Examples

# example 1 (run)
sim1  <- multivar_sim(
  k = 2,  # individuals
  d = 5,  # number of variables
  n = 20, # number of timepoints
  prop_fill_com = 0.1, # proportion of paths common
  prop_fill_ind = 0.05, # proportion of paths unique
  lb = 0.1,  # lower bound on coefficient magnitude
  ub = 0.5,  # upper bound on coefficient magnitude
  sigma = diag(5) # noise
)

model1 <- constructModel(data = sim1$data, weightest = "ols")
fit1 <- canonical.multivar(model1)

Description

Core computation of the Extended Bayesian Information Criterion (EBIC) for each scenario in a beta array. Used internally by cv.multivar when selection = "ebic" and by select_by_ebic for post-hoc reselection.

Usage

compute_ebic(beta_array, Z, Y, d, gamma = 0.5, spec = NULL)

Arguments

Details

With spec = NULL (legacy behavior), EBIC counts every nonzero entry in B as one degree of freedom:

$\text{EBIC}_i = \sum_{j=1}^{d} n \log(\text{RSS}_{ij} / n) + k_i \log(n) + 2 \gamma \, k_i \log(p)$

With a spec object, degrees of freedom are counted per position: for each d x d coefficient position (i,j), the position contributes 1 df if any block (common, any unique, any TVP) has a nonzero entry at that position. The eBIC extension uses p = d (predictors per equation) instead of the full design matrix width, and employs the Chen & Chen (2008) lchoose(p, df) formulation per equation.

Value

Numeric vector of EBIC values, one per scenario.

Description

Construct an object of class multivar

Usage

constructModel(
  data = NULL,
  lag = 1,
  horizon = 0,
  t1 = NULL,
  t2 = NULL,
  lambda1 = NULL,
  nlambda1 = 30,
  n_ratios_subgroup = 30,
  depth = NULL,
  tol = 1e-04,
  window = 1,
  standardize = TRUE,
  weightest = "maity",
  canonical = FALSE,
  threshold = FALSE,
  lassotype = "adaptive",
  intercept = FALSE,
  W = NULL,
  ratios_unique = NULL,
  ratios_subgroup = NULL,
  ratios_unique_tvp = NULL,
  cv = "blocked",
  nfolds = 10,
  lamadapt = FALSE,
  subgroup_membership = NULL,
  subgroup = FALSE,
  B = NULL,
  pendiag = TRUE,
  tvp = FALSE,
  inittvpcoefs = list(),
  breaks = list(),
  lambda_choice = "lambda.min",
  common_effects = TRUE,
  common_tvp_effects = NULL,
  save_beta = TRUE,
  ncores = 1,
  max_grid_size = NULL,
  eps = 0.001,
  warmstart = TRUE,
  stopping_crit = "absolute",
  selection = "cv",
  ebic_gamma = 0.5,
  weight_type = "standard",
  maity_opts = list()
)

Arguments

Examples

sim  <- multivar_sim(
  k = 2,  # individuals
  d = 3,  # number of variables
  n = 20, # number of timepoints
  prop_fill_com = 0.1, # proportion of paths common
  prop_fill_ind = 0.1, # proportion of paths unique
  lb = 0.1,  # lower bound on coefficient magnitude
  ub = 0.9,  # upper bound on coefficient magnitude
  sigma = diag(3) # noise
)

plot_sim(sim, plot_type = "common")

model <- constructModel(data = sim$data, weightest = "ols")

Description

Performs k-fold blocked cross-validation for time series data.

Usage

cv_blocked(
  B,
  Z,
  Y,
  W,
  Ak,
  k,
  d,
  lambda1,
  t1,
  t2,
  eps,
  intercept = FALSE,
  cv,
  nfolds,
  tvp = FALSE,
  breaks = NULL,
  spec = NULL,
  ncores = 1,
  warmstart = FALSE,
  stopping_crit = 0L
)

Arguments

Value

List with beta coefficients and MSFE matrix

Description

Routes to the appropriate CV method (blocked or rolling).

Usage

cv_multivar(
  B,
  Z,
  Y,
  W,
  Ak,
  bk,
  k,
  d,
  lambda1,
  t1,
  t2,
  eps,
  intercept = FALSE,
  cv,
  nfolds,
  tvp = FALSE,
  breaks = NULL,
  spec = NULL,
  ncores = 1,
  warmstart = FALSE,
  stopping_crit = 0L
)

Arguments

Value

List with beta coefficients and MSFE matrix

Description

Cross Validation for multivar

Usage

cv.multivar(object)

Arguments

Details

The main function of the multivar package. Performs cross validation to select penalty parameters over a training sample and evaluates them over a test set.

Value

An object of class multivar.results.

Examples

# example 1 (run)
sim1  <- multivar_sim(
  k = 2,  # individuals
  d = 5,  # number of variables
  n = 20, # number of timepoints
  prop_fill_com = 0.1, # proportion of paths common
  prop_fill_ind = 0.05, # proportion of paths unique
  lb = 0.1,  # lower bound on coefficient magnitude
  ub = 0.5,  # upper bound on coefficient magnitude
  sigma = diag(5) # noise
)

model1 <- constructModel(data = sim1$data)
fit1 <- multivar::cv.multivar(model1)

Description

This dataset contains multivariate time series data for k = 9 individuals with $d = 10$ variables collected at $t = 100$ equidistant time points. The data was generated such that each individual's VAR(1) transition matrix has 20 percent nonzero entries. This means, for example, each individual has 20 nonzero directed relationships in their data generating model. The position of non-zero elements in each individual's transition matrix was selected randomly given the following constraints: 2/3 of each individual's paths are shared by all individuals, and 1/3 are unique to each individual. For each individual, coefficient values between U(0,1, 0.9) were randomly drawn until stability conditions for the VAR model were satisfied.

Usage

dat_multivar_sim

Format

A list containing

mat_com

a common effects transition matrix

mat_ind_unique

a list of unique (individual-specific) effect matrices

mat_ind_final

a list of total (common + individual-specific) effect matrices

data

a list of multivariate time series for all subjects

...

Description

'estimate_initial_coefs' returns consistent initial estimates to be used for calculating adaptive weights in adaptive LASSO regularization.

Usage

estimate_initial_coefs(
  Ak,
  bk,
  d,
  k,
  lassotype,
  weightest,
  subgroup_membership,
  subgroup,
  nlambda1,
  tvp,
  breaks,
  intercept,
  nfolds,
  lambda_choice = "lambda.min",
  common_effects = TRUE,
  common_tvp_effects = TRUE,
  maity_opts = list()
)

Arguments

Value

A list containing:

• common_effects: d x d matrix of common effects

• subgroup_effects: list of d x d matrices per subgroup (if subgroup=TRUE)

• unique_effects: list of d x d matrices per subject

• tvp_effects: list of lists of d x d matrices per subject per period (if tvp=TRUE)

• common_tvp_effects: list of d x d matrices per period (if tvp=TRUE and k>1)

• total_effects: list of d x d matrices per subject

Description

Compare estimated transition matrices to simulation truth, but only for components that actually exist in both 'sim_obj' and 'fit_obj'. For example, if 'fit_obj$mats' only contains 'total', only "total" rows are returned.

Usage

eval_multivar_performance(
  sim_obj,
  fit_obj,
  eps = 1e-08,
  reduced.output = TRUE,
  averages.only = TRUE,
  label = NULL,
  intercept = FALSE
)

Arguments

Description

Converts a c(start, end) range to a full sequence.

Usage

expand_range(range)

Arguments

Value

Integer vector or NULL

Description

Creates a tidy dataframe with all tested hyperparameter combinations and their cross-validation error (MSFE). Each row represents one (lambda1, ratio) combination.

Usage

extract_multivar_hyperparams(object, fit)

Arguments

Details

The function averages MSFE across cross-validation folds and returns one row per hyperparameter combination. The dataframe can be used for:

• Finding best hyperparameters: df[which.min(df$MSFE), ]

• Plotting MSFE surface

• Filtering/analyzing hyperparameter performance

Value

A data frame with columns:

ratio_index

Index of ratio (1 to n_ratios)

lambda1_index

Index of lambda1 (1 to n_lambda)

slice_index

Flattened index using C++ ordering

lambda1_value

Actual lambda1 penalty value

ratio_value

Actual ratio value

lambda2_value

Computed lambda2 = lambda1 * ratio

MSFE

Mean squared forecast error averaged across CV folds

Description

Get Common Effect Columns

Usage

get_common_cols(spec)

Arguments

Value

Integer vector c(start, end) or NULL

Description

Get Common TVP Columns for a Period

Usage

get_common_tvp_cols(spec, period)

Arguments

Value

Integer vector c(start, end) or NULL

Description

Simplified functions for extracting dynamics matrices from fitted multivar models. These handle both TVP and non-TVP models transparently.

Usage

get_dynamics(
  fit,
  subject = NULL,
  period = NULL,
  type = c("total", "common", "unique")
)

get_common_effects(fit, period = NULL)

get_unique_effects(fit, subject = NULL, period = NULL)

get_total_effects(fit, subject = NULL, period = NULL)

get_dynamics_all_subjects(
  fit,
  period = NULL,
  type = c("total", "common", "unique")
)

Arguments

Details

These helper functions simplify access to dynamics matrices by handling the different structures of TVP and non-TVP models automatically.

**For non-TVP models:** - 'fit$mats$common' is a single matrix (shared across all subjects) - 'fit$mats$unique' is a list of K matrices (one per subject) - 'fit$mats$total' is a list of K matrices (common + unique per subject)

**For TVP models:** - 'fit$mats$common' is a single matrix (time-invariant, shared across all subjects) - 'fit$mats$tvp' is a list of K lists, each containing P period matrices (time-varying per subject) - 'fit$mats$total' is a list of K lists of P matrices (common + tvp per subject per period)

The helper functions hide this complexity from the user.

Examples

## Not run: 
# Non-TVP model
fit <- cv.multivar(object)
get_dynamics(fit, subject = 1)              # Subject 1 total effects
get_common_effects(fit)                     # Common effects
get_dynamics_all_subjects(fit)              # All subjects

# TVP model
fit_tvp <- cv.multivar(object_tvp)
get_dynamics(fit_tvp, subject = 1, period = 2)    # Subject 1, Period 2
get_common_effects(fit_tvp, period = 2)           # Common for Period 2
get_dynamics_all_subjects(fit_tvp, period = 2)    # All subjects, Period 2

## End(Not run)

Description

Returns the number of time-varying periods in a TVP model. Returns 1 for non-TVP models.

Usage

get_n_periods(fit)

Arguments

Value

Integer number of periods

Examples

## Not run: 
n_periods <- get_n_periods(fit)

## End(Not run)

Description

Get Row Range for a Subject's Period

Usage

get_period_rows(spec, subject, period)

Arguments

Value

Integer vector c(start, end) or NULL

Description

Get Subgroup Effect Columns

Usage

get_subgrp_cols(spec, subgroup)

Arguments

Value

Integer vector c(start, end) or NULL

Description

Get Subgroup Columns for a Subject

Usage

get_subgrp_cols_for_subject(spec, subject)

Arguments

Value

Integer vector c(start, end) or NULL

Description

Get Row Range for a Subject

Usage

get_subject_rows(spec, subject)

Arguments

Value

Integer vector c(start, end)

Description

Get Unique Effect Columns for a Subject

Usage

get_unique_cols(spec, subject)

Arguments

Value

Integer vector c(start, end) or NULL (for k=1 TVP)

Description

Get Unique TVP Columns for a Subject and Period

Usage

get_unique_tvp_cols(spec, subject, period)

Arguments

Value

Integer vector c(start, end) or NULL

Description

Returns TRUE if the model has time-varying parameters (TVP).

Usage

is_tvp(fit)

Arguments

Value

Logical indicating if model is TVP

Examples

## Not run: 
if (is_tvp(fit)) {
  cat("Model has time-varying parameters\n")
}

## End(Not run)

Description

Construct the sequence (path) of lambda1 values that will be used for cv.

Usage

lambda_grid(
  depth = 1000,
  nlam = 30,
  Y,
  Z,
  W = NULL,
  tol = 1e-04,
  intercept = FALSE,
  lamadapt = FALSE,
  k = NULL
)

Arguments

Details

In this package, lambda1 is the main sparsity / regularization strength applied to the common effects. Laer, we scaled these into the unique effects via ratios_unique. Importantly, we don't just pick a single lambda1 — we build a path of candidate lambda1 values, from very large (forces almost all penalized coefficients to zero) down to much smaller (less shrinkage).

This path is defined by a decreasing set of values, generated from lambda_max down to lambda_max / depth, so'depth' is used to control that span and the default value is 1,000.

Note: the models have multiple penalty scenarios as alluded to earlier. We call these scenarios, where each scenario encodes a different way to penalize the common and unique effects. Those scenarios live along the 3rd dimension of the array W are typically of length length(ratios_unique)*legnth(lambda1), such that

dim(W) = [ n_outcomes , n_predictors , n_scenarios ]

We build one lambda1 path per scenario. That’s why the return value is a matrix with:

rows = number of lambda1 values along the path columns = number of scenarios

- cv.multivar() stores this matrix in object@lambda1. - During CV, for each scenario column i and each row r in that column, the solver fits the model at that lambda1 and scores prediction error.

- We also have object@ratios_unique, where each ratio says how strong to penalize deviation blocks (individual/subgroup/time-varying) relative to the common effects. Conceptually this means lambda2 = ratio * lambda1.

- We use the ingredients: * lambda1 grid (here) * ratios_unique (constructModel) to construct the lambda2 grid..

Some notes on the arguments:

Value

lambda1 matrix, dimension [nlam x n_scenarios].

Column i is the lambda1 path for scenario i. Row r is a particular lambda1 magnitude along that path

Later: - cv.multivar() stores this in object@lambda1, and cross-validates across [row r, column i]. - Combined with object@ratios_unique, you can reconstruct lambda2 via lambda2 = ratio * lambda1 and plot CV error surfaces.

Description

Provides debiased, robustly-aggregated, analytically-thresholded initial estimates for adaptive LASSO weights. Replaces the default median-based decomposition when weightest = "maity" in constructModel.

References

Maity, S., Sun, Y., & Banerjee, M. (2022). Meta-analysis of heterogeneous data: integrative sparse regression in high-dimensions. JMLR, 23, 1-50.

Description

Standard S3 methods for fitted multivar models.

Usage

## S3 method for class 'multivar_fit'
print(x, ...)

## S3 method for class 'multivar_fit'
summary(object, ...)

## S3 method for class 'multivar_fit'
plot(x, type = c("dynamics", "prevalence"), ...)

## S3 method for class 'multivar_fit'
coef(
  object,
  type = c("total", "common", "unique"),
  subject = NULL,
  period = NULL,
  ...
)

Arguments

Functions

• print(multivar_fit): Print a fitted multivar model.

• summary(multivar_fit): Summarize a fitted multivar model.

• plot(multivar_fit): Plot fitted dynamics or prevalence.

• coef(multivar_fit): Extract coefficient matrices.

Description

Simulate multivar data.

Usage

multivar_sim(
  k,
  d,
  n,
  prop_fill_com,
  prop_fill_ind,
  lb,
  ub,
  sigma,
  unique_overlap = FALSE,
  mat_common = NULL,
  mat_unique = NULL,
  mat_total = NULL,
  diag = FALSE,
  intercept = NULL,
  standardize_output = FALSE
)

Arguments

Examples

k <- 3
d <- 10
n <- 20
prop_fill_com <- .1
prop_fill_ind <- .05
lb <- 0.1
ub <- 0.5
sigma <- diag(d)
data <- multivar_sim(k, d, n, prop_fill_com, prop_fill_ind, lb, ub,sigma)$data

Description

Simulate multivar data with time-varying dynamics according to a latent growth curve.

Usage

multivar_sim_breaks(
  k,
  d,
  n,
  prop_fill_com,
  prop_fill_ind,
  prop_fill_com_growth,
  prop_fill_ind_growth,
  offset = 0,
  lb,
  ub,
  lb_int,
  ub_int,
  lb_slope,
  ub_slope,
  sigma,
  diag = FALSE,
  iters = 1,
  intercept = NULL
)

Arguments

Description

Simulate multivar data with time-varying dynamics according to a latent growth curve.

Usage

multivar_sim_growth(
  k,
  d,
  n,
  prop_fill_com,
  prop_fill_ind,
  prop_fill_com_growth,
  prop_fill_ind_growth,
  offset = 0,
  lb,
  ub,
  lb_int,
  ub_int,
  lb_slope,
  ub_slope,
  lb_quad = 0,
  ub_quad = 0,
  sigma,
  diag = FALSE,
  iters = 1,
  type = "linear",
  extreme_cutoff = 10,
  intercept = NULL
)

Arguments

Description

Generates VAR data with a hierarchical structure: common effects (shared by all), subgroup effects (shared within groups), and unique effects (individual-specific). Importantly, these three components occupy NON-OVERLAPPING positions in the coefficient matrices, ensuring clean identifiability.

Usage

multivar_sim_subgroups(
  k,
  d,
  n,
  subgroup,
  p_com = d/d^2,
  p_sub = 0.05,
  p_ind = 0.05,
  sigma = NULL,
  lb = 0,
  ub = 1,
  intercept = NULL,
  max_attempts = 100,
  max_eig_threshold = 0.99
)

Arguments

Details

The function generates a three-level hierarchical structure where:

Common effects

Shared by all subjects, typically on diagonal positions

Subgroup effects

Shared within subgroups, on distinct off-diagonal positions

Unique effects

Individual-specific, on remaining positions

The key feature is that these three components occupy NON-OVERLAPPING positions in the coefficient matrices, which ensures identifiability of the decomposition:

$total[i] = common + subgroup[g] + unique[i]$

The function repeatedly generates matrices until all subjects have stationary dynamics (maximum absolute eigenvalue < max_eig_threshold).

Value

A list with the following components:

data

List of k data matrices (each n x d)

subgroup

Vector of subgroup memberships

mat_com

Common effect matrix (d x d)

mat_sub_unique

List of k subgroup effect matrices (replicated within groups)

mat_ind_unique

List of k unique effect matrices

mat_ind_final

List of k total effect matrices

intercept

List of intercept vectors (if provided)

Examples

## Not run: 
# Simulate 6 subjects in 2 subgroups
sim <- multivar_sim_subgroups(
  k = 6,
  d = 4,
  n = 100,
  subgroup = c(1, 1, 1, 2, 2, 2),
  p_com = 0.2,
  p_sub = 0.1,
  p_ind = 0.05
)

# Fit model
object <- constructModel(sim$data, subgroup_membership = sim$subgroup)
fit <- cv.multivar(object)

# Evaluate including subgroup effects
perf <- eval_multivar_performance(sim, fit)
print(perf)

## End(Not run)

Description

Simulate multivar TVP data (piecewise or continuous)

Usage

multivar_sim_tvp(
  k,
  d,
  n,
  sigma,
  phi_common_list,
  phi_unique_list = NULL,
  T_per_period = NULL,
  mat_tvp = NULL,
  type = c("piecewise", "continuous"),
  burn = 500,
  intercept = NULL
)

Arguments

Value

List with components:

Examples

# Piecewise TVP with 2 subjects, 2 periods
d <- 3
phi_com1 <- matrix(0, d, d); phi_com1[1,2] <- 0.3
phi_com2 <- matrix(0, d, d); phi_com2[1,2] <- 0.3; phi_com2[2,3] <- 0.2

phi_uniq <- list(
  list(matrix(0, d, d), matrix(0, d, d)),  # subject 1: no unique effects
  list(matrix(0, d, d), matrix(0, d, d))   # subject 2: no unique effects
)
phi_uniq[[2]][[1]][3,1] <- 0.2  # subject 2, period 1 has unique edge

sim <- multivar_sim_tvp(
  k = 2, d = d, n = 200, sigma = diag(d),
  phi_common_list = list(phi_com1, phi_com2),
  phi_unique_list = phi_uniq,
  T_per_period = c(100, 100),
  type = "piecewise"
)

Description

An object class to be used with cv.multivar

Details

To construct an object of class multivar, use the function constructModel

Slots

k

Numeric. The number of subjects (or groupings) in the dataset.

n

Numeric Vector. Vector containing the number of timepoints for each dataset.

d

Numeric Vector. Vector containing the number of variables for each dataset.

Ak

List. A list (length = k) of lagged (T-lag-horizon) by d multivariate time series.

bk

List. A list (length = k) of (T-lag-horizon) by d multivariate time series.

Ak_orig

List. A list (length = k) of lagged (T-lag-horizon) by d unscaled multivariate time series.

bk_orig

List A list (length = k) of (T-lag-horizon) by d unscaled multivariate time series.

Hk

List. A list (length = k) of (horizon) by d multivariate time series.

A

Matrix. A matrix containing the lagged ((T-lag-horizon)k) by (d+dk) multivariate time series.

b

Matrix. A matrix containing the non-lagged ((T-lag-horizon)k) by (d) multivariate time series.

H

Matrix. A matrix containing the non-lagged (horizon k) by d multivariate time series.

data_means

List. A list (length = k) of column means for each group's data, used for intercept recovery when intercept=TRUE.

lag

Numeric. The VAR order. Currently only lag 1 is supported.

horizon

Numeric. Forecast horizon.

t1

Numeric vector. Index of time series in which to start cross validation for individual k.

t2

Numeric vector. Index of time series in which to end cross validation for individual k.

t1k

Numeric vector. Index of time series in which to start cross validation for individual k.

t2k

Numeric vector. Index of time series in which to end cross validation for individual k.

ntk

Numeric. Number of usable timepoints (rows of b) per individual k

ndk

Numeric. Number of variables (cols of b) per individual k

lambda1

Numeric vector. Regularization parameter grid.

nlambda1

Numeric. Number of lambda1 values to search over. Default is 30.

n_ratios_subgroup

Numeric. Number of ratios_subgroup values to search over. Default is 30.

gamma

Numeric. Power for adaptive weighting (default 1). Controls how strongly initial estimates influence the penalty: weight = 1/|coef|^gamma.

tol

Numeric. Convergence tolerance.

depth

Numeric. Depth of grid construction. Default is 1000.

window

Numeric. Size of rolling window.

standardize

Logical. Default is true. Whether to standardize the individual data.

weightest

Character. How to estimate initial coefficients for adaptive weights. Default is "lasso". Other options include "ridge" and "ols".

canonical

Logical. Default is false. If true, individual datasets are fit to a VAR(1) model.

threshold

Logical. Default is false. If true, and canonical is true, individual transition matrices are thresholded based on significance.

lassotype

Character. Default is "adaptive". Choices are "standard" or "adaptive" lasso.

intercept

Logical. Default is FALSE.

W

Matrix. Default is NULL.

ratios_unique

Numeric vector. Penalty ratio for unique effects. Default is NULL.

ratios_subgroup

Numeric vector. Penalty ratio for subgroup effects. Default is NULL.

ratios_unique_tvp

Numeric vector. Default is NULL.

ratios_common_tvp

Numeric vector. Penalty ratio for common TVP effects. Default is NULL.

cv

Character. Default is "blocked" for k-folds blocked cross-validation. rolling window cross-validation also available using "rolling". If "blocked" is selected the nfolds argument should be specified.

nfolds

Numeric. The number of folds for use with "blocked" cross-validation.

lamadapt

Logical. Should the lambdas be calculated adaptively. Default is FALSE.

subgroup_membership

Numeric. Vector of subgroup assignments.

subgroup

Logical. Argument whether to run subgrouping algorithm.

B

Matrix. Default is NULL.

initcoefs

List. A list of initial consistent estimates for the total, subgroup, unique and common effects.

pendiag

Logical. Logical indicating where autoregressive paramaters should be penalized. Default is true.

tvp

Logical.

inittvpcoefs

List. A list of initial tvp estimates.

breaks

List. A list of length K indicating structural breaks in the time series.

lambda_choice

Character. Which lambda to use for initial coefficient estimation ("lambda.min" or "lambda.1se").

common_effects

Logical. Whether to include common effects in TVP models. Only applies when tvp = TRUE.

common_tvp_effects

Logical. Whether to include common TVP effects (shared time-varying patterns) in TVP models. Only applies when tvp = TRUE.

save_beta

Logical. Whether to retain the full beta array in the cv.multivar result. Default is TRUE.

ncores

Numeric. Number of cores for parallel computation. Default is 1.

spec

List. Design matrix specification object created by build_matrix_spec. Single source of truth for column/row indices.

eps

Numeric. FISTA convergence tolerance. Default is 1e-3.

warmstart

Logical. Whether to use warm starts in the FISTA solver. Default is TRUE.

stopping_crit

Character. FISTA convergence criterion. One of "absolute", "relative", or "objective".

selection

Character. Model selection criterion: "cv" (default) for cross-validated MSFE, or "ebic" for Extended BIC.

ebic_gamma

Numeric. EBIC tuning parameter, used when selection = "ebic". Default is 0.5.

weight_type

Character. Adaptive weight function type: "standard" (default) uses 1/|coef|^gamma, "bounded" uses 1/(1+|coef|/tau)^gamma.

maity_opts

List. Options for Maity et al. (2022) MrLasso initial estimation when weightest = "maity". Supports eta (scale parameter for re-descending loss, default 2.0) and thresh_const (threshold multiplier, default 1.0).

See Also

constructModel

Description

Plot CV error over (lambda1, lambda2)

Usage

plot_cv_lambda_grid(fit, object = NULL, log_scale = TRUE, show_best = TRUE)

Arguments

Description

Creates heatmap visualizations of transition matrices from fitted multivar models. Supports plotting common effects, unique (subject-specific) effects, total effects, subgroup effects, and time-varying parameter (TVP) effects.

Usage

plot_results(
  x,
  plot_type = c("common", "unique", "total", "subgrp", "tvp", "common_tvp", "tvp_total"),
  facet_ncol = 3,
  subjects = "all",
  periods = "all",
  ub = 1,
  lb = -1,
  palette = "default",
  show_zeros = FALSE,
  ...
)

Arguments

Value

A ggplot2 object that can be further customized

See Also

plot_sim, plot_transition_mat

Examples

## Not run: 
sim <- multivar_sim(k = 3, d = 5, n = 50, prop_fill_com = 0.2,
                    prop_fill_ind = 0.1, sigma = diag(5))
model <- constructModel(sim$data)
fit <- cv.multivar(model)

# Plot common effects
plot_results(fit, plot_type = "common")

# Plot unique effects for first 2 subjects
plot_results(fit, plot_type = "unique", subjects = 1:2)

# For TVP models
fit_tvp <- cv.multivar(constructModel(data, tvp = TRUE, breaks = list(c(50, 100))))
plot_results(fit_tvp, plot_type = "tvp")  # TVP deviations by period
plot_results(fit_tvp, plot_type = "tvp_total")  # Total effects by period

# Customize the plot
p <- plot_results(fit, plot_type = "total")
p + ggplot2::ggtitle("My Custom Title") + ggplot2::theme_minimal()

## End(Not run)

Description

Creates heatmap visualizations of true transition matrices from multivar simulations. Supports plotting common effects, unique (subject-specific) effects, total effects, and subgroup effects.

Usage

plot_sim(
  x,
  plot_type = c("common", "unique", "total", "subgrp"),
  facet_ncol = 3,
  subjects = "all",
  ub = 1,
  lb = -1,
  palette = "default",
  show_zeros = FALSE,
  ...
)

Arguments

Value

A ggplot2 object that can be further customized

See Also

plot_results, plot_transition_mat

Examples

## Not run: 
# Simulate data
sim <- multivar_sim(k = 3, d = 5, n = 50,
                    prop_fill_com = 0.2, prop_fill_ind = 0.1,
                    lb = 0.1, ub = 0.9, sigma = diag(5))

# Plot true common effects
plot_sim(sim, plot_type = "common")

# Plot true unique effects for subjects 1-2
plot_sim(sim, plot_type = "unique", subjects = 1:2)

# Simulate with subgroups
sim_sub <- multivar_sim_subgroups(k = 4, d = 3, n = 40,
                                  subgroup = c(1, 1, 2, 2),
                                  p_com = 0.25, p_sub = 0.10, p_ind = 0.05,
                                  sigma = diag(3))

# Plot subgroup effects
plot_sim(sim_sub, plot_type = "subgrp")

## End(Not run)

Description

Creates a heatmap visualization of an arbitrary transition matrix. Useful for plotting custom matrices or comparing different estimates.

Usage

plot_transition_mat(
  x,
  title = NULL,
  subtitle = NULL,
  ub = 1,
  lb = -1,
  legend = TRUE,
  dimnames = TRUE,
  palette = "default",
  show_zeros = FALSE,
  ...
)

Arguments

Value

A ggplot2 object that can be further customized

See Also

plot_results, plot_sim

Examples

## Not run: 
# Plot random matrix
plot_transition_mat(matrix(rnorm(25), 5, 5), title = "Random Matrix")

# Plot with custom bounds
mat <- matrix(runif(25, -0.5, 0.5), 5, 5)
plot_transition_mat(mat, title = "Small Coefficients", lb = -0.5, ub = 0.5)

# Without legend or dimension names
plot_transition_mat(mat, legend = FALSE, dimnames = FALSE)

## End(Not run)

Description

Fits the common block using cv.glmnet on a block-diagonal design. Returns A_pre (d x (d+1)); if include_intercept = FALSE, the first column is 0.

Usage

pretrained_var_stage1_glmnet(
  object,
  include_intercept = FALSE,
  standardize = FALSE,
  foldid = NULL,
  lambda = NULL,
  lambest = "min"
)

Arguments

Value

list(A_pre, lambda_best, cvfit, include_intercept)

Description

Prints a dynamics matrix (or list of matrices) with clear row/column labels and formatting options. Handles matrices of different sizes gracefully. Can also accept a fit object directly (e.g., from cv.multivar()).

Usage

print_dynamics(
  dynamics,
  digits = 3,
  zero_char = ".",
  highlight_nonzero = TRUE,
  max_print = 15,
  time_labels = TRUE,
  subject_label = NULL,
  period_labels = NULL
)

Arguments

Value

Invisibly returns the input (for chaining)

Examples

# Simple 4x4 matrix
phi <- matrix(c(0.3, 0, 0.4, 0, 0.4, 0.5, 0, 0.3, 0, 0.3, 0.2, 0, 0, 0, 0, 0.4), 4, 4)
print_dynamics(phi)

# List of matrices
phi_list <- list(phi, phi * 0.8)
print_dynamics(phi_list)

# From fit object (automatic extraction)
# fit <- cv.multivar(object)
# print_dynamics(fit)  # Automatically extracts and prints total effects

# TVP fit object with period labels
# fit_tvp <- cv.multivar(object_tvp)
# print_dynamics(fit_tvp, period_labels = c("Pre", "During", "Post"))

Description

Prints dynamics for TVP models where dynamics vary across time periods. This is a specialized wrapper around print_dynamics() for TVP structures.

Usage

print_dynamics_tvp(dynamics_tvp, period_labels = NULL, ...)

Arguments

Value

Invisibly returns the input (for chaining)

Examples

# Three periods with different dynamics
phi_period1 <- matrix(c(0.3, 0, 0.4, 0, 0.4, 0.5, 0, 0.3, 0, 0.3, 0.2, 0, 0, 0, 0, 0.4), 4, 4)
phi_period2 <- matrix(c(0.4, 0.3, 0, 0.2, 0, 0.6, 0.4, 0, 0.3, 0, 0.3, 0, 0, 0.4, 0, 0.5), 4, 4)
phi_period3 <- matrix(c(0.5, 0.2, 0, 0, 0, 0.4, 0.5, 0, 0.5, 0, 0.4, 0.2, 0, 0.5, 0, 0.3), 4, 4)

# For single subject
print_dynamics_tvp(list(phi_period1, phi_period2, phi_period3))

# With custom period labels
print_dynamics_tvp(list(phi_period1, phi_period2, phi_period3),
                   period_labels = c("Baseline", "Treatment", "Follow-up"))

Description

For models with multiple subjects (K>1), displays a matrix showing the proportion/percentage of subjects that have each edge (non-zero entry).

Usage

print_edge_prevalence(
  fit,
  type = c("proportion", "count"),
  digits = 2,
  zero_char = ".",
  threshold = 1e-10
)

Arguments

Value

Invisibly returns the prevalence matrix

Examples

# fit <- cv.multivar(object)  # K=2 or more subjects
# print_edge_prevalence(fit)

Description

Print method for matrix_spec

Usage

## S3 method for class 'matrix_spec'
print(x, ...)

Arguments

Description

After fitting the model on centered data, this function recovers the intercepts using the formula: c = mean(b) - Phi * mean(A). For multi-group models, it also decomposes intercepts into common and group-specific components.

Usage

recover_intercepts(
  mats,
  data_means,
  k,
  d,
  subgroup = FALSE,
  subgroup_membership = NULL,
  tvp = FALSE,
  breaks = NULL
)

Arguments

Details

For each group k, the total intercept is recovered using: $c_total^(k) = mean(b^(k)) - Phi^(k) * mean(A^(k))$ where $b^{(k)}$ are the outcomes $(Y_t)$ and $A^{(k)}$ are the predictors $(Y_{t-1})$. This is the same approach used by glmnet for intercept recovery.

For multi-group models (k > 1), intercepts are decomposed as: $c = (1/K) * sum_k(c_total^(k))$ # Common intercept $c^(k) = c_total^(k) - c$ # Group-specific intercept

This decomposition ensures the constraint: sum_k(c^(k)) = 0

Note: This formula c = mean(b) - Phi * mean(A) is different from the steady-state formula $c = (I - Phi) * mean(Y)$, which assumes mean(b) = mean(A) = mean(Y). In finite samples they give similar but not identical results. The formula used here is preferred as it matches standard practice in penalized regression (e.g., glmnet).

For TVP models, this function recovers time-varying intercepts c_t for each period. When data_means contains per-period means (mean_b_t, mean_A_t), time-specific intercepts are computed. Otherwise, falls back to base intercepts.

Value

List containing: - intercepts_total: List of total intercepts for each group (c_total^(k)) - intercept_common: Common intercept (c) - intercepts_unique: List of group-specific intercepts (c^(k)) - intercepts_subgrp: List of subgroup-specific intercepts (if applicable) - intercepts_tvp: List of time-varying intercepts for TVP models (if applicable)

Description

Takes the output of cv.multivar() and reselects the best penalty scenario using the Extended Bayesian Information Criterion (EBIC) or standard BIC. This can improve structure recovery (fewer false positives) compared to prediction-based cross-validation, especially when n >> p.

Usage

select_by_ebic(fit, gamma = 0.5)

Arguments

Value

A modified copy of fit with:

The original MSFE results are preserved.

Description

Default show method for an object of class multivar

Usage

## S4 method for signature 'multivar'
show(object)

Arguments

Value

Displays the following information about the multivar object:

• To do.

See Also

constructModel

Description

Prints a summary of a fitted multivar model, including dataset information, model specifications, and key results.

Usage

summary_multivar(fit)

Arguments

Value

Invisibly returns the input object

Examples

# fit <- cv.multivar(object)
# summary_multivar(fit)

Description

Checks internal consistency of a matrix_spec object.

Usage

validate_matrix_spec(spec)

Arguments

Value

TRUE if valid, otherwise throws an error

Description

Simulate a time-varying VAR model (piecewise or continuous)

Usage

var_sim_tvp(
  n,
  sigma,
  phi_list,
  T_per_period = NULL,
  mat_tvp = NULL,
  type = c("piecewise", "continuous"),
  burn = 500,
  intercept = 0
)

Arguments

Value

List with components:

Examples

# Piecewise TVP
phi1 <- matrix(c(0.3, 0, 0, 0.3), 2, 2)
phi2 <- matrix(c(0.3, 0.2, 0, 0.3), 2, 2)
sim <- var_sim_tvp(n = 200, sigma = diag(2), phi_list = list(phi1, phi2),
                   T_per_period = c(100, 100), type = "piecewise")

# Piecewise with period-specific intercepts
sim <- var_sim_tvp(n = 200, sigma = diag(2), phi_list = list(phi1, phi2),
                   T_per_period = c(100, 100), type = "piecewise",
                   intercept = list(c(1, 0), c(2, 1)))