Lebel2008Model

%Redirecting folder to correct path. clear.
clear all; clc;

%Define measure
measure = 'fa';

%Insert local path of Tshort.csv and Tlong.csv file
Tshort = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/Tshort.csv';
Tlong = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/Tlong.csv';
Diagnosis = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/Diagnosis.csv';
ASDdiagnosis = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/ASDSubjects.csv';
Motordiagnosis = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/MotorSubjects.csv';
colorProfiles = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/colorProfiles.csv';
rsqTableAdj = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/rsqTableAdj.csv';
rsqTableOrd = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/rsqTableOrd.csv';
inflecTable = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/inflecTable.csv';

differenceTable = table('Size', [0, 2], 'VariableTypes', {'string', 'double'}, 'VariableNames', {'TractID', 'Difference'});
anovaBootTable = table('Size', [0, 6], 'VariableTypes', {'double', 'string', 'string', 'string', 'double', 'double'}, ...
    'VariableNames', ...
    {'SampleNum', 'TractIDs', 'Hemisphere', 'Tract', 'AsymptoteY', 'AsymptoteX'});
% Initialize a table to store the tract name and AIC value for each tract
aicTable = table('Size', [0, 2], 'VariableTypes', {'string', 'double'}, 'VariableNames', {'TractName', '2008_AIC'});
outputTable = table('Size', [0, 6], 'VariableTypes', {'string', 'double', 'double', 'double'...
    , 'double', 'double'}, 'VariableNames', {'TractName', 'yRange', ...
    'minValue', 'maxValue', 'yat90', 'x_intercept_yat90'});

%Convert csv into a table.
Tshort = readtable(Tshort); 
Tlong = readtable(Tlong); 
DiagnosisT = readtable(Diagnosis); 
ASDT = readtable(ASDdiagnosis);
MotorDiagnosisT = readtable(Motordiagnosis); 
colorProfiles = readtable(colorProfiles);
rsqTableAdj = readtable(rsqTableAdj);
rsqTableOrd = readtable(rsqTableOrd);
inflecTable = readtable(inflecTable);
anovaBootTable.Hemisphere = categorical(anovaBootTable.Hemisphere);
anovaBootTable.Tract = categorical(anovaBootTable.Tract);
asymptoteTable = table('Size', [0, 3], 'VariableTypes', {'string', 'double', 'double'}, 'VariableNames', {'tractName', 'asymptote', 'xIntersect'});

summaryTable = table('Size', [0, 5], ...
    'VariableTypes', {'string', 'double', 'double', 'double', 'double'}, ...
    'VariableNames', {'TractName', 'N', 'AdjustedR2', 'FStatistic', 'timeConstant'});

%For ANOVA table
lastRow = 1; 

% Filter Tshort and Tlong to include only rows within the specified age range
ageMin = 0; % Minimum age
ageMax = 25; % Maximum age
Tshort = Tshort(Tshort.Age >= ageMin & Tshort.Age <= ageMax, :);
Tlong = Tlong(Tlong.Age >= ageMin & Tlong.Age <= ageMax, :);

%============== Generate Plots ==============
%Filter for subjects with ASD 
subjectsWithASD = ASDT.Identifiers(strcmp(ASDT.ASD, 'Yes') ...
    & startsWith(ASDT.Identifiers, 'N'));

%Filter for subjects with Motor Disorders
subjectsWithMotor = MotorDiagnosisT.Identifiers(strcmp(MotorDiagnosisT.Motor, 'Yes') ...
    & startsWith(MotorDiagnosisT.Identifiers, 'N'));

% Filter for subjects with "No Diagnosis" and IDs that start with "N"
subjectsWithNoDiagnosisN = DiagnosisT.Identifiers(strcmp(DiagnosisT.Diagnosis_ClinicianConsensus_DX_01, 'No Diagnosis Given') ...
    & startsWith(DiagnosisT.Identifiers, 'N'));

% Get all subjects whose IDs start with "P" in Tshort and Tlong
subjectsWithP = Tshort.subjectID(startsWith(Tshort.subjectID, 'P'));

% Combine both groups: subjects with "No Diagnosis" and IDs starting with "N", and all "P" subjects
validSubjects = [subjectsWithNoDiagnosisN; subjectsWithP]; %No Diagonosis 
%validSubjects = [subjectsWithASD]; %ASD
%validSubjects = [subjectsWithMotor]; %MotorDisorder

% Filter Tshort and Tlong to include only the valid subjects
Tshort = Tshort(ismember(Tshort.subjectID, validSubjects), :);
Tlong = Tlong(ismember(Tlong.subjectID, validSubjects), :);

%Generate column of tracts of interest ids
mask = ismember(Tlong.structureID, colorProfiles{:, 1});
tractIDs = Tlong(mask, :);
tractIDs = unique(tractIDs.structureID);
rsqSimpleLin = table(tractIDs); 

%Close all previous plots
close all

for t = 1:length(tractIDs)

    f = figure(t);

    %f.Position = [startingx startingy width height];
    f.Position = [1000 1000 800 700];

    hold on 

    %Defining model variables
    Age = Tshort.Age; 
    
    %define sex as a categorical variable.
    Sex = categorical(Tshort.Sex); 
    yVar = Tshort.(char(tractIDs(t)));
   
    %Create table
    tbl = table(Age, Sex, yVar); 
    
    %Remove any rows with missing data (NaN or empty character) from tbl 
    tbl(any(ismissing(tbl), 2), :) = [];

    %Assign Age, Sex, and yVar to new tbl columns with removed row that had missing values
    Age = tbl.Age; 
    Sex = categorical(tbl.Sex); 
    yVar = tbl.yVar; 
    
    tbl = table(Age, yVar);
    
    %======================================================================
    % Checking the effect of sex
    
    tblSex = table(Age, Sex, yVar); 
    
    % Fit a linear model to include Age and Sex as predictors
    lm_sex = fitlm(tblSex, 'yVar ~ Age + Sex');

    % Display the model summary
    disp(lm_sex);

    % Extract beta coefficient for Sex
    beta_sex = lm_sex.Coefficients.Estimate(strcmp(lm_sex.Coefficients.Row, 'Sex_1'));
    
    % Adjust FA values by removing the effect of Sex
    yVar_adjusted = tbl.yVar - (double(tblSex.Sex) * beta_sex);

    % Add the adjusted FA column to the table
    tbl.yVar = yVar_adjusted;

    % Verify the adjusted FA values
    %disp(head(Tshort)); % Display the first few rows to confirm

    %======================================================================
    % Outliers Identification 
    
    %Replace all outliers with zero
    %Removes outliers from yVar that is more than 3 sd from the mean
    yVar = filloutliers(yVar, 0, "mean"); 

    %Delete rows with outliers (identified because all outliers have been
    %replaced with zero)
    tbl(any(ismissing(tbl), 2), :) = [];
    tbl(~yVar, :) = []; 

    %======================================================================
    % Generating the model
    
    % Define the new model function: C + A * exp(-Age / t)
    modelFun = @(b, Age) b(1) + b(2) * exp(-Age / b(3));

    % Set initial guesses for parameters A, B, and C
    initialGuess = [0.5, 1, 10];  % Adjust these based on data

	%Fit the nonlinear model to the data
    mdl = fitnlm(tbl, modelFun, initialGuess);
    
    %Get appropriate RGB color for tract by indexing into colorProfiles.csv
    idx = find(strcmp(colorProfiles.NameOfTrack, char(tractIDs(t))) == 1);
    markerColor = [str2double(colorProfiles.Red{idx})/255, ...
               str2double(colorProfiles.Green{idx})/255, ...
               str2double(colorProfiles.Blue{idx})/255];
 
   % Define a range of ages for plotting the fitted response
    ageRange = linspace(min(tbl.Age), max(tbl.Age), 100)';

    % Convert ageRange to a table with the same column name as used in mdl
    ageRangeTable = table(ageRange, 'VariableNames', {'Age'});

    % Get the fitted response and prediction intervals
    [fittedResponse, confidenceInterval] = predict(mdl, ageRangeTable, 'Prediction', 'curve');

    % Extract the lower and upper bounds of the confidence interval
    ciLower = confidenceInterval(:, 1);
    ciUpper = confidenceInterval(:, 2);

    % Plot the original data
    figure;
    hold on
    scatter(tbl.Age, tbl.yVar, 100, 'MarkerEdgeColor', 'white', 'MarkerFaceColor', markerColor, 'Marker', 'o');

    % Plot the fitted poisson model
    plot(ageRange, fittedResponse, 'LineWidth', 2, 'Color', markerColor, 'DisplayName', 'Fitted Curve');

    % Plot the confidence intervals
    fill([ageRange; flipud(ageRange)], [ciLower; flipud(ciUpper)], markerColor, 'FaceAlpha', 0.2, 'EdgeColor', 'none', 'DisplayName', 'Confidence Interval');

    % Add labels and title
    xlabel('Age');
    ylabel('yVar');
    
    %================= Adding to Summary Table ===================
    % Number of subjects (N)
    N = height(tbl); 

    % Adjusted R-squared
    adjustedR2 = mdl.Rsquared.Adjusted;

    % Extract values from the model
    SSR = mdl.SSR; % Sum of Squares for Regression
    SSE = mdl.SSE; % Sum of Squares for Error
    n = mdl.NumObservations; % Total number of observations
    p = numel(mdl.Coefficients.Estimate); % Number of estimated parameters

    % Compute F-statistic
    FStatistic = (SSR / p) / (SSE / (n - p));
    
    % Extract the time constant (b(3))
    timeConstant = mdl.Coefficients.Estimate(3);

    % Add these to the summary table
    newRow = {char(tractIDs(t)), N, adjustedR2, FStatistic, timeConstant};
    summaryTable = [summaryTable; newRow];
    
    %================= Calculating Min-Max Metric ===================
    % Extract the estimated parameters from mdl
    %C = mdl.Coefficients.Estimate(1); % Asymptote value
    %A = mdl.Coefficients.Estimate(2);
    %T = mdl.Coefficients.Estimate(3);

    % Calculate the predicted values (y)
    x = linspace(min(tbl.Age), max(tbl.Age), 1000); 
    x_table = table(x', 'VariableNames', {'Age'});
    y_pred = predict(mdl, x_table); % Get predicted values using the fitted model
    
    % Extract minimum value from mdl 
    [minValue, minIndex] = min(y_pred);
    
    % Extract minimum value from mdl 
    [maxValue, maxIndex] = max(y_pred);
    
    % Calculate yRange 
    yRange = maxValue - minValue; 
    
    % Calculate yat90
    yat90 = 0.90*(yRange) + minValue; 
    
    % Find the x-intercept of yat90
    % This is the x-value (or Age) where y_pred is closest to yat90
    [~, intersectIndex] = min(abs(y_pred - yat90));
    x_intercept_yat90 = x(intersectIndex);
    
    % Plot this intersection point on the graph
    plot(x_intercept_yat90, yat90, 'ro', 'MarkerSize', 25, 'MarkerFaceColor', 'r');
    
    % Add to output table
    newRow = {char(tractIDs(t)), yRange, minValue, maxValue, yat90, x_intercept_yat90};

    % Store values in outputTable
    outputTable = [outputTable; newRow];
    
    %======= Calculating Inflection & Fastest Rate of Change =======
    %{
    % Define a high-resolution age range for interpolation
    x = linspace(min(tbl.Age), max(tbl.Age), 1000); 

    % Calculate the predicted values (y)
    x_table = table(x', 'VariableNames', {'Age'});
    y = predict(mdl, x_table); % Get predicted values using the fitted model

    % 1. Interpolate and smooth the data for higher resolution
    x_interp = linspace(min(x), max(x), 10000); % High-resolution x values
    y_interp = interp1(x, y, x_interp, 'pchip'); % Smooth interpolation
    y_smooth = smoothdata(y_interp, 'gaussian', 50); % Stronger Gaussian smoothing

    % 2. Calculate first and second derivatives
    dy_dx = gradient(y_smooth) ./ gradient(x_interp); % First derivative
    d2y_dx2 = gradient(dy_dx) ./ gradient(x_interp); % Second derivative

    % 3. Detect zero-crossings in the second derivative
    sign_changes = find(diff(sign(d2y_dx2)) ~= 0); % Detect strict sign changes

    % 4. Initialize arrays for inflection points, ensuring unique detections
    inflection_points_x = [];
    inflection_points_y = [];

    for i = 1:length(sign_changes)
        inflection_x = x_interp(sign_changes(i));
        inflection_y = y_smooth(sign_changes(i));
        % Ignore points near the edges, e.g., within 5% of each end of the x range
        if inflection_x > x_interp(round(0.05 * length(x_interp))) && ...
                inflection_x < x_interp(round(0.95 * length(x_interp)))
    
        % Ensure unique inflection points based on proximity
        if isempty(inflection_points_x) || all(abs(inflection_x - inflection_points_x) > 1e-3)
            inflection_points_x = [inflection_points_x, inflection_x];
            inflection_points_y = [inflection_points_y, inflection_y];
            
            % Display each unique inflection point
            disp(['Inflection point: (x, y) = (' num2str(inflection_x) ', ' num2str(inflection_y) ')']);
        end
         end
    end

    % 5. If no inflection points are found, display a message
    if isempty(inflection_points_x)
        disp({'No inflection points found for', char(tractIDs(t))});
    end

    %Plot the inflection points
    scatter(inflection_points_x, inflection_points_y, 100, 'red', 'filled');
    %} 
    
    hold off

    %===========================================================================
    %Style Settings for the Plot

    %Style the plot
    pltLeg = legend('', '', '');
    set(pltLeg,'visible','off')

    %Add title and color to the model
    plotTitle = {char(tractIDs(t))};
    plotTitle = strjoin(['Exponential Fitting Equation (Lebel 2008) Model for', plotTitle]);
    title(plotTitle);
    xlabel('Age (years)');
    ylabel(measure);

    % Set up plot and measure-specific details.
    capsize = 0;
    marker = 'o';
    linewidth = 1.5;
    linestyle = 'none';
    markersize = 100;
    xtickvalues = [1 2 3 4];
    xlim_lo = min(xtickvalues)-0.5; xlim_hi = max(xtickvalues)+0.5;
    fontname = 'Arial';
    fontsize = 30;
    fontangle = 'italic';
    yticklength = 0;
    xticklength = 0.02;

    % xaxis
    xax = get(gca, 'xaxis');
    xax.TickDirection = 'out';
    xax.TickLength = [xticklength xticklength];
    set(gca, 'XLim', [3 22], 'XTick', [3 12.5 22]);
    xax.FontName = fontname;
    xax.FontSize = fontsize;

    % yaxis
    yax = get(gca,'yaxis');
    yax.TickDirection = 'out';
    yax.TickLength = [yticklength yticklength];
    set(gca, 'YLim', [0.2 0.6], 'YTick', [0.2 0.4 0.6]);
    yax.FontName = fontname;
    yax.FontSize = fontsize;
    yax.FontAngle = fontangle;

    %change figure background to white
    set(gcf, 'color', 'w')
    
    %===========================================================================
    % Extract the AIC value for the fitted model
    aicValue = mdl.ModelCriterion.AIC;

    % Store the tract name and AIC value in the aicTable
    newRow = {char(tractIDs(t)), aicValue};
    aicTable = [aicTable; newRow];
    %=========================================================================================
    % Define the custom path and filename
    customPath = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/finalPlots'; % Replace with your desired folder path
    fileName = [char(tractIDs(t)), '.png']; % Specify the file name
    fullFilePath = fullfile(customPath, fileName); % Combine path and file name
    
    % Save as a png image
    saveas(gcf, fullFilePath);
    
    %=========================================================================================
     %{
    % Calculate the difference between the first and last points of mdlci
    first_point = mdlci.Variables.y(1);
    last_point = mdlci.Variables.y(end);
    difference = last_point - first_point;
    
     tractIDx = string(tractIDs{t});
    % Append the new row to the table
    newRow = {tractIDx, difference};
    differenceTable = [differenceTable; newRow];

    hold off

    %add adjusted inflection point and fastest rate to table.
    inflecTable.MultInflecX(t) = inflptx; 
    inflecTable.MultInflecY(t) = inflpty; 
    inflecTable.MultFastRateX(t) = fr.x(1); 
    inflecTable.MultFastRateY(t) = fr.y(1); 

    %}     
end

%============== Parametric Bootstrapping (Omnibus Anova) ==============
observedXint = outputTable.x_intercept_yat90;
numTracts = length(tractIDs);

% Pool all data points across tracts into a single dataset

% Define the tracts of interest
tractsOfInterest = {'leftMDLFang', 'leftMDLFspl', 'leftTPC', 'leftpArc', 'leftVOF'...
    'rightMDLFang', 'rightMDLFspl', 'rightTPC', 'rightpArc', 'rightVOF', 'Age'};

% Define the columns to exclude
excludeColumns = {'subjectID', 'Sex', 'Handedness'}; 

% Get a list of all columns in the table
allColumns = Tshort.Properties.VariableNames;

% Select only the columns for the tracts of interest, excluding unwanted columns
selectedColumns = setdiff(intersect(allColumns, tractsOfInterest), excludeColumns);

% Filter the table to include only the selected columns
pooledData = Tshort(:, selectedColumns);

% Stack the PVP and VOF columns into a single variable
T_stacked = stack(pooledData, pooledData.Properties.VariableNames(2:end), ...
                  'NewDataVariableName', 'FA_Value', ...
                  'IndexVariableName', 'Tract_Name');

% The new table will have columns: 'Age', 'Tract_Name', and 'FA_Value'
T_stacked = T_stacked(:, {'Age', 'FA_Value', 'Tract_Name'});

pooledData = table(T_stacked.Age, T_stacked.FA_Value);

%Remove any rows with missing data (NaN or empty character) from tbl 
pooledData(any(ismissing(pooledData), 2), :) = [];

% Set the size of each sample
N = 400; 

% Set the number of permutations
numPermutations = 10000;

% Define model function and initial guess for parameters
modelFun = @(b, Age) b(1) + b(2) * exp(-Age / b(3));
initialGuess = [0.5, 1, 10]; % Initial guess for model parameters

% Initialize vector to store permuted x_intersect values
permutedXIntersects = zeros(numPermutations, 1);

% Permutation testing loop
for i = 1:numPermutations
    % Randomly sample N points from the pooled data with replacement
    sampleIndices = randsample(height(pooledData), N, true);
    sampleData = pooledData(sampleIndices, :);
    sampleData.Properties.VariableNames = {'Age', 'FA_Value'}; 

    % Fit the nonlinear model to the sampled data
    tbl = table(sampleData.Age, sampleData.FA_Value, 'VariableNames', {'Age', 'FA_Value'});
    
    if height(tbl) > 1  % Ensure there's enough data to fit the model
        % Calculate the predicted values (y)
        x = linspace(min(tbl.Age), max(tbl.Age), 1000); 
        x_table = table(x', 'VariableNames', {'Age'});
        y_pred = predict(mdl, x_table); % Get predicted values using the fitted model
    
        % Extract minimum value from mdl 
        [minValue, minIndex] = min(y_pred);
    
        % Extract minimum value from mdl 
        [maxValue, maxIndex] = max(y_pred);
    
        % Calculate yRange 
        yRange = maxValue - minValue; 
    
        % Calculate yat90
        yat90 = 0.90*(yRange) + minValue; 
    
        % Find the x-intercept of yat90
        % This is the x-value (or Age) where y_pred is closest to yat90
        [~, intersectIndex] = min(abs(y_pred - yat90));
        x_intercept_yat90 = x(intersectIndex);
        
        % Store the x_intersect for this permutation
        permutedXIntersects(i) = x_intercept_yat90;
    else
        % If not enough data, set to NaN
        permutedXIntersects(i) = NaN;
    end
end

% Calculate p-values for each tract's observed x_intersect
pValues = zeros(numTracts, 1);
significanceLevel = 0.05; % Define the significance level
significantLabels = strings(numTracts, 1); % Initialize significance labels

for t = 1:numTracts
    % Remove NaNs from the permuted x_intersects
    validPermutedXIntersects = permutedXIntersects(~isnan(permutedXIntersects));
    
    % Calculate p-value by comparing observed x_intersect to permuted distribution
    pValues(t) = mean(abs(validPermutedXIntersects - mean(validPermutedXIntersects)) >= abs(observedXint(t) - ...
        mean(validPermutedXIntersects)));
    
    % Determine significance based on the p-value
    if pValues(t) < significanceLevel
        significantLabels(t) = "Yes";
    else
        significantLabels(t) = "No";
    end
end

% Calculate the mean and standard deviation of the permuted x-intersects
meanPermutedXIntersect = mean(validPermutedXIntersects);
stdPermutedXIntersect = std(validPermutedXIntersects);

% Calculate z-scores for each observed x-intercept
zScores = (observedXint - meanPermutedXIntersect) / stdPermutedXIntersect;

% Append z-scores to the output table
outputTable.zScore = zScores;

% Append p-values label to the output table
outputTable.pValue = pValues;

% Add a new column to outputTable that is pValue divided by 10
outputTable.pValue_div10 = outputTable.pValue / 10;

% Append significance label to the output table
outputTable.Significant = significantLabels;

% Filter outputTable to include only rows with tract names in tractsOfInterest
outputTable = outputTable(ismember(outputTable.TractName, tractsOfInterest), :);

% Plot the histogram of permuted x-intersects
figure;
histogram(validPermutedXIntersects, 30, 'FaceColor', [0.7, 0.7, 0.7], 'EdgeColor', 'k');
hold on;

% Loop over each tract to overlay observed x-intersect with label only if it is in tractsOfInterest
for t = 1:numTracts
    if ismember(tractIDs{t}, tractsOfInterest)
        % Plot a vertical line for the observed x-intersect
        xline(observedXint(t), 'LineWidth', 2, 'Color', 'r');
        
        % Add text label for the tract name near the vertical line
        text(observedXint(t) + 0.1, max(ylim) * 0.95, tractIDs{t}, ...
             'Color', 'r', 'FontSize', 10, 'HorizontalAlignment', 'left', ...
             'Rotation', 90);
    end
end
% Customize plot
xlabel('X-Intercept at 90% Y Level');
ylabel('Frequency');
title('Distribution of Permuted X-Intercepts with Observed Values');
legend('Permuted X-Intercepts', 'Observed X-Intercepts');
hold off;

% Rank order the z-scores in descending order (highest z-score gets rank 1)
[~, zScoreRank] = sort(outputTable.zScore, 'descend');
outputTable = sortrows(outputTable, 'zScore', 'descend');

% Display and save results
disp(outputTable);
writetable(outputTable, 'PermutationTest_XIntersectResults.csv');

% Calculate the range of the permutation distribution
minValue = min(validPermutedXIntersects);
maxValue = max(validPermutedXIntersects);
distributionRange = maxValue - minValue;

% Display the range
disp(['Range of the permutation distribution: ', num2str(distributionRange)]);
disp(['Minimum value: ', num2str(minValue)]);
disp(['Maximum value: ', num2str(maxValue)]);

% Calculate the mean and standard deviation of the permutation distribution
distributionMean = mean(validPermutedXIntersects);
distributionStd = std(validPermutedXIntersects);

% Display the results
disp(['Mean of the permutation distribution: ', num2str(distributionMean)]);
disp(['Standard deviation of the permutation distribution: ', num2str(distributionStd)]);

%================== Follow-up T-Tests (2 Tracts / Time) ===================

% Define unique tract names
tractNames = unique(outputTable.TractName);

% Initialize variables to store t-test results
numPairs = nchoosek(length(tractNames), 2); % Number of pairs
tractPairNames = strings(numPairs, 1);
pValues = zeros(numPairs, 1);
testStatistics = zeros(numPairs, 1);
significanceLabels = strings(numPairs, 1); 

% Pair counter
pairIdx = 1;

% Loop through all pairs of tracts
for i = 1:length(tractNames) - 1
    for j = i + 1:length(tractNames)
        % Get tract names
        tract1 = tractNames(i);
        tract2 = tractNames(j);
        
        % Extract FA values for each tract
        values1 = T_stacked.FA_Value(T_stacked.Tract_Name == tract1);
        values2 = T_stacked.FA_Value(T_stacked.Tract_Name == tract2);
        
        % Perform a two-sample t-test
        [~, p, ~, stats] = ttest2(values1, values2);
        
        % Determine significance and store label
        if p < significanceLevel
            significanceLabels(pairIdx) = "Yes";
        else
            significanceLabels(pairIdx) = "No";
        end
        
        % Store the results
        tractPairNames(pairIdx) = tract1 + " vs " + tract2;
        pValues(pairIdx) = p;
        testStatistics(pairIdx) = stats.tstat;
        
        % Increment pair counter
        pairIdx = pairIdx + 1;
    end
end

% Add a new column to tTestResults that is pValue divided by the number of
% permutations
pValue_div = pValues / length(pValues);

% Create a table for the t-test results with significance labels
tTestResults = table(tractPairNames, pValues, pValue_div, testStatistics, significanceLabels, ...
                     'VariableNames', {'Tract_Pair', 'pValue', 'pValue_div', 'TStatistic', 'Significant'});

% Sort tTestResults by TStatistic in ascending order
tTestResults = sortrows(tTestResults, 'TStatistic', 'ascend');

% Display the t-test results
disp(tTestResults);

%========================= Summary Table =========================

% Display the summary table
disp(summaryTable);

%================== Hemispheric Differences T-Tests ===================
% Initialize table to store results
tractComparison = table('Size', [0, 5], ...
    'VariableTypes', {'string', 'double', 'double', 'double', 'double'}, ...
    'VariableNames', {'TractName', 'MeanDifference', 'StandardDeviation', 'tStatistic', 'pValue'});

% Extract unique tract suffixes (e.g., MDLFang, MDLFspl, TPC)
uniqueSuffixes = unique(erase(tractIDs, ["left", "right"]));

% Loop through each suffix to compute paired t-test
for t = 1:length(uniqueSuffixes)
    % Get suffix (e.g., "MDLFang")
    tractSuffix = uniqueSuffixes{t};

    % Generate left and right tract names
    leftTractName = strcat("left", tractSuffix);
    rightTractName = strcat("right", tractSuffix);

    % Extract diffusion metrics for left and right tracts
    leftXIntercept = Tshort.(char(leftTractName));
    rightXIntercept = Tshort.(char(rightTractName));

    % Ensure paired comparison by removing missing values
    validIndices = ~isnan(leftXIntercept) & ~isnan(rightXIntercept);
    leftXIntercept = leftXIntercept(validIndices);
    rightXIntercept = rightXIntercept(validIndices);

    % Check if data is valid for comparison
    if isempty(leftXIntercept) || isempty(rightXIntercept)
        disp(['Skipping comparison for ', tractSuffix, ' due to missing data.']);
        continue;
    end

    % Calculate differences
    differences = leftXIntercept - rightXIntercept;

    % Check for zero variance (identical values)
    if all(differences == differences(1))
        disp(['All differences for ', tractSuffix, ' are identical. No variability.']);
        meanDifference = differences(1); % Mean of identical values
        sdDifference = 0;                % Standard deviation is zero
        tStatistic = NaN;                % t-statistic is undefined
        pValue = NaN;                    % p-value is undefined
    else
        % Compute mean and standard deviation
        meanDifference = mean(differences);
        sdDifference = std(differences);

        % Compute sample size
        n = length(differences);

        % Compute standard error
        seDifference = sdDifference / sqrt(n);

        % Compute t-statistic
        tStatistic = meanDifference / seDifference;

        % Compute p-value (two-tailed test)
        pValue = 2 * (1 - tcdf(abs(tStatistic), n - 1));
    end

    % Add results to the table
    newRow = {tractSuffix, meanDifference, sdDifference, tStatistic, pValue};
    tractComparison = [tractComparison; newRow];
end

% Display results
disp('Per-Tract Paired t-Test Results:');
disp(tractComparison);

%============== Export rsqTable as a csv ==============
%local path to save table: 
mainpath = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/';

table_path_format_rsqTAdj = fullfile(mainpath, 'rsqTableAdj.csv');
table_path_format_rsqTOrd = fullfile(mainpath, 'rsqTableOrd.csv');
table_path_format_inflecT = fullfile(mainpath, 'inflecTable.csv');
table_path_format_aicTable = fullfile(mainpath, 'aicTable.csv');
table_path_format_anovaBootTable = fullfile(mainpath, 'anovaBootTable.csv');

%finally, save tables
writetable(rsqTableAdj, table_path_format_rsqTAdj);
writetable(rsqTableOrd, table_path_format_rsqTOrd);
writetable(inflecTable, table_path_format_inflecT);
writetable(aicTable, table_path_format_aicTable);
writetable(anovaBootTable, table_path_format_anovaBootTable);