ridgeVAR1: Ridge ML estimation of the VAR(1) model
In wvanwie/ragt2ridges: Ridge Estimation of Vector Auto-Regressive (VAR) Processes

Description Usage Arguments Details Value Note Author(s) References See Also Examples

Ridge penalized maximum likelihood estimation of the parameters of the VAR(1), first-order vector auto-regressive, model. The VAR(1) model explains the current vector of observations \mathbf{Y}_{\ast,t+1} by a linear combination of the previous observation vector: \mathbf{Y}_{\ast,t+1} = \mathbf{A} \mathbf{Y}_{\ast,t} + \mathbf{\varepsilon}_{\ast,t+1}, where \mathbf{A} is the autoregression coefficient matrix and \mathbf{\varepsilon}_{\ast,t+1} the vector of errors (or innovations). The VAR(1)-process is assumed to have mean zero. The experimental design is allowed to be unbalanced.

ridgeVAR1(Y, lambdaA=0, lambdaP=0, 
          targetA=matrix(0, dim(Y)[1], dim(Y)[1]), 
          targetP=matrix(0, dim(Y)[1], dim(Y)[1]), targetPtype="none",
          fitA="ml", zerosA=matrix(nrow=0, ncol=2), 
          zerosAfit="sparse", zerosP=matrix(nrow=0, ncol=2), 
          cliquesP=list(), separatorsP=list(), 
          unbalanced=matrix(nrow=0, ncol=2), diagP=FALSE, 
          efficient=TRUE, nInit=100, minSuccDiff=0.001)

`Y`	Three-dimensional `array` containing the data. The first, second and third dimensions correspond to covariates, time and samples, respectively. The data are assumed to be centered covariate-wise.
`lambdaA`	Ridge penalty parameter (positive `numeric` of length 1) to be used in the estimation of \mathbf{A}, the matrix with autoregression coefficients.
`lambdaP`	Ridge penalty parameter (positive `numeric` of length 1) to be used in the estimation of the inverse error covariance matrix (\mathbf{Ω}_{\varepsilon} (=\mathbf{Σ_{\varepsilon}^{-1}})): the precision matrix of the errors.
`targetA`	Target `matrix` to which the matrix \mathbf{A} is to be shrunken.
`targetP`	Target `matrix` to which the in the inverse error covariance matrix, the precision matrix, is to be shrunken.
`fitA`	A `character`. If `fitA="ml"` the parameter \mathbf{A} is estimate by (penalized) maximum likelihood. If `fitA="ss"` the parameter \mathbf{A} is estimate by (penalized) sum of squares. The latter being much faster as it need not iterate.
`targetPtype`	A `character` indicating the type of target to be used for the precision matrix. When specified it overrules the `targetP`-option. See the `default.target`-function for the options.
`zerosA`	A `matrix` with indices of entries of \mathbf{A} that are constrained to zero. The matrix comprises two columns, each row corresponding to an entry of \mathbf{A}. The first column contains the row indices and the second the column indices.
`zerosAfit`	A `character`, either "sparse" or "dense". With "sparse", the matrix \mathbf{A} is assumed to contain many zeros and a computational efficient implementation of its estimation is employed. If "dense", it is assumed that \mathbf{A} contains only few zeros and the estimation method is optimized computationally accordingly.
`zerosP`	A `matrix`-object with indices of entries of the precision matrix that are constrained to zero. The matrix comprises two columns, each row corresponding to an entry of the adjacency matrix. The first column contains the row indices and the second the column indices. The specified graph should be undirected and decomposable. If not, it is symmetrized and triangulated (unless `cliquesP` and `seperatorsP` are supplied). Hence, the employed zero structure may differ from the input `zerosP`.
`cliquesP`	A `list`-object containing the node indices per clique as object from the `rip`-function.
`separatorsP`	A `list`-object containing the node indices per clique as object from the `rip`-function.
`unbalanced`	A `matrix` with two columns, indicating the unbalances in the design. Each row represents a missing design point in the (time x individual)-layout. The first and second column indicate the time and individual (respectively) specifics of the missing design point.
`diagP`	A `logical`, indicates whether the inverse error covariance matrix is assumed to be diagonal.
`efficient`	A `logical`, affects estimation of \mathbf{A}. Details below.
`nInit`	Maximum number of iterations (positive `numeric` of length 1) to be used in maximum likelihood estimation.
`minSuccDiff`	Minimum distance (positive `numeric` of length 1) between estimates of two successive iterations to be achieved.

The ridge ML estimator employs the following estimator of the variance of the VAR(1) process:

\frac{1}{n (\mathcal{T} - 1)} ∑_{i=1}^{n} ∑_{t=2}^{\mathcal{T}} \mathbf{Y}_{\ast,i,t} \mathbf{Y}_{\ast,i,t}^{\mathrm{T}}.

This is used when efficient=FALSE. However, a more efficient estimator of this variance can be used

\frac{1}{n \mathcal{T}} ∑_{i=1}^{n} ∑_{t=1}^{\mathcal{T}} \mathbf{Y}_{\ast,i,t} \mathbf{Y}_{\ast,i,t}^{\mathrm{T}},

which is achieved by setting when efficient=TRUE. Both estimators are adjusted accordingly when dealing with an unbalanced design.

A list-object with slots:

`A`	Ridge ML estimate of the matrix \mathbf{A}, the `matrix` with lag one auto-regressive coefficients.
`P`	Ridge ML estimate of the inverse error covariance `matrix` \mathbf{Ω}_{\varepsilon} (=\mathbf{Σ_{\varepsilon}^{-1}}).
`lambdaA`	Positive `numeric` of length one: ridge penalty used in the estimation of \mathbf{A}.
`lambdaP`	Positive `numeric` of length one: ridge penalty used in the estimation of inverse error covariance matrix \mathbf{Ω}_{\varepsilon} (=\mathbf{Σ_{\varepsilon}^{-1}}).

When the target of the precision matrix is specified through the targetPtype-argument, the target is data-driven (for both fitA="ss" and fitA="ml"). In particular, it is updated at each iteration when fitA="ml".

Wessel N. van Wieringen <w.vanwieringen@vumc.nl>

Miok, V., Wilting, S.M., Van Wieringen, W.N. (2017), “Ridge estimation of the VAR(1) model and its time series chain graph from multivariate time-course omics data”, Biometrical Journal, 59(1), 172-191.

loglikLOOCVVAR1, ridgeP, default.target, ridgePchordal.

# set dimensions (p=covariates, n=individuals, T=time points)
p <- 3; n <- 4; T <- 10

# set model parameters
SigmaE <- diag(p)/4
A      <- createA(p, "chain")

# generate data
Y <- dataVAR1(n, T, A, SigmaE)

# fit VAR(1) model
ridgeVAR1(Y, 1, 1)$A