Lebel2010Model

%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';
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';
anovaBootTable = '/Volumes/LANDLAB/projects/hbn/projectTrackProfiles/supportFiles/anovaBootTable.csv';

differenceTable = table('Size', [0, 2], 'VariableTypes', {'string', 'double'}, 'VariableNames', {'TractID', 'Difference'});
% Initialize a table to store the tract name and AIC value for each tract
aicTable = table('Size', [0, 2], 'VariableTypes', {'string', 'double'}, 'VariableNames', {'TractName', '2010_AIC'});

%Convert csv into a table.
Tshort = readtable(Tshort); 
Tlong = readtable(Tlong); 
DiagnosisT = readtable(Diagnosis); 
colorProfiles = readtable(colorProfiles);
rsqTableAdj = readtable(rsqTableAdj);
rsqTableOrd = readtable(rsqTableOrd);
inflecTable = readtable(inflecTable);
anovaBootTable = readtable(anovaBootTable); 
anovaBootTable.Hemisphere = categorical(anovaBootTable.Hemisphere);
anovaBootTable.Tract = categorical(anovaBootTable.Tract);

%For ANOVA table
lastRow = 1; 

%============== Generate Plots ==============
% 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];

% 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 = tbl.Sex; 
    yVar = tbl.yVar; 
    
    tbl = table(Age, yVar);

    %======================================================================
    % 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 poisson model function: A*age*exp(-B*age) + C
    modelFun = @(b, Age) b(1) * Age .* exp(-b(2) * Age) + b(3);

    % Set initial guesses for parameters A, B, and C
    initialGuess = [1, 0.01, 0.1];  % Adjust these based on your 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');
    
    %======= 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(['Poisson (Lebel 2010) 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 = 50;
    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];

    %===========================================================================
    % 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];
    
    %===========================================================================
    %{
    %Curvature Analysis
    
    % 1. Calculate curvature
    curvature = abs(d2y_dx2) ./ (1 + dy_dx.^2).^(3/2);

    % 2. Curvature threshold
    curvature_threshold = 0.02; % Adjust this based on your data’s scale

    % 3. Define the margin as a percentage of the total length (e.g., 2% of each endpoint)
    n_points = length(x_interp);
    exclude_idx = round(0.02 * n_points); % Exclude 5% of points at each endpoint

    % 4. Define the analysis region, excluding the endpoints
    analysis_region = exclude_idx + 1 : n_points - exclude_idx;

    % 5. Find points within the analysis region where curvature exceeds the threshold
    high_curvature_idx = analysis_region(curvature(analysis_region) > curvature_threshold);

    % 6. Get the corresponding x and y values for high curvature points
    high_curvature_x = x_interp(high_curvature_idx);
    high_curvature_y = y_smooth(high_curvature_idx);

    % 7. Plot the results
    hold on;
    scatter(high_curvature_x, high_curvature_y, 300, 'm', 'filled'); % High curvature points

    % 8. Plot curvature as a secondary axis for visualization
    yyaxis right
    plot(x_interp, curvature, 'm-', 'LineWidth', 1.2); % Curvature plot
    ylabel('Curvature (\kappa)');
    
    % 9. Display detected maximum and minimum curvature points
    disp('High Curvature Points (x, y):');
    disp([high_curvature_x, high_curvature_y]);
    
    hold off; 
    %} 
    %===========================================================================
    %Add adjusted R-squared & AIC values to table.
    %rsqTableAdj.MultNonLin(t) = mdlci.Rsquared.Adjusted;
    %rsqTableOrd.MultNonLin(t) = mdlci.Rsquared.Ordinary;
    %aicTable.MultNonLin(t)= mdlci.ModelCriterion.AIC;

    %legend
    %lgd = legend(f, {"Inflection Point","Fastest Rate of Change"});
    %lgd.FontName = 'Arial';
    %lgd.FontSize = 18;
    %legend box off;
    %pbaspect([1 1 1]);
    
    %=========================================================================================
     %{
    % 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); 

    %============== Bootstrapping for ANOVA ==============
    hold on 

    figure(5*(t + length(tractIDs)))

    N = 100; 

    for q = 1:10000

        %Select random sample of N

        %define Age, Sex, and measurement variable
        x = tbl.Age; 
        y = tbl.yVar;
        z = tbl.Sex; 

        %randomly select x, y, and z values of size N
        msize = size(x);
        idx = randperm(msize(1), N);
        x = x(idx(1,:));
        y = y(idx(1,:));
        z = z(idx(1,:));

        tbl2 = table(x, y, z);

        %generating the model
        Q2 = 'y ~ x^2 + z'; 
        mdl2 = fitlm(tbl2, Q2);

        %plotting the model
        h2 = plotAdjustedResponse(mdl2, 'x', 'visible', 'off');  
        set(h2, 'DefaultFigureVisible', 'off');

        %get data for plotting the confidence intervals and add CI to plot.
        j2 = array2table(cat(2, h2(1).XData', h2(1).YData')); j2.Properties.VariableNames =  {'x', 'y'};
        mdlci2 = fitlm(j2, 'y~x^2');

        %Find the inflection points
        %ydt = detrend(y,1); 
        % Detrend 'y' To Facilitate Analysis
      

        %plot(x, y, '-*'); hold on;

        %======= Calculating Inflection & Fastest Rate of Change =======
        x = unique(mdlci2.Variables.x); 
        y = predict(mdlci2, x);

        % Detrend 'y' To Facilitate Analysis
        y = detrend(y,1);                                     

        %Store x and y values in ratetbl
        ratetbl = table(x, y);
        %remove duplicate points
        ratetbl = unique(ratetbl, 'rows');
     
        %Assign x and y variables to new x and y vectors without duplicate
        %points
        y = ratetbl.y; 
        x = ratetbl.x;
    
        % Calculate Numerical Derivative
        dydx = gradient(y) ./ gradient(x);
        dyydx  = gradient(dydx) ./ gradient(x);

        %Save unordered derivatives in ratetbl                              
        ratetbl.dydx = dydx;
        ratetbl.dyydx = dyydx;

        %Sort derivatives 
        ratetbl = sortrows(ratetbl, 'dydx');                                      
        [~, ind] = unique(ratetbl(:,1), 'first');
        ratetbl = ratetbl(ind, :);
     
        [~, ind] = unique(ratetbl(:,"dydx"), 'first');
        ratetbl = ratetbl(ind, :);
        dydx = ratetbl.dydx; 
    
        %Interpolation Index Lower Limit
        [maxdydx, idxmax] = max(dydx);
        [mindydx, idxmin] = min(dydx);
        idxrng = idxmin: idxmax; 
    
        %inflection_idx = find(diff(sign(diff(y)))) + 1;
    
        %Find Inflection Point X-Value
        inflptx = interp1(dydx(idxrng), ratetbl.x(idxrng), 0, 'linear');           

        %Find Inflection Point Y-Value
        x = unique(mdlci.Variables.x);
        y = predict(mdlci, x);
        inflpty = interp1(x, y, inflptx, 'linear');   

        %plot(inflptx, inflpty, '*'); hold off;

        anovaBootTable.SampleNum(lastRow) = q; 
        anovaBootTable.TractIDs(lastRow) = tractIDs(t); 
        anovaBootTable.MultInflecX(lastRow) = inflptx; 
        anovaBootTable.MultInflecY(lastRow) = inflpty; 
        
        tractIDsTemp = string(tractIDs{t});
        if  tractIDsTemp(1) == "l"
            anovaBootTable.Hemisphere(lastRow) = 'left';
            anovaBootTable.Tract(lastRow) = string(extractAfter(tractIDs(t), 4));
        end

        if  tractIDsTemp(1) == "r"
            anovaBootTable.Hemisphere(lastRow) = 'right';
            anovaBootTable.Tract(lastRow) = string(extractAfter(tractIDs(t), 5));
        end

        lastRow = lastRow + 1; 

    end
    
    hold off 
    %} 
    
         
end

%============== 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);