function result = mlr(x,y,varargin)
%MLR is the classical least squares estimator for multivariate multiple
% linear regression. It can handle both one or several ('multiple')
% predictor variables and one or several ('multivariate') response variables.
% If the regression model contains an intercept, the x-matrix may not contain
% a column with ones.
% If there is only one response variable, the function ols.m is called.
%
% See for example:
% R.A. Johnson and D.W. Wichern,
% "Applied multivariate statistical analysis (Fifth Edition)",
% Prentice Hall, chapter 7.
%
% Required input arguments:
% x : Data matrix of the explanatory variables
% (n observations in rows, p variables in columns)
% y : Data matrix of the response variables
% (n observations in rows, q variables in columns)
%
% Optional input arguments:
% intercept : logical flag: if 1, a model with constant term will be
% fitted; if 0, no constant term will be included. (default: 1)
% plots : If equal to one, a menu is shown which allows to draw several plots,
% such as residual plots and a regression outlier map. (default)
% If 'plots' is equal to zero, all plots are suppressed.
% See also makeplot.m
%
% I/O: result=mlr(x,y,'plots',0);
% The user should only give the input arguments that have to change their default value.
%
% The output is a structure containing:
%
% result.slope : Slope estimate
% result.int : Intercept estimate
% result.fitted : Fitted values
% result.res : Residuals
% result.stdres : Standardized residuals
% result.cov : Estimated variance-covariance matrix of the residuals
% result.rsquared : R-squared value
% result.md : Score distances (Mahalanobis distances in x-space)
% result.resd : Residual distances (when there are several response variables).
% If univariate regression is performed, it contains the standardized residuals.
% result.cutoff : Cutoff values for the score and residual distances
% result.flag : The observations whose residual distance is larger than result.cutoff.resd
% receive a flag equal to zero. The other observations receive a flag 1.
% result.class : 'MLR' (when q > 1) or 'LS' (when q = 1)
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
% http://wis.kuleuven.be/stat/robust.html
%
% Written by Nele Smets on 13/12/2003
% Last update: 05/04/2004
%
q=size(y,2);
p=size(x,2);
geg=[x,y];
[n,m]=size(geg);
intercept=ones(n,1);
default=struct('plots',1,'intercept',1);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
counter=1;
if nargin>2
%
% placing inputfields in array of strings
%
for j=1:nargin-2
if rem(j,2)~=0
chklist{i}=varargin{j};
i=i+1;
end
end
% Checking which default parameters have to be changed
% and keep them in the structure 'options'.
while counter<=IN
index=strmatch(list(counter,:),chklist,'exact');
if ~isempty(index) % in case of similarity
for j=1:nargin-2 % searching the index of the accompanying field
if rem(j,2)~=0 % fieldnames are placed on odd index
if strcmp(chklist{index},varargin{j})
I=j;
end
end
end
options=setfield(options,chklist{index},varargin{I+1});
index=[];
end
counter=counter+1;
end
end
%%%%%%%%%Main part%%%%%%%
if q==1
result=libra_ols(x,y,'intercept', options.intercept,'plots',options.plots);
else
cmean=mean(geg);
ccovar=cov(geg);
meanx=cmean(1:p)';
meany=cmean((p+1):m)';
cy=mcenter(y);
covarx=ccovar(1:p,1:p);
covary=ccovar((p+1):m,(p+1):m);
covarxy=ccovar(1:p,(p+1):m);
covaryx=covarxy';
res.betas=[inv(covarx)*covarxy; (meany-(covaryx*inv(covarx)*meanx))'];
res.fitted=[x intercept]*res.betas;
res.residuals=y-res.fitted;
res.cov=covary-res.betas(1:p,1:q)'*covarx*res.betas(1:p,1:q);
res.stdresid= res.residuals./repmat(diag(res.cov)',n,1);
res.rsquared= 1-(det(res.residuals'*res.residuals)/det(cy'*cy));
res.class='MLR';
% x-distances (md) needed in diagnostic regression plot
if (-log(det(ccovar))/m) > 50
res.md='singularity';
res.resd='singularity';
else
res.md=sqrt(libra_mahalanobis(x,cmean(1:p),'cov',covarx))';
cen=zeros(q,1)';
res.resd=sqrt(libra_mahalanobis(res.residuals,cen,'cov',res.cov))';
end
%cutoff values
res.cutoff.md=sqrt(chi2inv(0.975,p));
res.cutoff.resd=sqrt(chi2inv(0.975,q));
res.flag=(abs(res.resd)<=res.cutoff.resd);
result=struct('slope',{res.betas(1:p,:)},'int',{res.betas(p+1,:)}, ...
'fitted',{res.fitted},'res',{res.residuals},'stdres',{res.stdresid},...
'cov',{res.cov},'rsquared',{res.rsquared},'md',{res.md},'resd',{res.resd},...
'cutoff',{res.cutoff},'flag',{res.flag},'class',{res.class});
if ~intercept
result=setfield(result, 'int', 0)
end
try
if options.plots==1
makeplot(result)
end
catch %output must be given even if plots are interrupted
%> delete(gcf) to get rid of the menu
end
end