%% Behaviour of the routines fitlm and stepwiselm
% Batch de démonstration pour le cours de plans d'expériences. Le batch a pour but d'illustrer les méthodes reliées aux fonctions fitlm et  stepwiselm. Cela se fait au travers de l'analyse de mesures publiées par Hydro Quebec pour un brevet d'acier inoxydable austénitique au cobalt ultra résistant à la cavitation érosive.
% Batch demonstrating the behaviour of the least square fit routines for the course of design of experiments. This is done analysing some public data on austenic steel. Austenitic cobalt stainless steel is ultra resistant to erosive cavitation
% © Jean-Marie Fürbringer, EPFL, 2023
%
% Table of Contents
% Load of the data from an external file
% Analysis of the hardness (dureté)
%     All possible factors
%     Use of ANOVA
%     Use of coefCI()
%     Use of coefTest()
%     Use of plot(mdl)
%     Use of plotAdded(mdl, coef)
%     Use of plotDiagnostics
% Only the significant factors of the first round
%     Plot of the effects with their confidence intervals
%     List of confidence intervals
%     Plot of residues vs  model
%     Plot of the residues by order of experiments
%         NormalPlot of residuals
%         Cook distances
%         Cook's distance vs residuals
% Without the oulier
%         Confidence intervals
% Step wise regression
% _______________________________________________________________________
%
%% Load of the data from an external file
% The data is loaded from an Excel file  called dataDOE2019master.xlsx and put in a table called Steel_tbl

Steel_tbl = readtable('dataDOE2019master.xlsx',...
    'Sheet','acier');

%
% Here is the data:

disp(Steel_tbl)

%
%% Analysis of the hardness (dureté)
%
%% All possible factors
%
% Linear regression of the hardness of the steel in function of the content in cupper (C), cobalt (Co), Manganese (Mn), Scilicium (Si), Chromium (Cr), Molibdene (Mo)  and Nitrogen (N). The model is described with the Wilkinson notation 'Durete~C+Mn+Si+Cr+Co+Mo+N' equivalent to
% Durete = a_o + a_1\ C + a_2\ Mn +a_3\ Si + a_4\ Cr + a_5\ Co +a_6\ Mo + a_7\ N + \epsilon
%

mdl_hardness_1 = fitlm(Steel_tbl,'Durete~C+Mn+Si+Cr+Co+Mo+N');
disp(mdl_hardness_1)

% ___________________________________________________________________
% ANALYSIS
% The linear coefficients for Mn (p=24%) Si (p=11%), Cr (p=53%), Co (p=81%) are clearly over the usual limit of 5%.
% The coefficient for Mo is close of the limit with p\approx7\%.
% Equivalently the corresponding values of tStat are lower than the criterion t(5%/2,N-P=20)=2.08.
% The F statistic of 8.87 is higher than the criterion F(5%,1,N-P=20)=4.35 indicating that the model is globaly acceptable.
% It makes sense to test a simpler model without the questionable terms.
% _______________________________________________________________________
%%
%
%% Use of ANOVA
%
% First try with the routine with default parameters:

anova(mdl_hardness_1)

% The results are equivalent to the fitlm() output, exept that the test is now based on the Fisher distribution.
% The linear coefficients for Mn (p=24%) Si (p=11%), Cr (p=53%), Co (p=81%) and Mo (p=7%) are over the limit of 5%
% The p-values are slightly different, but the figure is strickly the same
% In the next step, the parameter of the routine allow to precise that all the components and the sum of squares of type I (sequential) are expected:

anova(mdl_hardness_1,'component',1)

% In this case, we observe strong differences in the sums of squares and the p-values. If the VIF are computed, it becomes evident that the data is not orthogonal with several VIF coefficients that are higher than 10. This is a situation that requires SS of type II due to the non orthoganality of the regressors.

% building of the model matrix
X=x2fx(Steel_tbl{:,1:7});

% compute the dispersion matrix
D=inv(X'*X);

% compute the correlation matrix
C=corrcov(D);

% compute the VIF
VIF=diag(inv(C));

% plot the VIF
bar(VIF(2:8))
xticklabels({'C','Mn','Si','Cr','Co','Mo','N'})
set(gca,"YScale","log")
grid on
title(' Variance inflation factors',"FontSize",16)
hold on
plot([0.25,7.75],[10,10],'r--',"LineWidth",3)
hold off

% The correlation plot unfortunately does not help to see precisely why the estimation of the coefficient for the molibden is sensitive to the type of sum of squares. Nevertheless it is evident that the data is not sufficient to make a safe estimation of the linear effects of the different elements.

corrplot(Steel_tbl{:,1:7})

% The anova computed with the parameter 'summary' confirms that the model is globaly ok with a p-value of 56 ppm

anova(mdl_hardness_1,'summary')

%
%% Use of coefCI()
% Let us now compute the confidence intervals of the coefficients with the routine coefCI()

CI_1=coefCI(mdl_hardness_1);
disp(array2table(CI_1))

% Let us produce a graphical representation of this data:

coef=mdl_hardness_1.Coefficients.Estimate;
bar(1:8,coef)
title("Estimated coefficients and CI","FontSize",16)
xticklabels({'a_o','C','Mn','Si','Cr','Co','Mo','N'})
hold on
H=errorbar(1:8,coef,CI_1(:,1),CI_1(:,2));
H.Color=[0,0,0];
H.LineStyle="none";
hold off

%
% The bad quality of the estimate is evident in this bar chart
% _____________________________________________________________________
%
%% Use of coefTest()
% This routine gives the possibility to test the values of the coefficients. We must be carefull between the difference of the ANOVA results giving the probability to have a null model and the test of the values of the coefficient.
% The routine is not straight forward... so be careful to what is really computed.
% In the first case, the routine compute the probability corresponding to the hypothesis of the ANOVA  H_o = "All coefficients are zero" against H_1= "not all the coefficients are zero".
% Use H to specify the coefficients to include in the F-test. So in this case thera are all used as H is the unity matrix:

H=eye(8);
[p,F,r]=coefTest(mdl_hardness_1,H)

% We see the H_1 wins as p< 0.05. The routine can also returns the F-test statistic F and the numerator degrees of freedom for the test , r.
% In this second case, the routine compute the probability corresponding to the hypothesis of the ANOVA  H_o = "All the coefficients have the value computed by the LSF algorithm" against H_1= "not all the coefficients have the estimated value".

H=eye(8);
[p,F,r]=coefTest(mdl_hardness_1,H,coef)

% We see that the probability is now 1, H_o wins. This result is not of interest, because obvious.
% Now let's test some specific values. We have forced the non-significant coefficients to zero and put some arbitrary values for the other ones, H_o = ( a_o=12,\ a_1=12,\ a_2=a_3=a_4=a_5=0,\ a_6=1.5,\ a_7=80 )

H=eye(8);
C=[11.8,12,0,0,0,0,1.5,78];
[p,F,r]=coefTest(mdl_hardness_1,H,C)

% The result shows that this hypothesis H_o is rejected. To have an acceptable model with only 4 factors it would be necessary to compute again the coefficients with less coefficients.
% Alternatively ,we can test the values of the significant coefficients with H_o = ( a_o=12, a_1=12, a_6=1.5, a_7=80 )

H=[1 0 0 0 0 0 0 0
   0 1 0 0 0 0 0 0
   0 0 0 0 0 0 1 0
   0 0 0 0 0 0 0 1];

C=[12,12,1.5,80];
[p,F,r]=coefTest(mdl_hardness_1,H,C)

%
%% Use of plot(mdl)
% This function is a version of "model vs data" plot. It compares the data with the model, without the constant a_o.  Is is very basic. A sophistication is possible with the next routine plotadded(mdl) which will let us to select the regressors.

plot(mdl_hardness_1)

%% Use of plotAdded(mdl, coef)
% Without specification this routine is equivalent to plot(mdl). But its value lies in the possibility to include only some of the regressors.
% In the example below, only the second regressor, C, and the seventh, N, have been considered and is is possible to observe that the figure is quite equivalent to the whole model and showing that there is an important part of the data that is not taken in consideration by the model (the points outside the confidence hyperbole).
%

plotAdded(mdl_hardness_1,[2,8])
plotAdded(mdl_hardness_1,8)
plotAdded(mdl_hardness_1,2)

%
%% Use of plotDiagnostics
%
% The routine plotDiagnostics(mdl,type) gives access to several graphics putting in evidence points whose leverage on the data is important and that would merit a more intensive scrutiny to be sure that the conclusions are not driven by erroneous data. Cook’s distance is the scaled change in fitted values, which is useful for identifying outliers in the X values (observations for predictor variables). Cook’s distance shows the influence of each observation on the fitted response values. An observation with Cook’s distance larger than three times the mean Cook’s distance might be an outlier. We see in the first graphic that two points (2, 28) have effectively a Cook's distance which is higher than 2 (the Cook's distance  is 0.65), in fact 8 and 9

M=mean(mdl_hardness_1.Diagnostics.CooksDistance)
plotDiagnostics(mdl_hardness_1,"cookd")
plotDiagnostics(mdl_hardness_1,"contour")
plotDiagnostics(mdl_hardness_1,"leverage")

%
%
%% Only the significant factors of the first round
%
% Inference of a linear model of the hardness based on factors that show p-values close or smaller than 5%

mdl_durete_2=fitlm(Steel_tbl,'Durete~C+Mn+Mo+N');
disp(mdl_durete_2)

%
% _______________________________________________________________________
% ANALYSIS
% The remaining four linear coefficients now have p-values below the usual 5% threshold.
% The pValue of the constant has also improved significantly.
% The values of tStat are all greater than the criterion t(5%/2,N-P=23)=2.07.
% A model based on the concentration of C, Mn, Mo and N seems acceptable: a p-value for the model of about 4\ 10^{-6} .
% At this stage it is interesting to look at the diagnostics.
% _______________________________________________________________________
%% Plot of the effects with their confidence intervals

figure
plotEffects(mdl_durete_2)

%
%% List of confidence intervals

CI=table(coefCI(mdl_durete_2),...
    'VariableNames',{'CI'},...
    'RowNames',{'const','C','Mn','Mo','N'});
disp(CI)

%
% _______________________________________________________________________
% ANALYSE
% The four linear coefficients have relatively wide confidence intervals with respect to their value.
% The prediction capacity of the model is therefore very limited.
% At this stage, it is interesting to look at the diagnostics.
% _______________________________________________________________________
%
%% Plot of residues vs  model
%

figure
plotResiduals(mdl_durete_2,'fitted')

%
%% Plot of the residues by order of experiments
%

figure
plotResiduals(mdl_durete_2,'caseorder')

%
%% NormalPlot of residuals

figure
plotResiduals(mdl_durete_2,'probability')

%
%% Cook distances

figure
plotDiagnostics(mdl_durete_2,'cookd')

%
%% Cook's distance vs residuals

figure
plotDiagnostics(mdl_durete_2,'contour')

%
% _______________________________________________________________________
% ANALYSIS
% No specific pattern in residual analysis.
% Case 28 (the last) at a distance from Cook large enough to suggest that it could be an outlier.
% On the other hand, the contour plot does not seem to confirm this hypothesis.
% Let's see what happens when we remove this point.
% _______________________________________________________________________
%% Without the oulier

[~,larg] = max(mdl_durete_2.Diagnostics.CooksDistance);

mdl_durete_3=fitlm(Steel_tbl,'Durete~C+Mn+Mo+N','Exclude',larg);
disp(mdl_durete_3)
figure
plotEffects(mdl_durete_3)

%
%% Confidence intervals

CI3=table(coefCI(mdl_durete_3),...
    'VariableNames',{'CI'},...
    'RowNames',{'const','C','Mn','Mo','N'});
disp(CI3)

%
% _______________________________________________________________________
% ANALYSIS
% It does not significantly improve the model.
% _______________________________________________________________________
%% Step wise regression

mdl_durete_4 = stepwiselm(Steel_tbl,'Durete~C+Mn+Si+Cr+Co+Mo+N',...
    'PredictorVars',{'C','Mn','Si','Cr','Co','Mo','N'})

%
% _______________________________________________________________________
% ANALYSIS
% The method quickly leads to a better result by integrating interaction term
% _______________________________________________________________________