function result = cvRpcr(x,y,kmax,rmsecv,h,k)
%CVRPCR calculates the robust RMSECV (root mean squared error of cross-validation) curve
% for RPCR or the robust RMSEP (root mean squared error of prediction) value in a fast way.
% The R-RMSECV curve can be used to make a selection of the optimal number of
% components to include in the regression model. The function is used in rpcr.m.
%
% Input arguments:
% x : the explanatory variables
% y : the response variables
% kmax : the maximal number of components to be considered in the model.
% rmsecv : Optional. If equal to 1 (default), the rmsecv is computed.
% Else, rmsecv = 0 and then the rss and R2 are computed.
% h : the quantile used in RPCR
% k : optional, if equal to zero (default), the RMSECV is calculated.
% If different from 0, the RMSEP value is calculated.
% Output:
% If rmsecv = 1
% result.rmsecv : the r-rmsecv values (only if the RMSECV is computed)
% result.rmsep : the rmsep value (only if RMSEP is computed)
% result.residu : the residuals for every k = 1,...,kmax
% result.outWeights : the fixed weights used in the robust version of the RMSE.
% result.R2 : the coefficient of determination for every k.
% result.rss : the RSS values for every k.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
% http://wis.kuleuven.be/stat/robust.html
%
% Written by Sanne Engelen
% Last Update: 08/10/2004, 03/07/2006
% Last Revision: 03/07/2006
% some initialisations
n = size(x,1);
p = size(x,2);
q = size(y,2);
teller_if_lus = 0;
cutoffWeights = sqrt(chi2inv(0.975,q));
if nargin < 4
alfa=0.75;
h=floor(2*floor((n+kmax+2)/2)-n+2*(n-floor((n+kmax+2)/2))*alfa);
k = 0;
rmsecv = 1;
elseif nargin == 4
alfa=0.75;
h=floor(2*floor((n+kmax+2)/2)-n+2*(n-floor((n+kmax+2)/2))*alfa);
k = 0;
elseif nargin == 5
k = 0;
end
if q == 1
FixedWeights = weightscvLtsregres(x,y,kmax,h,k);
else
FixedWeights = weightscvMcdregres(x,y,kmax,h,k,cutoffWeights);
end
kmax=FixedWeights.kmax;
if rmsecv
weights = FixedWeights.weightsk;
resrob = FixedWeights.resrob;
if k == 0
w_min = FixedWeights.w_min;
end
if size(resrob.Hsubsets.H0,2)==1
resrob.Hsubsets.H0=resrob.Hsubsets.H0';
end
Hsets = [FixedWeights.Hsets;resrob.Hsubsets.H0;resrob.Hsubsets.H1;resrob.Hsubsets.Hfreq];
for i = 1:n
disp(['observation ',num2str(i),' is left out'])
x_min_i = removal(x,i,0);
y_min_i = removal(y,i,0);
same.value = 0;
if isempty(find(resrob.Hsubsets.H0 == i))
if teller_if_lus >= 1
same.value = 1;
end
teller_if_lus = teller_if_lus + 1;
end
% constructing Hsets of right size: h - 1.
Hsets_min_i = RemoveObsHsets(Hsets,i);
if k == 0
res = removeObsRobpca(x,i,kmax,Hsets_min_i,same,1);
else
res = removeObsRobpca(x,i,k,Hsets_min_i,same);
end
if isempty(find(resrob.Hsubsets.H0 == i))
same.res = res;
end
Prob_min_i = res.Pk_min_i;
murob_min_i = res.muk_min_i;
Tk_min_i = [];
Lk_min_i = [];
coeffk_min_i = [];
if k == 0
j_ind = kmax;
else
j_ind = k;
end
Tkmax_min_i = (x_min_i - repmat(murob_min_i,n-1,1))*Prob_min_i(:,1:j_ind);
if q == 1
[resLts2,rawLts2] = ltsregres(Tkmax_min_i,y_min_i,'plots',0,'Hsets',Hsets_min_i,'h',h-1);
else
if k == 0
Mcdres = mcdcov([Tkmax_min_i,y_min_i],'h',h-1,'plots',0,'Hsets',Hsets_min_i);
% the full mcdregres is not needed, we only need the results of mcdcov.
resMcd2.Mu = Mcdres.center';
resMcd2.Sigma = Mcdres.cov;
end
end
if k == 0
j_start1 = 1;
j_end1 = kmax;
else
j_start1 = k;
j_end1 = k;
end
for j = j_start1:j_end1
if k == 0
if q == 1
[Bk,intk,weightslts] = extractlts(rawLts2,x_min_i, y_min_i, murob_min_i, Prob_min_i,kmax,j,h-1,n-1);
else
Tk_min_i = (x_min_i - repmat(murob_min_i,n-1,1))*Prob_min_i(:,1:j);
[Bk,intk,sigmayykmaxrew_k,sigmattkmaxrew_k] = extractmcdregres(resMcd2,Tk_min_i,y_min_i,kmax,n-1,q,j,h-1,cutoffWeights);
coeffk = [Bk;intk];
end
else
if q == 1
Bk = resLts2.slope;
intk = resLts2.int;
weightslts = resLts2.flag;
else
resMcdregres = mcdregres(Tkmax_min_i,y_min_i,'Hsets',Hsets_min_i,'plots',0,'h',h-1);
Bk = resMcdregres.slope;
intk = resMcdregres.int;
coeffk = [Bk;intk];
end
end
finalB = Prob_min_i(:,1:j)*Bk;
finalInt = intk - murob_min_i*finalB;
yhat_min_i = x(i,:)*finalB + finalInt;
residk_min_i(i,(j - 1)*q + 1:j*q) = y(i,:) - yhat_min_i;
% calculation of the resd:
if q > 1
if k == 0
rewE2=sigmayykmaxrew_k-coeffk(1:j,1:q)'*sigmattkmaxrew_k*coeffk(1:j,1:q);
cen=zeros(q,1);
resd(i,j)=sqrt(libra_mahalanobis(residk_min_i(i,(j - 1)*q + 1:j*q),cen','cov',rewE2))'; %robust distances of residuals
else
resd(i,j) = sqrt(libra_mahalanobis(residk_min_i(i,(j - 1)*q + 1:j*q),zeros(q,1),'cov',resMcdregres.cov));
end
else
if k == 0
scale=sqrt(sum(weightslts.*residk_min_i(i,(j - 1)*q + 1:j*q).^2)/(sum(weightslts)-1));
resd(i,j) = residk_min_i(i,(j-1)*q + 1:j*q)/scale;
else
resd(i,j) = residk_min_i(i,(j-1)*q + 1:j*q)/resLts2.scale;
end
end
end
end
if k == 0
for j = 1:kmax
resk = residk_min_i(:,(j-1)*q + 1:j*q);
if q == 1
rmsecv_min(j) = sqrt(1/sum(w_min)*w_min*(resk).^2);
else
rmsecv_min(j) = sqrt(1/sum(w_min)*w_min*(mean((resk').^2))');
end
end
result.rmsecv = rmsecv_min;
result.residu = residk_min_i;
else
weights = weights(:,k);
if q == 1
rmsep = sqrt(1/sum(weights)*weights'*(residk_min_i(:,k)).^2);
else
rmsep = sqrt(1/sum(weights)*weights'*(mean((residk_min_i(:,(k-1)*q + 1:k*q)').^2))');
end
result.rmsep = rmsep;
result.residu = residk_min_i(:,(k-1)*q + 1:k*q);
end
end
result.outWeights = FixedWeights;
result.rss = FixedWeights.rss;
result.R2 = FixedWeights.R2;
%-------------------------------------------------------------------------
function Hsets_min_i = RemoveObsHsets(Hsets,i)
% Removes the correct index from the $h$-subsets in Hsets to
% obtain (h - 1)-subsets.
% Every h-set is put as a row in Hsets.
% i is the index of the observation that is removed from the whole data.
for r = 1:size(Hsets,1)
if ~isempty(find(Hsets(r,:)== i))
Hsets_min_i(r,:) = removal(Hsets(r,:),0,find(Hsets(r,:) == i));
else
Hsets_min_i(r,:) = Hsets(r,1:(end-1));
end
for j = 1:length(Hsets_min_i(r,:))
if Hsets_min_i(r,j) > i
Hsets_min_i(r,j) = Hsets_min_i(r,j) - 1;
end
end
end
%------------------------------------------------------------------------
function out = weightscvLtsregres(x,y,kmax,h,k)
% Computes the weights needed in the R-RMSECV function.
%
% Input arguments:
% x : the independent variables
% y : the response variables
% kmax : the maximal number of components to be considered.
% h : the number of observations on which the calculations are based.
% k : if zero (default), the kmax approach is used (for RMSECV).
% Else robpca is calculated for a specific k (for RMSEP).
%
% Output arguments:
% if k = 0 (RMSECV):
% out.w_min : the weights obtained by taking the minimum over all k
% out.weightsk : the weights for every observation and every k (n x kmax).
% out.resrob : the results of robpca on x for kmax components.
% if k ~= 0 (RMSEP):
% out.weightsk : the weights for a specific k
% out.resrob : the results of robpca on x for k components.
% out.R2 : the weighted Rsquared for each value of k
% out.rss : the weighted rss for each value of k
n = size(x,1);
p = size(x,2);
q = size(y,2);
if nargin < 5
k = 0;
end
if k == 0
ResRobWhole = robpca(x,'plots',0,'kmax',kmax,'h',h,'k',kmax);
else
ResRobWhole = robpca(x,'plots',0,'kmax',kmax,'h',h,'k',k);
end
kmax = ResRobWhole.kmax;
Trob = ResRobWhole.T;
Mrob = ResRobWhole.M;
Prob = ResRobWhole.P;
HsetsH0 = ResRobWhole.Hsubsets.H0;
HsetsH1 = ResRobWhole.Hsubsets.H1;
HsetsHfreqrob = ResRobWhole.Hsubsets.Hfreq;
[ResLtsWhole,RawLtsWhole]=ltsregres(Trob,y,'plots',0,'h',h);
HsetsHfreq = ResLtsWhole.Hsubsets.Hfreq;
if k == 0
j_start = 1;
j_end = kmax;
else
j_start = k;
j_end = k;
end
for j = j_start:j_end
if k == 0
if j ~= kmax
[Bk,intk,weightslts,HsetsHopt(j,:)] = extractlts(RawLtsWhole,x,y,Mrob,Prob,kmax,j,h,n);
else
[Bk,intk,weightslts] = extractlts(RawLtsWhole,x,y,Mrob,Prob,kmax,j,h,n);
HsetsHopt(kmax,:) = ResLtsWhole.Hsubsets.Hopt';
end
else
Bk = ResLtsWhole.slope;
intk = ResLtsWhole.int;
weightslts = ResLtsWhole.flag;
HsetsHopt = ResLtsWhole.Hsubsets.Hopt';
end
coeffk=[Bk;intk];
finalB = Prob(:,1:j)*Bk;
finalInt = intk - Mrob*finalB;
fitted=x*finalB + repmat(finalInt,n,1);
residuals(:,j)=y-fitted;
scale=sqrt(sum(weightslts.*residuals(:,j).^2)/(sum(weightslts)-1));
s0=scale;
weightsk(:,j)=abs(residuals(:,j)/scale)<=sqrt(chi2inv(0.975,1));
end
% Determining fixed weights.
if k == 0
if kmax == 1
w_min = weightsk';
else
w_min = min(weightsk');
end
out.w_min = w_min;
yw = mean(y(w_min == 1,:));
y2=(y-repmat(yw,n,1)).^2;
R = residuals.^2;
D=sum(y2(w_min==1,:));
for j = 1:kmax
R1=R(w_min==1,j);
rss(j) = sum(sum(R1));
R2(j)=1-rss(j)/sum(D);
end
out.rss = 1/sum(w_min)*rss;
out.R2 = R2;
else
yw = mean(y(weightsk(:,k) == 1,:));
y2=(y-repmat(yw,n,1)).^2;
R = residuals.^2;
D=sum(y2(weightsk(:,k)==1,:));
R1=R(weightsk(:,k)==1,:);
rss(k) = sum(sum(R1));
R2(k)=1-rss(k)/sum(D);
out.rss = 1/sum(weightsk(:,k))*rss(k);
out.R2 = R2(k);
end
out.weightsk = weightsk;
out.resrob = ResRobWhole;
if size(HsetsH0,2)==1
HsetsH0=HsetsH0';
end
out.kmax=kmax;
out.Hsets = [HsetsHopt;HsetsHfreq;HsetsH0;HsetsH1;HsetsHfreqrob];
%-----------------------------------------------------------------------------
function [Bk,intk,weightslts,Hopt] = extractlts(rawltskmax, x_min_i, y_min_i, murob_min_i, Prob_min_i,kmax,k,h,n)
% This function extracts lts-regression coefficients for a certain k based on
% the slope and intercept of the regression with kmax coefficients.
%
% Input arguments:
% rawltskmax : the raw results of ltsregression with kmax components.
% x_min_i : the regressors without observation i.
% y_min_i : the dependent variable without observation i.
% murob_min_i : the center estimated using robpca on the data minus observation i.
% Prob_min_i : the loadings matrix estimated using robpca on the data minus observation i.
% k : the number of components
%
% Output arguments:
% Bk : the slope for k components.
% intk : the intercept for k components.
% weightslts : the weights needed for the reweighting.
% Hopt : the optimal h-subset. Not availabe for kmax.
j = k;
coeffkraw_min_i = rawltskmax.coefficients;
sloperaw = coeffkraw_min_i(1:j,:);
intraw = coeffkraw_min_i(end,:);
Tk_min_i = (x_min_i - repmat(murob_min_i,n,1))*Prob_min_i(:,1:j);
% perform csteps on the raw estimates of the parameters:
prevobj = 0;
if j ~= kmax
for noCsteps = 1:10
residu = y_min_i - Tk_min_i*sloperaw - repmat(intraw,n,1);
[sortresid,indsr] = sort(residu.^2);
obj = sum(sortresid(1:h));
[Q,R] = qr([Tk_min_i(indsr(1:h),:) ones(h,1)],0);
z = R\(Q'*y_min_i(indsr(1:h),1));
sloperaw = z(1:j,:);
intraw = z(j+1,:);
if noCsteps >= 20 | abs(obj - prevobj) < 10^(-4)
if noCsteps >= 20
disp('no convergence in Csteps')
end
break
end
prevobj = obj;
end
Hopt = indsr(1:h)';
end
% reweighting:
factor =rawconsfactorlts(h,n);
residu = y_min_i - Tk_min_i*sloperaw - repmat(intraw,n,1);
[sortresid,indsr] = sort(residu.^2);
sh0 = sqrt(1/(h)*sum(sortresid(1:h)));
s0 = sh0*factor;
m = 2*(n)/asvarscalekwad((h),(n));
quantile = tinv(0.9875,m);
weightslts = abs(residu/s0) <= quantile;
[Q,R] = qr([Tk_min_i(weightslts == 1,:) ones(sum(weightslts),1)],0);
z = R\(Q'*y_min_i(weightslts == 1));
Bk = z(1:j,:);
intk = z(j+1,:);
%----------------------------------------------------------------------------------------------------------------------
function out = weightscvMcdregres(x,y,kmax,h,k,cutoffWeights)
% computes the weights needed in the R-PRESS function.
%
% Input arguments:
% x : the independent variables
% y : the response variables
% kmax : the maximal number of components to be considered.
% h : the number of observations on which the calculations are based.
% k : if zero (default), the kmax approach is used then (for RMSECV).
% Else robpca is calculated for a specific k (for RMSEP).
%
% Output arguments:
% if k = 0 (RMSECV):
% out.w_min : the weights obtained by taking the minimum over all k
% out.weightsk : the weights for every observation and every k (n x kmax).
% out.resrob : the results of robpca on [X,Y] for kmax components.
% if k ~= 0 (RMSEP):
% out.weightsk : the weights for a specific k
% out.resrob : the results of robpca on [X,Y] for k components.
% out.R2 : the weighted Rsquared for each value of k
% out.rss : the weighted rss for each value of k
n = size(x,1);
p = size(x,2);
q = size(y,2);
if nargin < 5
k = 0;
end
if k == 0
ResRobWhole = robpca(x,'plots',0,'k',kmax,'kmax',kmax,'h',h);
else
ResRobWhole = robpca(x,'plots',0,'k',k,'kmax',kmax,'h',h);
end
kmax=ResRobWhole.kmax;
Trob = ResRobWhole.T;
Prob = ResRobWhole.P;
Mrob = ResRobWhole.M;
HsetsH0 = ResRobWhole.Hsubsets.H0;
HsetsH1 = ResRobWhole.Hsubsets.H1;
HsetsHfreqrob = ResRobWhole.Hsubsets.Hfreq;
ResMCDWhole = mcdregres(Trob,y,'h',h,'plots',0);
HsetsHfreq = ResMCDWhole.Hsubsets.Hfreq;
if k == 0
for j = 1:kmax
if j ~= kmax
if k == 0
Tk = (x - repmat(Mrob,n,1))*Prob(:,1:j);
[Bk,intk,sigmayykmaxrew_k,sigmattkmaxrew_k,HsetsHopt(j,:)] = extractmcdregres(ResMCDWhole,Tk,y,kmax,n,q,j,h,cutoffWeights);
else
Tk = (x - repmat(Mrob,n,1))*Prob(:,1:j);
[Bk,intk,sigmayykmaxrew_k,sigmattkmaxrew_k,HsetsHopt(j,:)] = extractmcdregres(ResMCDWhole,Tk,y,j,n,q,j,h,cutoffWeights);
end
else
[Bk,intk,sigmayykmaxrew_k,sigmattkmaxrew_k] = extractmcdregres(ResMCDWhole,Trob,y,kmax,n,q,j,h,cutoffWeights);
HsetsHopt(kmax,:) = ResMCDWhole.Hsubsets.Hopt;
end
coeffk = [Bk;intk];
finalB = Prob(:,1:j)*Bk;
finalInt = intk - Mrob*finalB;
yhat = x*finalB + repmat(finalInt,n,1);
residu(:,(j-1)*q+1:j*q) = y - yhat;
cen=zeros(q,1)';
cov=sigmayykmaxrew_k - coeffk(1:j,1:q)'*sigmattkmaxrew_k*coeffk(1:j,1:q);
[nn,pp]=size(residu(:,(j-1)*q+1:j*q));
resd = sqrt(libra_mahalanobis(residu(:,(j-1)*q+1:j*q),cen,'cov',cov))';
weightsk(:,j) = (abs(resd)<=cutoffWeights);
end
else
ResMCDWhole = mcdregres(Trob(:,1:k),y,'h',h,'plots',0);
if size(ResMCDWhole.flag,2)~=1
ResMCDWhole.flag = ResMCDWhole.flag';
end
weightsk(:,k) = ResMCDWhole.flag;
HsetsHfreq = ResMCDWhole.Hsubsets.Hfreq;
HsetsHopt = ResMCDWhole.Hsubsets.Hopt;
end
% Determining fixed weights.
if k == 0
if kmax == 1
w_min = weightsk';
else
w_min = min(weightsk');
end
out.w_min = w_min;
yw = mean(y(w_min == 1,:));
y2=(y-repmat(yw,n,1)).^2;
R = residu.^2;
D=sum(y2(w_min==1,:));
for j = 1:kmax
R1=R(w_min==1,(j-1)*q+1:j*q);
rss(j) = sum(sum(R1));
R2(j)=1-rss(j)/sum(D);
end
out.rss = 1/(q*sum(w_min))*rss;
out.R2 = R2;
else
yw = mean(y(weightsk(:,k) == 1,:));
y2=(y-repmat(yw,n,1)).^2;
R = ResMCDWhole.res.^2;
D=sum(y2(weightsk(:,k)==1,:));
R1=R(weightsk(:,k)==1,:);
rss = sum(sum(R1));
R2=1-rss/sum(D);
out.rss = 1/(q*sum(weightsk(:,k)))*rss;
out.R2 = R2;
end
out.weightsk = weightsk;
out.resrob = ResRobWhole;
if size(HsetsH0,2)==1
HsetsH0=HsetsH0';
end
out.kmax=kmax;
out.Hsets = [HsetsHopt;HsetsHfreq;HsetsH0;HsetsH1;HsetsHfreqrob];
%---------------------------------------------------------------
function rawconsfaclts=rawconsfactorlts(quan,n)
rawconsfaclts=(1/sqrt(1-((2*n)/(quan*(1/norminv((quan+n)/(2*n)))))*...
normpdf(1/(1/(norminv((quan+n)/(2*n)))))));
%--------------------------------------------------------------
function asvar=asvarscalekwad(quan,n)
alfa=quan/n;
alfa=1-alfa;
qalfa=chi2inv(1-alfa,1);
c2=gamcdf(qalfa/2,1/2+1);
c1=1/c2;
c3=3*gamcdf(qalfa/2,1/2+2);
asvar=qalfa*(1-alfa)-c2;
asvar=asvar^2;
asvar=(c3-2*qalfa*c2+(1-alfa)*(qalfa^2))-asvar;
asvar=c1^2*asvar;