Simultaneous envelope (SimulEnv)

Here we explain the basics of fitting simultaneous envelope model.

Contents

Load the package and data sets

clear all;
cd D:\EnvelopeComputing\EnvelopeMLM; % change directory to the package folder
setpaths;   % load all the functions in the package
rng(2016)   % set random seed

Simulate a data set

We generate a data set from the multivariate linear model, where (X,Y) are simulated from a joint multivariate normal distribution

p = 10; % number of predictors
r = 15; % number of responses
dx = 3; % dimension of X-envelope
dy = 3; % dimension of Y-envelope
rr = 3; % rank of the regression coefficient
N = 600; % number of observations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Envelopes model parameters          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
trueG = orth(rand(p,dx));   % X-envelope basis
trueG0 = null(trueG');
trueH = orth(rand(r,dy));   % Y-envelope basis
trueH0 = null(trueH');
Eta = rand(dx,dy);          % regression coefficient matrix of reduced Y on reduced X

Omega = eye(dx);    % material variation in X
Phi = eye(dy);      % material variation of Y|X
Omega0 = 0.1*eye(p-dx);    % immaterial variation in X
Phi0 = 10*eye(r-dy);       % immaterial variation of Y|X

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Regression coefficient and covariance matrices      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta = trueH*Eta'*trueG';
sigmaXY = trueG*Omega*Eta*trueH';
sigmaX = trueG*Omega*trueG' + trueG0*Omega0*trueG0';
sigmaY = trueH*(Phi+Eta'*Omega*Eta)*trueH' + trueH0*Phi0*trueH0';
sigmaC = [sigmaX, sigmaXY; sigmaXY', sigmaY];
sigmaD = blkdiag(sigmaX,sigmaY);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Data generation          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
datavec = mvnrnd(zeros(p+r,1),sigmaC,N);
X = datavec(:,1:p);
Y = datavec(:,(p+1):(p+r));

Fit the ordinary least squares (OLS) estimator

Fit the standard model with ordinary least squares. The output is just the regression coefficient matrix, since we are not interested in the slope.

betahat_ols = OLS(X,Y); % r-by-p regression coefficient matrix
% Frobenius norm of true beta minus estimator
disp(['OLS estimation error is: ', num2str(norm(beta-betahat_ols,'fro'))]);

OLS estimation error is: 3.8747

Fit the response envelope (Y-env) estimator

Fit the response envelope model (Y-env) by Cook, Li and Chiaromonte (2010; Statistica Sinica), for reducing the multivariate response Y in multivariate linear model. We offer two options for obtaining the response envelope estimator: FG = 0 or 1 indicating either the 1D algortihm estimation (Cook and Zhang 2016; JCGS) or the Full Grassmannian Optimization, which uses the 1D algorithm to get initial value.

FG = 1; % Full Grassmannian Optimization, the MLE under normality assumption
dy = 3; % specify the Y-envelope dimension
out_yenv = YEnv(X,Y,dy,FG);
% output: an estimated envelope basis matrix (r-by-dy), and the estimated
% regresson coefficient matrix beta (r-by-p)
disp(out_yenv);
betahat_yenv = out_yenv.beta;
% Frobenius norm of true beta minus estimator
disp(['Response envelope (FG) estimation error is: ', num2str(norm(beta-betahat_yenv,'fro'))]);

FG = 0; % switch to the 1D algorithm estimator (faster but not MLE)
out_yenv = YEnv(X,Y,dy,FG);
betahat_yenv = out_yenv.beta;
% Frobenius norm of true beta minus estimator
disp(['Response envelope (1D) estimation error is: ', num2str(norm(beta-betahat_yenv,'fro'))]);

    Gyhat: [15x3 double]
     beta: [15x10 double]

Response envelope (FG) estimation error is: 1.6447
Response envelope (1D) estimation error is: 1.6444

Fit the predictor envelope (X-env) estimator

Fit the predictor envelope model (X-env) by Cook, Helland and Su (2013), for reducing the multivariate predictor vector in multivariate linear models. Similar to the response envelope estimation, we also implemented both FG and 1D estimators.

FG = 1; % FG estimation
dx = 3; % specify the X-envelope dimension
out_xenv = XEnv(X,Y,dx,FG);
betahat_xenv = out_xenv.beta;
disp(['Predictor envelope (FG) estimation error is: ', num2str(norm(beta-betahat_xenv,'fro'))]);

FG = 0; % 1D estimation
out_xenv = XEnv(X,Y,dx,FG);
betahat_xenv = out_xenv.beta;
disp(['Predictor envelope (1D) estimation error is: ', num2str(norm(beta-betahat_xenv,'fro'))]);

Predictor envelope (FG) estimation error is: 2.2712
Predictor envelope (1D) estimation error is: 2.2866

Fit the simultaneous envelope (SimulEnv) estimator

Fit the simultaneous envelope model (SimulEnv) by Cook and Zhang (2015; Technometrics), for simultaneously reducing the predictor and response. We recommend using the 1D algorithm estimator as the final estimator. An alternative is to use full Grassmannian optimization with 1D algorithm as initial value, similar to the X- and Y-envelope models.

FG = 1; % FG estimation
dx=3; dy=3; % specify the simultaneous envelope dimension
out_senv = SimulEnv(X,Y,dx,dy,FG);
betahat_senv = out_senv.beta;
disp(['Simultaneous envelope (FG) estimation error is: ', num2str(norm(beta-betahat_senv,'fro'))]);

FG = 0;
out_senv = SimulEnv(X,Y,dx,dy,FG);
betahat_senv = out_senv.beta;
disp(['Simultaneous envelope (1D) estimation error is: ', num2str(norm(beta-betahat_senv,'fro'))]);

Simultaneous envelope (FG) estimation error is: 0.52822
Simultaneous envelope (1D) estimation error is: 0.57601

Envelope dimensions determined by 1D algorithm

Simultaneous envelope dimensions selected by AIC, BIC and likelihood-ratio test (LRT). Because both the X- and Y-envelope dimensions are greater than or equal to the rank of the regression coefficient matrix beta, so if the rank of beta is know, we can set it as the lower bound for the envelope dimensions. To speed up the computation, we use 1D algortihm estimators as an approximation, instead of the MLE. To our experience, it is often useful to look at all three selected dimensions from AIC, BIC and LRT; and BIC selection is the most favorable method.

alpha = 0.05; % significance level for LRT
mindim = 0; % lower bound for the envelope dimensions (e.g. rank of beta, if known)
[dxs,dys] = tritests(X,Y,mindim,alpha);   %<-separate estimate dx and dy
disp(['Response envelope dimension is selected as: ', ...
    num2str(dys(1)),'(AIC), ',num2str(dys(2)),'(BIC), ',num2str(dys(3)),'(LRT at 0.05 level) ',]);
disp(['Predictor envelope dimension is selected as: ', ...
    num2str(dxs(1)),'(AIC), ',num2str(dxs(2)),'(BIC), ',num2str(dxs(3)),'(LRT at 0.05 level) ',]);

Response envelope dimension is selected as: 5(AIC), 3(BIC), 2(LRT at 0.05 level) 
Predictor envelope dimension is selected as: 4(AIC), 3(BIC), 2(LRT at 0.05 level) 

K-fold cross-validation prediction error

Given the envelope dimensions, compute the k-fold cross-validation prediction errors for various estimators.

k = 10;  % number of folds
nsim = 1;   % number of repeated k-fold cross-validation
dx = 3;
dy = 3;
FG=0; % use 1D estimators
mse = kFoldCV_senv(X,Y,dx,dy,FG,nsim,k);
disp(['OLS 10-fold cross-validation MSE is: ',num2str(mse.ols)])
disp(['Response envelope 10-fold cross-validation MSE is: ',num2str(mse.yenv)])
disp(['Predictor envelope 10-fold cross-validation MSE is: ',num2str(mse.xenv)])
disp(['Simultaneous envelope 10-fold cross-validation MSE is: ',num2str(mse.senv)])

OLS 10-fold cross-validation MSE is: 124.5324
Response envelope 10-fold cross-validation MSE is: 123.0393
Predictor envelope 10-fold cross-validation MSE is: 122.9031
Simultaneous envelope 10-fold cross-validation MSE is: 122.6969

Published with MATLAB® R2015a