cv.SiER {SiER}R Documentation

Cross-validation for high-dimensional multivariate regression

Description

Conduct the cross-validation and build the final model for the following high dimensional multivariate regression model:

Y=μ+Xβ+ϵ,Y= \mu+X\beta+\epsilon,

where YY is the n×qn\times q response matrix with q1q\ge 1, XX is the n×pn\times p predictor matrix, and ϵ\epsilon is the noise matrix. The coefficient matrix β\beta is p×qp\times q and μ\mu is the intercept. The number of predictors pp can be much larger than the sample size nn. The response is univariate if q=1q=1 and multivariate if q>1q>1.

Usage

cv.SiER(X, Y, K.cv = 5, upper.comp = 10, thresh = 0.01)

Arguments

X

the n×pn\times p predictor matrix.

Y

the n×qn\times q response matrix, where q1q\ge 1 is the number of response variables.

K.cv

the number of CV sets. Default is 5.

upper.comp

the upper bound for the maximum number of components to be calculated. Default is 10.

thresh

a number between 0 and 1 specifying the minimum proportion of variation to be explained by each selected component relative to all the selected components. It is used to determine the maximum number of components to be calculated in the CV procedure. The optimal number of components will be selected from the integers from 1 to the minimum of upper.comp and this maximum number. A smaller thresh leads to a larger maximum number of components and a longer running time. A larger thresh value leads to a smaller running time, but may miss some important components and lead to a larger prediction error. Default is 0.01.

Details

Based on the best rank KK approximation to XβX\beta, the coefficient matrix has decomposition β=αkwkT\beta=\sum \alpha_k w_k ^T, where αk\alpha_k is the vector so that XαkX\alpha_k has the maximum correlation with YY under the restriction that XαkX\alpha_k has unit variance and is uncorrelated with Xα1X\alpha_1,..., Xαk1X\alpha_{k-1}. We estimate αk\alpha_k by solving a penalized generalized eigenvalue problem with penalty ταkλ2\tau||\alpha_k||_{\lambda}^2 where αkλ2=(1λ)αk22+λαk12||\alpha_k||_{\lambda}^2=(1-\lambda)||\alpha_k||_2^2+\lambda||\alpha_k||_1^2 is a mixture of the squared l2l_2 and squared l1l_1 norms. The wkw_k is estimated by regressing YY on XαkX\alpha_k.

Value

A fitted CV-object, which is used in the function pred.SiER() for prediction and getcoef.SiER() for extracting the estimated coefficient matrix.

mu

the estimated intercept vector.

beta

the estimated slope coefficient matrix.

min.error

minimum CV error.

scale.x

the maximum absolute value of X used to scale X.

X

the input X.

Y

the input Y.

params.set

a 9*2 matrix specifying the set of values of tautau and lambdalambda used in CV.

error

a list for CV errors.

opt.K

optimal number of components to be selected.

opt.tau

optimal value for tautau.

opt.lambda

optimal value for lambdalambda.

Author(s)

Ruiyan Luo and Xin Qi

References

Ruiyan Luo and Xin Qi (2017) Signal extraction approach for sparse multivariate response regression, Journal of Multivariate Statistics. 153: 83-97.

Examples

# q=1
library(MASS)
nvar=100
nvarq <- 1
sigmaY <- 0.1
sigmaX=0.1
nvar.eff=15
rho <- 0.3
Sigma=matrix(0,nvar.eff,nvar.eff)
for(i in 1:nvar.eff){
    for(j in 1:nvar.eff){
        Sigma[i,j]=rho^(abs(i-j))
    }
}

betas.true <- matrix(0, nvar, 1)
betas.true[1:15,1]=rep(1,15)/sqrt(15)

ntest <- 100
ntrain <- 90
ntot <- ntest+ntrain
X <- matrix(0,ntot,nvar)
X[,1:nvar.eff] <- mvrnorm(n=ntot, rep(0, nvar.eff), Sigma)
X[,-(1:nvar.eff)] <- matrix(sigmaX*rnorm((nvar-nvar.eff)*dim(X)[1]),
                            dim(X)[1],(nvar-nvar.eff))
Y <- X%*%betas.true
Y <- Y+rnorm(ntot, 0, sigmaY)

X.train <- X[1:ntrain,]
Y.train <- Y[1:ntrain,]
X.test <- X[-(1:ntrain),]
Y.test <- Y[-(1:ntrain),]

cv.fit <- cv.SiER(X.train,Y.train, K.cv=5)

Y.pred=pred.SiER(cv.fit, X.test)
error=sum((Y.pred-Y.test)^2)/ntest
print(c("predict error=", error))
coefs=getcoef.SiER(cv.fit)


#q>1
library(MASS)
total.noise <- 0.1
rho <- 0.3
rho.e <- 0.2
nvar=500
nvarq <- 3
sigma2 <- total.noise/nvarq
sigmaX=0.1
nvar.eff=150

Sigma=matrix(0,nvar.eff,nvar.eff)
for(i in 1:nvar.eff){
    for(j in 1:nvar.eff){
        Sigma[i,j]=rho^(abs(i-j))
    }
}
Sigma2.y <- matrix(sigma2*rho.e,nvarq, nvarq)
diag(Sigma2.y) <- sigma2

betas.true <- matrix(0, nvar, 3)
betas.true[1:15,1]=rep(1,15)/sqrt(15)
betas.true[16:45,2]=rep(0.5,30)/sqrt(30)
betas.true[46:105,3]=rep(0.25,60)/sqrt(60)

ntest <- 500
ntrain <- 90
ntot <- ntest+ntrain
X <- matrix(0,ntot,nvar)
X[,1:nvar.eff] <- mvrnorm(n=ntot, rep(0, nvar.eff), Sigma)
X[,-(1:nvar.eff)] <- matrix(sigmaX*rnorm((nvar-nvar.eff)*dim(X)[1]),
                           dim(X)[1],(nvar-nvar.eff))
Y <- X%*%betas.true
Y <- Y+mvrnorm(n=ntot, rep(0,nvarq), Sigma2.y)

X.train <- X[1:ntrain,]
Y.train <- Y[1:ntrain,]
X.test <- X[-(1:ntrain),]
Y.test <- Y[-(1:ntrain),]

cv.fit <- cv.SiER(X.train,Y.train, K.cv=5)

Y.pred=pred.SiER(cv.fit, X.test)
error=sum((Y.pred-Y.test)^2)/ntest
print(c("predict error=", error))


[Package SiER version 0.1.0 Index]