Package 'MRCE'

Multivariate regression with covariance estimation

Description

Computes the MRCE estimators (Rothman, Levina, and Zhu, 2010) and has the dataset stock04 used in Rothman, Levina, and Zhu (2010), originally analyzed in Yuan et al. (2007).

Details

The primary function is mrce. The dataset is stock04.

Author(s)

Adam J. Rothman

Maintainer: Adam J. Rothman <[email protected]>

References

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

Yuan, M., Ekici, A., Lu, Z., and Monteiro, R. (2007). Dimension reduction and coefficient estimation in multivariate linear regression. Journal of the Royal Statistical Society Series B 69(3):329–346.

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.

---

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 $i$th 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 $i$th 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 $i$th row $(y_{i}-\bar y)'$, and $X\in R^{n\times p}$ has $i$th 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 $i$th 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 $i$th 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 $q$th column of $Y$ is equal to the $q$th 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

---

log-returns of 9 stocks from 2004

Description

Weekly log-returns of 9 stocks from 2004, analyzed in Yuan et al. (2007)

Usage

data(stock04)

Format

The format is: num [1:52, 1:9] 0.002275 -0.003795 0.012845 0.017489 -0.000369 ... - attr(*, "dimnames")=List of 2 ..$ : NULL ..$ : chr [1:9] "Walmart" "Exxon" "GM" "Ford" ...

Source

Yuan, M., Ekici, A., Lu, Z., and Monteiro, R. (2007). Dimension reduction and coefficient estimation in multivariate linear regression. Journal of the Royal Statistical Society Series B, 69(3):329–346.

References

Yuan, M., Ekici, A., Lu, Z., and Monteiro, R. (2007). Dimension reduction and coefficient estimation in multivariate linear regression. Journal of the Royal Statistical Society Series B, 69(3):329–346.

---