Problems with implementing cross-validation for OLS, PLSR and PCR

Question

I am new to regression methods. I am creating Multiple Linear Regression, Partial Least Squares Regression and Principal Component Regression models for my dataset, and I am a bit confused with the results.

I have 60 predictors (or features) in my data, and 71 observations. Based on what I haver read (and understood), PLSR should give the best results. PCR should give results that are close to PLSR, and MLR is the most simplistic and which should give the worst results. Is my conclusion correct?

Because in my dataset, MLR and PLSR work very very similarly whereas PCR gives a much better prediction result. I have done a comparison of the methods, watching how the RMSE of the cross-validated methods changes with the number of features introduced for the MLR and PLSR methods, and with the number of PCs used for the PCR method. While MLR and PLSR suffer from overfitting very early, PCR keeps improving its prediction accuracy while the number of PCs is increased.

What do you think about this? Can this have a logical explanation?

---

This is what I am doing to compute PCR fit and prediction errors (Matlab code):

PCALoadings_all=PCALoadings;
PCAScores_all=PCAScores;

for ik=1:length(PCAVar) %add PCs one by one

PCALoadings=PCALoadings_all(:,1:ik);
PCAScores= PCAScores_all(:,1:ik);

%Function 'regress' needs a constant column so add it
PCAScores = [PCAScores, ones(length(PCAScores(:,1)),1)];
[betaPCR,~,~,~,statsPCR]  = regress(labels-mean(labels), PCAScores); %Create regression model


betaPCR=betaPCR(1:end-1,:); %remove constant effect


size_PCR= find(betaPCR~=0, 1, 'last'); %find 0-value coefficients

%To make the PCR results easier to interpret in terms of the original
%spectral data, transform to regression coefficients for the original, uncentered variables.
betaPCR = PCALoadings*betaPCR;

%Create model
betaPCR_fit = [mean(labels) - mean(norm_feature_matrix)*betaPCR; betaPCR];
yfitPCR = [ones(d1,1) norm_feature_matrix]*betaPCR_fit; %compute fitted labels

%Cross-validate
predictMat=[];
for test_idx=1:length(labels)

    testix = test_idx; %One value to test
    trainix = setdiff(1:length(labels),testix); %All the other values to train

    %Feature matrix and labels to train
    X=norm_feature_matrix(trainix,:);
    y=labels(trainix);

    %Feature matrix and labels to test
    X_test=norm_feature_matrix(testix,:);
    y_test=labels(testix);


    %compute beta with training data
    betaPCR_CV = [mean(y) - mean(X)*betaPCR; betaPCR];
    %prediction with test data
    yCV_PCR = [1 X_test]*betaPCR_CV;


    predictrow = [y_test, yCV_PCR]; % actual vs predicted
    predictMat = [predictMat; predictrow];

end

RMSE_fit_PCR(ik) = sqrt(mean((yfitPCR - labels).^2)); %Compute RMSE for regression fit 
RMSE_prediction_PCR(ik) = sqrt(mean((predictMat(:,1) - predictMat(:,2)).^2)); %Compute RMSE for prediction
%plot regressions

end

Do you think I'm doing something wrong? Did I misunderstood something?

Answer 1

Ok, I think I've got it...

These are my new curves (which I think are much more logical):

In fact, there were several mistakes in my code:

1.-In the PLSR part, I was changing the number of features instead of the number of PCs in each iteration. Nonsense.

2.-My Cross-Validation process in all of them was wrong. I was implicitly using the values to test when creating the Cross-Validation model, as @amoeba remarked.

The PCR code is now like this:

  for ik=1:length(PCAVar) %add PCs one by one

PCALoadings=PCALoadings_all(:,1:ik);
PCAScores= PCAScores_all(:,1:ik);

%Function 'regress' needs a constant column so add it
PCAScores = [PCAScores, ones(length(PCAScores(:,1)),1)];
[betaPCR,~,~,~,statsPCR]  = regress(labels-mean(labels), PCAScores); %Create regression model


betaPCR=betaPCR(1:end-1,:); %remove constant effect


size_PCR= find(betaPCR~=0, 1, 'last'); %find 0-value coefficients

%To make the PCR results easier to interpret in terms of the original
%spectral data, transform to regression coefficients for the original, uncentered variables.
betaPCR = PCALoadings*betaPCR;

%Create model
betaPCR_fit = [mean(labels) - mean(norm_feature_matrix)*betaPCR; betaPCR];
yfitPCR = [ones(d1,1) norm_feature_matrix]*betaPCR_fit; %compute fitted labels
RMSE_fit_PCR(ik) = sqrt(mean((yfitPCR - labels).^2)); %Compute RMSE for regression fit 
end

%Cross-validate
predictMat=[];
error_pcnum=[];
for test_idx=1:length(labels)

testix = test_idx; %One value to test
trainix = setdiff(1:length(labels),testix); %All the other values to train

%Feature matrix and labels to train
X=norm_feature_matrix(trainix,:);
y=labels(trainix);

%Feature matrix and labels to test
X_test=norm_feature_matrix(testix,:);
y_test=labels(testix);

[ PCALoadings_all,PCAScores_all,PCAVar,~] = compute_PCs( X );
predictrow_pcnum=[];
predictMat_pcnum=[];
for pcnum = 1:length(PCAVar)
    PCALoadings=PCALoadings_all(:,1:pcnum);
    PCAScores= PCAScores_all(:,1:pcnum);
    [betaPCR,~,~,~,~]  = regress(y-mean(y), PCAScores); %Create regression model
    betaPCR = PCALoadings*betaPCR;

    %compute beta with training data
    betaPCR_CV = [mean(y) - mean(X)*betaPCR; betaPCR];
    %prediction with test data
    yCV_PCR = [1 X_test]*betaPCR_CV;

    predictrow_pcnum = [y_test, yCV_PCR]; % actual vs predicted
    predictMat_pcnum = [predictMat_pcnum; predictrow_pcnum];
end
error_pcnum =[error_pcnum, (predictMat_pcnum(:,1) - predictMat_pcnum(:,2))];
predictMat = [predictMat; predictMat_pcnum'];
end
error_mean_bypcs=mean(error_pcnum.^2,2);
RMSE_prediction_PCR = sqrt(error_mean_bypcs); %Compute RMSE for prediction
[~,min_RMSE_nfeat]=min(RMSE_prediction_PCR); %find the minimum RMSE

%RMSE_prediction_PCR(ik) = sqrt(mean((predictMat(:,1) - predictMat(:,2)).^2)); %Compute RMSE for prediction
%plot regressions

predictMat=predictMat(:,min_RMSE_nfeat);
predictMat = reshape(predictMat,[2,length(predictMat)/2])';
with_vs_without_CV(labels, yfitPCR, predictMat, 3, sup_id);
legend(strcat(num2str(min_RMSE_nfeat) ,' PCs used in the PCR model'),'Location','southeast')

Thank you to everyone for your collaboration!

(And if you think I am still wrong, please, let me know...)