loglikLOOCVcontourVARX1: Contourplot of the LOOCV log-likelihood of VARX(1) model

Description

Evaluates the leave-one-out cross-validated log-likelihood of the VARX(1) model over a grid of the ridge penalty parameters (\(\lambda_a\) and \(\lambda_b\)) for the autoregression and time-varying covariate regression coefficient matrices \(\mathbf{A}\) and \(\mathbf{B}\), respectively, while keeping \(\lambda_{\omega}\), the penalty parameter of the inverse error covariance matrix \(\mathbf{\Omega}_{\varepsilon} (=\mathbf{\Sigma_{\varepsilon}^{-1}})\), fixed at a user-specified value. The result is plotted as a contour plot, which facilitates the choice of optimal penalty parameters. The function also works with a (possibly) unbalanced experimental set-up. The VARX(1)-process is assumed to have mean zero.

Usage

loglikLOOCVcontourVARX1(lambdaAgrid, lambdaBgrid, lambdaPgrid, Y, X,
                        lagX=0, figure=TRUE, verbose=TRUE, ...)

Arguments

lambdaAgrid

A numeric of length larger than one, comprising positive numbers only. It contains the grid points corresponding to the \(\lambda_a\) (the penalty parameter for the lag one autoregression coefficient matrix \(\mathbf{A}\)).

lambdaBgrid

A numeric of length larger than one, comprising positive numbers only. It contains the grid points corresponding to the \(\lambda_b\) (the penalty parameters for the regression coefficient matrix \(\mathbf{B}\) of the time-varying covariates.

lambdaPgrid

A numeric of length larger than one, comprising positive numbers only. It contains the grid points corresponding to the \(\lambda_{\omega}\) (the penalty parameters for the inverse error covariance matrix \(\mathbf{\Omega}_{\varepsilon} (=\mathbf{\Sigma_{\varepsilon}^{-1}})\)).

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 variate-wise.

Three-dimensional array containing the time-varying covariates. The first, second and third dimensions correspond to covariates, time and samples, respectively. The data are assumed to be centered covariate-wise.

lagX

An integer, either 0 or 1, specifying whether \(\mathbf{X}_t\) or \(\mathbf{X}_{t-1}\) affects \(\mathbf{Y}_t\), respectively.

figure

A logical indicating whether the contour plot should be generated.

verbose

A logical indicator: should intermediate output be printed on the screen?

...

Other arguments to be passed on (indirectly) to ridgeVARX1.

Value

A list-object with slots:

lambdaA

A numeric with the grid points corresponding to \(\lambda_a\) (the penalty parameter for the autoregression coefficient matrix \(\mathbf{A}\)).

lambdaB

A numeric with the grid points corresponding to \(\lambda_b\) (the penalty parameter for time-varying covariate regression coefficient matrix \(\mathbf{B}\).

lambdaP

A numeric with the grid points corresponding to \(\lambda_{\omega}\) (the penalty parameter for the inverse error covariance matrix \(\mathbf{\Omega}_{\varepsilon} (=\mathbf{\Sigma_{\varepsilon}^{-1}})\).

llLOOCV

A matrix of leave-one-out cross-validated log-likelihoods. Rows and columns correspond to \(\lambda_a\) and \(\lambda_{b}\) values, respectively.

References

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.

Examples

Run this code

# NOT RUN {
# set dimensions (p=covariates, n=individuals, T=time points)
p <- 3; n <- 12; T <- 10

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

# generate time-varying covariate data
X <- dataVAR1(n, T, Ax, SigmaE)

# regression parameter matrices of VARX(1) model
A <- createA(p, topology="clique", nonzeroA=0.1, nClique=1)
B <- createA(p, topology="hub", nonzeroA=0.1, nHubs=1)

# generate data
Y <- dataVARX1(X, A, B, SigmaE, lagX=0)

## plot contour of cross-validated likelihood
# }
# NOT RUN {
 lambdaAgrid <- seq(0.01, 1, length.out=20) 
# }
# NOT RUN {
 lambdaBgrid <- seq(0.01, 1000, length.out=20) 
# }
# NOT RUN {
 lambdaPgrid <- seq(0.01, 1000, length.out=20) 
# }
# NOT RUN {
 loglikLOOCVcontourVARX1(lambdaAgrid, lambdaBgrid, lambdaPgrid, Y, X) 
# }
# NOT RUN {
## determine optimal values of the penalty parameters
# }
# NOT RUN {
optLambdas <- constrOptim(c(1,1,1), loglikLOOCVVARX1, gr=NULL, 
# }
# NOT RUN {
              ui=diag(3), ci=c(0,0,0), Y=Y, X=X,
# }
# NOT RUN {
              control=list(reltol=0.01))$par 
# }
# NOT RUN {
## add point of optimum
# }
# NOT RUN {
 points(optLambdas[1], optLambdas[2], pch=20, cex=2, 
# }
# NOT RUN {
 col="red") 
# }