clear

% 1.Loading data

% Specify the filename of the Excel file to be read.
filename = 'crabsTitle.xlsx'
[a, name, aCell] = xlsread(filename);

% 'a' - numeric data, 'name'- the text data (like headers), 'aCell' will contain both text and numeric data in cell format

dataTable = a;
Y = dataTable(:, 1);         % Y - groupping variable
widths = dataTable(:, 4);    % widths
sattelts = dataTable(:, 5);  % sattelts
weights = dataTable(:, 6);   % 'weights'
catwidths = dataTable(:,7);  % 'catwidths'

% 2.Data presentation

% Y - groupping variable.
X = Y;
% Combine matrix 'Z' with multiple attributes for each entry in the dataset
Z = [widths sattelts weights catwidths];
XZ=[X,Z]; % it's useful for preprocessing

% gplotmatrix: a matrix of scatter plots and histograms for data (including heterogeneous data)
% Create an object  table for given data vectors:

tblmy = table(Y, widths, sattelts, weights, catwidths, ...
       'VariableNames', {'Y', 'Width', 'Sattelts', 'Weight', 'Catwidth'});
gplotmatrix(Z, [], Y, [], '^o*x+', [], 'on', '', ...
   {'Width', 'Sattelts', 'Weight', 'Catwidth'}, []);
title(['1.Width ', '2.Sattelts ', '3.Weight ','4.Catwidth']);

% 3.Estimate the normality of the data distribution
% a) by normplot

figure
nameT={'Width', 'Sattelts', 'Weight', 'Catwidth'}
for i = 1:4
    subplot(2,2,i);
    normplot(Z(:,i));title(nameT{i}) % sgtitle - header for total subplot
end

% For normally distributed data points should form an approximate straight line.  Empirical distribution function ecdf is constructed in accordance with empirical data, inverse of theoretical normal distribution function is determined analytically, based on an effective sample parameters and their superposition close to a straight line, since the empirical distribution is close to the theoretical one)
% b) by Statistical test (Kolmogorov-Smirnov) - is used for Standard Normal distribution,  control Width, but before let's do the standardization:

widthsS=(widths-mean(widths))/std(widths);
[h,p,KSstat,critValue] = kstest(widthsS)
figure
histfit(widthsS)

% 4.Preprocess the matrix form

XZinitial = [X, Z]; 

% Initialize an empty array 'NaNInd' to store indices of rows with NaNs.
NaNInd = [];
[n, k] = size(XZ);

% Loop through each column of 'XZ'.
for col = 1:k
    
        [i] = find(isnan(XZ(:,col)));  % indices of rows with NaNs 
    if ~isempty(i)
        % 'union' use because we can find the same row for other column 
        NaNInd = union(NaNInd, i); 
    end 
end

% Remove the rows with NaNs from 'XZ'.
XZ(NaNInd,:)=[];

% Standardization:

[p, q] = size(XZ);
m = mean(XZ); % vector of means 
s = std(XZ); % the standard deviation of each column
XZs = (XZ - repmat(m, p, 1)) ./ repmat(s,p,1);

% Two ways to get correlation coefficient matrix: (You cat compare using witout sign ;)

RhoXZs=XZs'*XZs/(p-1);    
[rho,P]=corr(XZs);  

% It can be seen that these five variables are related to each other, and the correlation between variables 'Width', 'Weight' and 'Catwidth' is the most significant.
% 5.Linear regression model fitting (matrix form); Weights - response:

% Rearrange the elements of the mean and standard deviation vector 
m = [m(1), m(2:5)];
s = [s(1), s(2:5)];

% The 4-th column of 'XZs' is the dependent variable in the regress.
X = XZs(:, 4);
XZs(:, 4) = [];
% Z now contains only the independent variables:
Z = XZs;
[n, k] = size(Z);

% Regression model is obtained with directly matrix approach, taking into account only variables with standard Normal distribution as response - weights and idependent variable - x2 (widths) - according to normplot ~ Z(:,2):

ZEd=Z(:,2); disp(' b2 (compare with  fitlm  0.83922) is equal:' );
bEd=Z(:,2)\X % XEd=ZEd*bEd; - our system; \ - such a division

XEd=ZEd*bEd;  % XEd predicted by regression; XEd(i) weight of i-th crab 
err=(X-XEd); figure; histfit(err);
disp('You can see outliers on histogramm!')
SqError=sqrt(var(X-XEd)) % standard deviation 
disp('SqError compare with result in fitlm');

VX=var(X);  disp('corrected R2 istead 0.787 (the same) is equal:'); 
R2matrixApproach=1-var(X-XEd)/VX 

%% 6. Fitting regression - fitlm
% Fit a linear regression model:  fitlm   is used  with Z as the independent variables and X as respnse (dependent variable) here.

disp('1-Intersept, x1-Y, x2-widths, x3 - sattelts, x4 - catwidths')
lrm = fitlm(Z,X,'linear')

%  7.Function regress; display model coefficient and confidence intervals:
% Perform linear regression using regress. This function estimates the linear regression coefficients ('b'),their confidence intervals (CI) ('bINT'),the residuals ('R'), the CI for residual  ('RINT') and other statistics ('STATS') from the data. The confidence level is set to 95%

[b, bINT, R, RINT, STATS] = regress(X, Z, 0.05);

% STATS=[ R2 stat, F-stat, F p-value,  error variance].


% 'errorbar' plots the estimated coefficients ('b') with their CI
% The points and error bars are colored red.
figure;
errorbar(1:4, b, bINT(:, 1), bINT(:, 2), 'color', 'red');

% R-error of regression:

figure;
errorbar(1:length(RINT(:, 1)), R, RINT(:, 1), RINT(:, 2), 'color', 'red');
err=R;
histfit(err, 11);
% Turn on the grid for the current figure.
grid on;

% The error distribution is close to the standard normal distribution, and the absolute value of most errors is less than 0.5
% If RINT(i,:) does not contain zero, then the i-th residual is larger than would be expected, at the 5% significance level. It means that the I-th observation is an outlier.
% 7.Remove outliers and compare original data with new data.

Xnew = X;
Znew = Z;

% Remove outliers:
  indicesRows = [];
  

% As before:

for row = 1:length(RINT(:,1))
    maxvalue = R(row) + RINT(row, 2);
    minvalue = R(row) + RINT(row, 1); % RINT(row, 1)-negative
    if maxvalue < 0 || minvalue > 0
        indicesRows(end+1) = row;     % outlier
    end
end
Xnew(indicesRows) = [];
Znew(indicesRows, :) = [];

% Reconstruct a new standardized matrix:
XZsnew = [Znew(:,1), Xnew, Znew(:,2), Znew(:,3),Znew(:,4)];

% Comparation:

lrmnew = fitlm(Znew, Xnew, 'linear')

% Coeffitient of determination - is the best!

