mrce {MRCE}

Do multivariate regression with covariance estimation (MRCE)

Description

Let S_{+}^q be the set of q by q symmetric and positive definite matrices and let y_i\in R^q be the measurements of the q responses for the ith subject (i=1,\ldots, n). The model assumes that y_i is a realization of the q-variate random vector

Y_i = \mu + \beta'x_i + \varepsilon_i, \ \ \ \ i=1,\ldots, n

where \mu\in R^q is an unknown intercept vector; \beta\in R^{p\times q} is an unknown regression coefficient matrix; x_i \in R^p is the known vector of values for ith subjects's predictors, and \varepsilon_1,\ldots, \varepsilon_n are n independent copies of a q-variate Normal random vector with mean 0 and unknown inverse covariance matrix \Omega \in S_{+}^q.

This function computes penalized likelihood estimates of the unknown parameters \mu, \beta, and \Omega. Let \bar y=n^{-1} \sum_{i=1}^n y_i and \bar{x} = n^{-1}\sum_{i=1}^n x_i. These estimates are

(\hat{\beta}, \hat\Omega) = \arg\min_{(B, Q)\in R^{p\times q}\times S_{+}^q} \left\{g(B, Q) +\lambda_1 \left(\sum_{j\neq k} |Q_{jk}| + 1(p\geq n) \sum_{j=1}^q |Q_{jj}| \right) + 2\lambda_{2}\sum_{j=1}^p\sum_{k=1}^q |B_{jk}|\right\}

and \hat\mu=\bar y - \hat\beta'\bar x, where

g(B, Q) = {\rm tr}\{n^{-1}(Y-XB)'(Y-XB) Q\}-\log|Q|,

Y\in R^{n\times q} has ith row (y_{i}-\bar y)', and X\in R^{n\times p} has ith row (x_{i}-\bar{x})'.

Usage

mrce(X,Y, lam1=NULL, lam2=NULL, lam1.vec=NULL, lam2.vec=NULL,
     method=c("single", "cv", "fixed.omega"),
     cov.tol=1e-4, cov.maxit=1e3, omega=NULL, 
     maxit.out=1e3, maxit.in=1e3, tol.out=1e-8, 
     tol.in=1e-8, kfold=5, silent=TRUE, eps=1e-5, 
     standardize=FALSE, permute=FALSE)

Arguments

X

An n by p matrix of the values for the prediction variables. The ith row of X is x_i defined above (i=1,\ldots, n). Do not include a column of ones.

Y

An n by q matrix of the observed responses. The ith row of Y is y_i defined above (i=1,\ldots, n).

lam1

A single value for \lambda_1 defined above. This argument is only used if method="single"

lam2

A single value for \lambda_2 defined above (or a p by q matrix with (j,k)th entry \lambda_{2jk} in which case the penalty 2\lambda_{2}\sum_{j=1}^p\sum_{k=1}^q |B_{jk}| becomes 2\sum_{j=1}^p\sum_{k=1}^q \lambda_{2jk}|B_{jk}|). This argument is not used if method="cv".

lam1.vec

A vector of candidate values for \lambda_1 from which the cross validation procedure searches: only used when method="cv" and must be specified by the user when method="cv". Please arrange in decreasing order.

lam2.vec

A vector of candidate values for \lambda_2 from which the cross validation procedure searches: only used when method="cv" and must be specified by the user when method="cv". Please arrange in decreasing order.

method

There are three options:

• method="single" computes the MRCE estimate of the regression coefficient matrix with penalty tuning parameters lam1 and lam2;

• method="cv" performs kfold cross validation using candidate tuning parameters in lam1.vec and lam2.vec;

• method="fixed.omega" computes the regression coefficient matrix estimate for which Q (defined above) is fixed at omega.

cov.tol

Convergence tolerance for the glasso algorithm that minimizes the objective function (defined above) with B fixed.

cov.maxit

The maximum number of iterations allowed for the glasso algorithm that minimizes the objective function (defined above) with B fixed.

omega

A user-supplied fixed value of Q. Only used when method="fixed.omega" in which case the minimizer of the objective function (defined above) with Q fixed at omega is returned.

maxit.out

The maximum number of iterations allowed for the outer loop of the exact MRCE algorithm.

maxit.in

The maximum number of iterations allowed for the algorithm that minimizes the objective function, defined above, with \Omega fixed.

tol.out

Convergence tolerance for outer loop of the exact MRCE algorithm.

tol.in

Convergence tolerance for the algorithm that minimizes the objective function, defined above, with \Omega fixed.

kfold

The number of folds to use when method="cv".

silent

Logical: when silent=FALSE this function displays progress updates to the screen.

eps

The algorithm will terminate if the minimum diagonal entry of the current iterate's residual sample covariance is less than eps. This may need adjustment depending on the scales of the variables.

standardize

Logical: should the columns of X be standardized so each has unit length and zero average. The parameter estimates are returned on the original unstandarized scale. The default is FALSE.

permute

Logical: when method="cv", should the subject indices be permutted? The default is FALSE.

Details

Please see Rothman, Levina, and Zhu (2010) for more information on the algorithm and model. This version of the software uses the glasso algorithm (Friedman et al., 2008) through the R package glasso. If the algorithm is running slowly, track its progress with silent=FALSE. In some cases, choosing cov.tol=0.1 and tol.out=1e-10 allows the algorithm to make faster progress. If one uses a matrix for lam2, consider setting tol.in=1e-12.

When p \geq n, the diagonal of the optimization variable corresponding to the inverse covariance matrix of the error is penalized. Without diagonal penalization, if there exists a \bar B such that the qth column of Y is equal to the qth column of X\bar B, then a global minimizer of the objective function (defined above) does not exist.

The algorithm that minimizes the objective function, defined above, with Q fixed uses a similar update strategy and termination criterion to those used by Friedman et al. (2010) in the corresponding R package glmnet.

Value

A list containing

Bhat

This is \hat\beta \in R^{p\times q} defined above. If method="cv", then best.lam1 and best.lam2 defined below are used for \lambda_1 and \lambda_2.

muhat

This is the intercept estimate \hat\mu \in R^q defined above. If method="cv", then best.lam1 and best.lam2 defined below are used for \lambda_1 and \lambda_2.

omega

This is \hat\Omega \in S_{+}^q defined above. If method="cv", then best.lam1 and best.lam2 defined below are used for \lambda_1 and \lambda_2.

mx

This is \bar x \in R^p defined above.

my

This is \bar y \in R^q defined above.

best.lam1

The selected value for \lambda_1 by cross validation. Will be NULL unless method="cv".

best.lam2

The selected value for \lambda_2 by cross validation. Will be NULL unless method="cv".

cv.err

Cross validation error matrix with length(lam1.vec) rows and length(lam2.vec) columns. Will be NULL unless method="cv".

Note

The algorithm is fastest when \lambda_1 and \lambda_2 are large. Use silent=FALSE to check if the algorithm is converging before the total iterations exceeds maxit.out.

Author(s)

Adam J. Rothman

References

Rothman, A. J., Levina, E., and Zhu, J. (2010) Sparse multivariate regression with covariance estimation. Journal of Computational and Graphical Statistics. 19: 947–962.

Jerome Friedman, Trevor Hastie, Robert Tibshirani (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3), 432-441.

Jerome Friedman, Trevor Hastie, Robert Tibshirani (2010). Regularization Paths for Generalized Linear Models via Coordinate Descent. Journal of Statistical Software, 33(1), 1-22.

Examples

set.seed(48105)
n=50
p=10
q=5

Omega.inv=diag(q)
for(i in 1:q) for(j in 1:q)
  Omega.inv[i,j]=0.7^abs(i-j)
out=eigen(Omega.inv, symmetric=TRUE)
Omega.inv.sqrt=tcrossprod(out$vec*rep(out$val^(0.5), each=q),out$vec)
Omega=tcrossprod(out$vec*rep(out$val^(-1), each=q),out$vec)

X=matrix(rnorm(n*p), nrow=n, ncol=p)
E=matrix(rnorm(n*q), nrow=n, ncol=q)%*%Omega.inv.sqrt
Beta=matrix(rbinom(p*q, size=1, prob=0.1)*runif(p*q, min=1, max=2), nrow=p, ncol=q)
mu=1:q

Y=rep(1,n)%*%t(mu) + X%*%Beta + E

lam1.vec=rev(10^seq(from=-2, to=0, by=0.5))
lam2.vec=rev(10^seq(from=-2, to=0, by=0.5))
cvfit=mrce(Y=Y, X=X, lam1.vec=lam1.vec, lam2.vec=lam2.vec, method="cv")
cvfit

fit=mrce(Y=Y, X=X, lam1=10^(-1.5), lam2=10^(-0.5), method="single")
fit

lam2.mat=1000*(fit$Bhat==0)
refit=mrce(Y=Y, X=X, lam2=lam2.mat, method="fixed.omega", omega=fit$omega, tol.in=1e-12) 
refit