Ridge penalized maximum likelihood estimation of the parameters of the first-order vector auto-regressive model with time-varying covariates, in shorthand VARX(1) model. The VARX(1) model explains the current vector of observations \(\mathbf{Y}_{\ast,t+1}\) by a linear combination of the previous observation endogeneous vector and an exogeneous time-varying covariate: \(\mathbf{Y}_{\ast,t+1} = \mathbf{A} \mathbf{Y}_{\ast,t} + \mathbf{B} \mathbf{X}_{\ast,t+1} + \mathbf{\varepsilon}_{\ast,t+1}\), where \(\mathbf{A}\) and \(\mathbf{B}\) are the lag one autoregression and time-varying regression coefficient matrix, respectively, and \(\mathbf{\varepsilon}_{\ast,t+2}\) the vector of errors (or innovations). The VARX(1)-process is assumed to have mean zero. The experimental design is allowed to be unbalanced.
ridgeVARX1(Y, X, lambdaA=-1, lambdaB=-1, lambdaP=-1, lagX,
targetA=matrix(0, dim(Y)[1], dim(Y)[1]),
targetB=matrix(0, dim(Y)[1], dim(X)[1]),
targetP=matrix(0, dim(Y)[1], dim(Y)[1]), targetPtype="none",
fitAB="ml", zerosA=matrix(nrow=0, ncol=2),
zerosB=matrix(nrow=0, ncol=2), zerosAfit="sparse",
zerosBfit="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)Three-dimensional array containing the response data. The first, second and third dimensions correspond to variates, time and samples, respectively. The data are assumed to be centered covariate-wise.
Three-dimensional array containing the time-varying covariate data. The first, second and third dimensions correspond to covariates, time and samples, respectively. The data are assumed to be centered covariate-wise.
Ridge penalty parameter (positive numeric of length 1) to be used in the estimation of \(\mathbf{A}\), the matrix with autoregression coefficients.
Ridge penalty parameter (positive numeric of length 1) to be used in the estimation of \(\mathbf{B}\), the matrix with regression coefficients of the time-varying covariates stored in array X.
Ridge penalty parameter (positive numeric of length 1) to be used in the estimation of the inverse error covariance matrix (\(\mathbf{\Omega}_{\varepsilon} (=\mathbf{\Sigma_{\varepsilon}^{-1}})\)): the precision matrix of the errors.
Integer, either 0 or 1, specifying whether \(\mathbf{X}_t\) or \(\mathbf{X}_{t-1}\) affects \(\mathbf{Y}_t\), respectively.
Target matrix to which the matrix \(\mathbf{A}\) is to be shrunken.
Target matrix to which the matrix \(\mathbf{B}\) is to be shrunken.
Target matrix to which the in the inverse error covariance matrix, the precision matrix, is to be shrunken.
A character. If fitAB="ml" the parameters \(\mathbf{A}\) and \(\mathbf{B}\) are estimated by (penalized) maximum likelihood. If fitAB="ss" the parameters \(\mathbf{A}\) and \(\mathbf{B}\) are estimated by (penalized) sum of squares.
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.
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.
A matrix with indices of entries of \(\mathbf{B}\) that are constrained to zero. The matrix comprises two columns, each row corresponding to an entry of \(\mathbf{B}\). The first column contains the row indices and the second the column indices.
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.
A character, either "sparse" or "dense". With "sparse", the matrix \(\mathbf{B}\) is assumed to contain many zeros and a computational efficient implementation of its estimation is employed. If "dense", it is assumed that \(\mathbf{B}\) contain only few zeros and the estimation method is optimized computationally accordingly.
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.
A list-object containing the node indices per clique as object from the rip-function.
A list-object containing the node indices per clique as object from the rip-function.
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.
A logical, indicates whether the inverse error covariance matrix is assumed to be diagonal.
A logical, affects estimation of \(\mathbf{A}\). Details below.
Maximum number of iterations (positive numeric of length 1) to be used in maximum likelihood estimation.
Minimum distance (positive numeric of length 1) between estimates of two successive iterations to be achieved.
A list-object with slots:
Ridge ML estimate of the matrix \(\mathbf{A}\), the matrix with lag one autoregression coefficients.
Ridge ML estimate of the matrix \(\mathbf{B}\), the matrix with regression coefficients of the time-varying covariates.
Ridge ML estimate of the inverse error covariance matrix \(\mathbf{\Omega}_{\varepsilon} (=\mathbf{\Sigma_{\varepsilon}^{-1}})\).
Positive numeric of length one: ridge penalty used in the estimation of \(\mathbf{A}\).
Positive numeric of length one: ridge penalty used in the estimation of \(\mathbf{B}\).
Positive numeric of length one: ridge penalty used in the estimation of inverse error covariance matrix \(\mathbf{\Omega}_{\varepsilon} (=\mathbf{\Sigma_{\varepsilon}^{-1}})\).
If diagP=TRUE, no penalization to estimation of the covariance matrix is applied. Consequently, the arguments lambdaP and targetP are ignored (if supplied).
The ridge ML estimator employs the following estimator of the variance of the VARX(1) process:
$$ \frac{1}{n (\mathcal{T} - 1)} \sum_{i=1}^{n} \sum_{t=2}^{\mathcal{T}} \mathbf{Y}_{\ast,i,t} \mathbf{Y}_{\ast,i,t}^{\top}. $$
This is used when efficient=FALSE. However, a more efficient estimator of this variance can be used
$$ \frac{1}{n \mathcal{T}} \sum_{i=1}^{n} \sum_{t=1}^{\mathcal{T}} \mathbf{Y}_{\ast,i,t} \mathbf{Y}_{\ast,i,t}^{\top}, $$
which is achieved by setting when efficient=TRUE. Both estimators are adjusted accordingly when dealing with an unbalanced design.
Miok, V., Wilting, S.M., Van Wieringen, W.N. (2019), ``Ridge estimation of network models from time-course omics data'', Biometrical Journal, 61(2), 391-405.
# NOT RUN {
# set dimensions (p=covariates, n=individuals, T=time points)
p <- 3; n <- 4; T <- 10
# set model parameters
SigmaE <- diag(p)/4
Ax <- createA(p, "chain", nBands=1)
# generate time-varying covariates in accordance with VAR(1) process
X <- dataVAR1(n, T, Ax, SigmaE)
# set model parameters
B <- createA(p, "clique", nCliques=1)
A <- createA(p, "hub", nHubs=1)
# generate time-varying covariates in accordance with VAR(1) process
Y <- dataVARX1(X, A, B, SigmaE, lagX=0)
# fit VARX(1) model
ridgeVARX1(Y, X, 1, 1, 1, lagX=0)
# }
Run the code above in your browser using DataLab