%% =========================================================================
%  PaperB_MasterScript.m
%
%  Callaway, Dunne & Williams
%  "The Choice of Asymmetry Index Affects Clinical Interpretation of Hip
%   Adduction Strength in Elite Soccer Players: A Season-Long Comparison
%   of Seven Formulas"
%
%  REQUIREMENTS
%  ------------
%  Place this script in the same folder as:
%      Data_SeasonLong_Anonymous.xlsx   (anonymised dataset, sheet: Raw)
%
%  The script uses fileparts(mfilename('fullpath')) to locate the data
%  file, so no path editing is required.
%
%  MATLAB R2021b or later recommended (uses 'omitnan' flag).
%  Toolboxes: Statistics and Machine Learning Toolbox.
%
%  OUTPUTS
%  -------
%  All figures are displayed and can be saved manually (File > Save As).
%  Summary statistics are printed to the Command Window.
%  Threshold sensitivity data exported to:  PaperB_ThresholdSensitivity.xlsx
%  Stabilisation estimates exported to:     PaperB_Stabilisation.xlsx
%  95% CI half-widths exported to:          PaperB_CI_HalfWidths.xlsx
%
%  SECTION INDEX
%  -------------
%  1.  Load data
%  2.  Compute seven asymmetry formulas
%  3.  Descriptive statistics (Table 1)
%  4.  Clinical flagging rates (10% and 15%)
%  5.  Inter-formula Pearson correlations (Figure 2)
%  6.  Bland-Altman analysis: |BSA| vs |BAI-2| (Figure 4)
%  7.  Threshold sensitivity analysis (Figure 5)
%  8.  Direction-of-asymmetry and weekly kappa (Results section)
%  9.  Rolling direction consistency (Figure 6)
%  10. Season-long spaghetti plots: BSA vs BAI-2 (Figure 7)
%  11. Weekly asymmetry by formula (Figure 1)
%  12. Per-player flagging heatmap (Figure 3)
%  13. Individual vs position-group representativeness (Figures 8 & 9)
%  14. Per-player cumulative-mean stabilisation (Figure 10)
%  15. 95% CI half-width analysis
%  16. Abduction and ratio supplementary results
%  17. Helper functions
%
%  =========================================================================

clc; close all; clear;

% -------------------------------------------------------------------------
%  1. LOAD DATA
% -------------------------------------------------------------------------
% Locate data file relative to this script — no hard-coded paths
scriptDir  = fileparts(mfilename('fullpath'));
dataFile   = fullfile(scriptDir, 'Data_SeasonLong_Anonymous.xlsx');

if ~isfile(dataFile)
    error(['Data file not found.\n' ...
           'Expected: %s\n' ...
           'Place Data_SeasonLong_Anonymous.xlsx in the same folder as this script.'], ...
           dataFile);
end

opts = spreadsheetImportOptions("NumVariables", 23);
opts.Sheet       = "Raw";
opts.DataRange   = "A2:W602";
opts.VariableNames = ["Player_ID","Date","Participant","Week","MD", ...
    "MD_Group","Mass","Height","Position","Player_Group","Group", ...
    "Add_NonDom","Add_Dom","Abd_NonDom","Abd_Dom", ...
    "DomVsNonDom_Add","DomVsNonDom_Abd", ...
    "Dom_AddAbd","NonDom_AddAbd", ...
    "Norm_Add_NonDom","Norm_Add_Dom", ...
    "Norm_Abd_NonDom","Norm_Abd_Dom"];
opts.VariableTypes = ["categorical","datetime","double","double","categorical", ...
    "double","double","double","categorical","categorical","double", ...
    "double","double","double","double","double","double","double","double", ...
    "double","double","double","double"];
opts = setvaropts(opts, ["Player_ID","MD","Position","Player_Group"], ...
    "EmptyFieldRule","auto");
opts = setvaropts(opts, "Date", "InputFormat", "");

data = readtable(dataFile, opts, "UseExcel", false);

% Keep only valid adduction observations (non-NaN in both limbs)
validIdx = ~isnan(data.Add_NonDom) & ~isnan(data.Add_Dom);
data     = data(validIdx, :);

fprintf('Data loaded: %d valid player-by-week observations, %d unique players.\n', ...
    height(data), numel(unique(data.Participant)));

% Position group codes: 1=GK, 2=DEF, 3=MID, 4=FWD
posLabels = {'GK','DEF','MID','FWD'};
nWeeks    = 35;

% -------------------------------------------------------------------------
%  2. COMPUTE SEVEN ASYMMETRY FORMULAS
%  All values in % units. Signed = positive when dominant limb is stronger.
% -------------------------------------------------------------------------
DL  = data.Add_Dom;
NDL = data.Add_NonDom;

% Group 1 — normalised to single limb (dominant or stronger)
LSI2  = (1 - NDL ./ DL) .* 100;                          % signed
BSA   = abs(DL - NDL) ./ max(DL, NDL) .* 100;            % unsigned
PctD  = BSA;                                              % identical to BSA for unilateral

% Group 2 — normalised to mean of both limbs
BAI2  = 2 .* (DL - NDL) ./ (DL + NDL) .* 100;           % signed

% Group 3 — normalised to sum (total output)
BAI1  = (DL - NDL) ./ (DL + NDL) .* 100;                 % signed
SI    = abs(DL - NDL) ./ (DL + NDL) .* 100;              % unsigned

% Group 4 — Symmetry Angle (arctangent transformation)
SA    = (45 - atand(NDL ./ DL)) ./ 90 .* 100;            % signed

% Collect into table for convenience
T = table(LSI2, BAI2, BSA, PctD, BAI1, SI, SA, ...
    'VariableNames', {'LSI2','BAI2','BSA','PctD','BAI1','SI','SA'});

formulaNames  = {'LSI-2','BAI-2','BSA','%Diff','BAI-1','SI','SA'};
formulaCodes  = {'LSI2','BAI2','BSA','PctD','BAI1','SI','SA'};
nFormulas     = numel(formulaNames);
isSigned      = logical([1 1 0 0 1 0 1]);  % 1 = signed, 0 = unsigned

% -------------------------------------------------------------------------
%  3. DESCRIPTIVE STATISTICS — TABLE 1
% -------------------------------------------------------------------------
fprintf('\n--- TABLE 1: Descriptive statistics (n = %d) ---\n', height(data));
fprintf('%-8s %8s %8s %8s %8s %8s %8s %8s %8s\n', ...
    'Formula','Mean','|Mean|','SD','Median','CoV','Min','Max','10%Flag');

for f = 1:nFormulas
    vals = T.(formulaCodes{f});
    absV = abs(vals);
    mn   = mean(vals,  'omitnan');
    absMn= mean(absV,  'omitnan');
    sd   = std(vals,   'omitnan');
    med  = median(vals,'omitnan');
    cov_ = std(absV,'omitnan') / mean(absV,'omitnan') * 100;
    mn_  = min(vals,  [],  'omitnan');
    mx   = max(vals,  [],  'omitnan');
    flag = mean(absV >= 10) * 100;
    fprintf('%-8s %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f%%\n', ...
        formulaNames{f}, mn, absMn, sd, med, cov_, mn_, mx, flag);
end

% Within-player CoV
fprintf('\nWithin-player CoV (mean across players, absolute values):\n');
participants = unique(data.Participant);
for f = 1:nFormulas
    covs = NaN(numel(participants),1);
    for p = 1:numel(participants)
        v = abs(T.(formulaCodes{f})(data.Participant == participants(p)));
        if sum(~isnan(v)) >= 3
            covs(p) = std(v,'omitnan') / mean(v,'omitnan') * 100;
        end
    end
    fprintf('  %-8s: %.1f%%\n', formulaNames{f}, mean(covs,'omitnan'));
end

% Dominant limb stronger in what proportion of observations?
domStronger = mean(DL > NDL) * 100;
fprintf('\nDominant limb stronger in %.1f%% of observations.\n', domStronger);

% -------------------------------------------------------------------------
%  4. CLINICAL FLAGGING RATES
% -------------------------------------------------------------------------
thresholds = [10, 15];
fprintf('\n--- Flagging rates by threshold ---\n');
fprintf('%-8s %10s %10s\n','Formula','10% flag','15% flag');
flagRate10 = NaN(1, nFormulas);
flagRate15 = NaN(1, nFormulas);
for f = 1:nFormulas
    absV  = abs(T.(formulaCodes{f}));
    fr10  = mean(absV >= 10) * 100;
    fr15  = mean(absV >= 15) * 100;
    flagRate10(f) = fr10;
    flagRate15(f) = fr15;
    fprintf('%-8s %9.1f%% %9.1f%%\n', formulaNames{f}, fr10, fr15);
end

% -------------------------------------------------------------------------
%  5. INTER-FORMULA PEARSON CORRELATIONS — FIGURE 2
% -------------------------------------------------------------------------
corrMat = zeros(nFormulas);
for i = 1:nFormulas
    for j = 1:nFormulas
        v1 = T.(formulaCodes{i});
        v2 = T.(formulaCodes{j});
        ok = ~isnan(v1) & ~isnan(v2);
        corrMat(i,j) = corr(v1(ok), v2(ok), 'Type','Pearson');
    end
end

figure('Name','Figure 2 - Correlation matrix','NumberTitle','off');
imagesc(corrMat, [-1 1]);
colormap(redblue_colormap());
colorbar;
set(gca,'XTick',1:nFormulas,'XTickLabel',formulaNames, ...
        'YTick',1:nFormulas,'YTickLabel',formulaNames,'FontSize',11);
title('Inter-formula Pearson correlations (n = 429)','FontSize',12);
for i = 1:nFormulas
    for j = 1:nFormulas
        text(j, i, sprintf('%.2f', corrMat(i,j)), ...
            'HorizontalAlignment','center','FontSize',9, ...
            'Color', ternary(abs(corrMat(i,j)) > 0.5, 'w', 'k'));
    end
end

% Print between-cluster correlations
fprintf('\n--- Between-cluster correlations (signed vs unsigned) ---\n');
signedIdx   = find(isSigned);
unsignedIdx = find(~isSigned);
for i = signedIdx
    for j = unsignedIdx
        fprintf('  %s vs %s: r = %.2f\n', formulaNames{i}, formulaNames{j}, corrMat(i,j));
    end
end

% -------------------------------------------------------------------------
%  6. BLAND-ALTMAN: |BSA| vs |BAI-2| — FIGURE 4
% -------------------------------------------------------------------------
bsa_abs  = abs(T.BSA);
bai2_abs = abs(T.BAI2);
ba_mean  = (bsa_abs + bai2_abs) / 2;
ba_diff  = bsa_abs - bai2_abs;

meanDiff = mean(ba_diff, 'omitnan');
sdDiff   = std(ba_diff,  'omitnan');
LoA_lo   = meanDiff - 1.96 * sdDiff;
LoA_hi   = meanDiff + 1.96 * sdDiff;

% Proportional bias regression
ok = ~isnan(ba_mean) & ~isnan(ba_diff);
mdl = fitlm(ba_mean(ok), ba_diff(ok));
slope = mdl.Coefficients.Estimate(2);
pSlope = mdl.Coefficients.pValue(2);
R2     = mdl.Rsquared.Ordinary;

fprintf('\n--- Bland-Altman: |BSA| vs |BAI-2| ---\n');
fprintf('  Mean difference:  %.2f%%\n', meanDiff);
fprintf('  SD difference:    %.2f%%\n', sdDiff);
fprintf('  LoA:              %.2f%% to %.2f%%\n', LoA_lo, LoA_hi);
fprintf('  Proportional bias: slope = %.2f, p = %.3f, R² = %.2f\n', slope, pSlope, R2);

figure('Name','Figure 4 - Bland-Altman','NumberTitle','off');
scatter(ba_mean, ba_diff, 20, [0.4 0.4 0.8], 'filled', 'MarkerFaceAlpha', 0.4);
hold on;
xRange = [min(ba_mean(ok)) max(ba_mean(ok))];
yFit   = predict(mdl, xRange');
plot(xRange, yFit,   'b-',  'LineWidth', 1.5, 'DisplayName', ...
    sprintf('Bias regression (slope = %.2f, R²=%.2f)', slope, R2));
yline(meanDiff, 'k-',  'LineWidth', 1.5, 'DisplayName', ...
    sprintf('Mean diff = %.2f%%', meanDiff));
yline(LoA_hi,   'r--', 'LineWidth', 1.2, 'DisplayName', ...
    sprintf('LoA = %.2f to %.2f%%', LoA_lo, LoA_hi));
yline(LoA_lo,   'r--', 'LineWidth', 1.2, 'HandleVisibility','off');
xlabel('Mean of |BSA| and |BAI-2| (%)','FontSize',11);
ylabel('Difference: |BSA| \minus |BAI-2| (%)','FontSize',11);
title('Bland-Altman: |BSA| vs |BAI-2|','FontSize',12);
legend('Location','northwest','FontSize',9);
grid on; box on;

% -------------------------------------------------------------------------
%  7. THRESHOLD SENSITIVITY ANALYSIS — FIGURE 5
% -------------------------------------------------------------------------
threshRange = 0:0.5:25;
flagCurves  = zeros(nFormulas, numel(threshRange));
for f = 1:nFormulas
    absV = abs(T.(formulaCodes{f}));
    for ti = 1:numel(threshRange)
        flagCurves(f, ti) = mean(absV >= threshRange(ti)) * 100;
    end
end

% Export
exportTbl = array2table(flagCurves', ...
    'VariableNames', strrep(formulaNames, '-','_'), ...
    'RowNames', arrayfun(@(x) sprintf('%.1f%%',x), threshRange,'UniformOutput',false));
writetable(exportTbl, fullfile(scriptDir,'PaperB_ThresholdSensitivity.xlsx'), ...
    'WriteRowNames',true);

lineStyles = {'-','-','--','--','-.','-.',':'};
colors6    = lines(nFormulas);
figure('Name','Figure 5 - Threshold sensitivity','NumberTitle','off');
hold on;
for f = 1:nFormulas
    plot(threshRange, flagCurves(f,:), lineStyles{f}, ...
        'Color', colors6(f,:), 'LineWidth', 1.8, 'DisplayName', formulaNames{f});
end
xline(10, 'k:', 'LineWidth', 1, 'HandleVisibility','off');
xline(15, 'k:', 'LineWidth', 1, 'HandleVisibility','off');
text(10.2, 95, '10%', 'FontSize', 9);
text(15.2, 95, '15%', 'FontSize', 9);
xlabel('Threshold (%)','FontSize',11);
ylabel('Flagging rate (%)','FontSize',11);
title('Sensitivity of flagging rate to threshold choice','FontSize',12);
legend('Location','northeast','FontSize',9);
xlim([0 25]); ylim([0 100]); grid on; box on;

% -------------------------------------------------------------------------
%  8. DIRECTION-OF-ASYMMETRY AND WEEKLY KAPPA
% -------------------------------------------------------------------------
% Use BAI-2 (signed) to determine direction each week
% Sort by participant and week to ensure consecutive pairing
data_sorted = sortrows(data, {'Participant','Week'});
BAI2_sorted = 2 .* (data_sorted.Add_Dom - data_sorted.Add_NonDom) ./ ...
              (data_sorted.Add_Dom + data_sorted.Add_NonDom) .* 100;

flipCount = 0;
pairCount = 0;
perPlayerFlip = NaN(numel(participants), 1);
weeklyKappas  = NaN(nWeeks-1, 1);

for p = 1:numel(participants)
    pIdx  = data_sorted.Participant == participants(p);
    pVals = BAI2_sorted(pIdx);
    pWeeks= data_sorted.Week(pIdx);
    [pWeeks, ord] = sort(pWeeks);
    pVals = pVals(ord);
    flips = 0; pairs = 0;
    for w = 2:numel(pVals)
        if ~isnan(pVals(w)) && ~isnan(pVals(w-1))
            pairs = pairs + 1;
            if sign(pVals(w)) ~= sign(pVals(w-1))
                flips = flips + 1;
            end
        end
    end
    perPlayerFlip(p) = ternary(pairs > 0, flips/pairs*100, NaN);
    flipCount = flipCount + flips;
    pairCount = pairCount + pairs;
end

overallFlipRate = flipCount / pairCount * 100;
fprintf('\n--- Direction of asymmetry ---\n');
fprintf('  Total consecutive pairs: %d\n', pairCount);
fprintf('  Direction flips:         %d (%.1f%%)\n', flipCount, overallFlipRate);
fprintf('  Per-player flip rate:    %.1f ± %.1f%%\n', ...
    mean(perPlayerFlip,'omitnan'), std(perPlayerFlip,'omitnan'));

% Week-to-week kappa across all players
% For each consecutive week pair, build 2×2 agreement table
for w = 1:nWeeks-1
    dirW  = NaN(numel(participants),1);
    dirW1 = NaN(numel(participants),1);
    for p = 1:numel(participants)
        pIdx = data.Participant == participants(p);
        v1 = T.BAI2(pIdx & data.Week == w);
        v2 = T.BAI2(pIdx & data.Week == (w+1));
        if ~isempty(v1) && ~isempty(v2) && ~isnan(v1) && ~isnan(v2)
            dirW(p)  = double(v1 > 0);
            dirW1(p) = double(v2 > 0);
        end
    end
    ok = ~isnan(dirW) & ~isnan(dirW1);
    if sum(ok) >= 3
        weeklyKappas(w) = computeKappa(dirW(ok), dirW1(ok));
    end
end

meanKappa = mean(weeklyKappas, 'omitnan');
sdKappa   = std(weeklyKappas,  'omitnan');
seKappa   = sdKappa / sqrt(sum(~isnan(weeklyKappas)));
df        = sum(~isnan(weeklyKappas)) - 1;
tCrit     = tinv(0.975, df);
CI_lo     = meanKappa - tCrit * seKappa;
CI_hi     = meanKappa + tCrit * seKappa;
fprintf('  Weekly kappa (mean ± SD): %.2f ± %.2f\n', meanKappa, sdKappa);
fprintf('  95%% CI: [%.2f, %.2f]\n', CI_lo, CI_hi);
fprintf('  Range: %.2f to %.2f\n', min(weeklyKappas,[],'omitnan'), max(weeklyKappas,[],'omitnan'));

% -------------------------------------------------------------------------
%  9. ROLLING DIRECTION CONSISTENCY — FIGURE 6
% -------------------------------------------------------------------------
windowSizes  = [3 5 8];
windowColors = {'b','r','g'};
rolling_consistency = NaN(nWeeks, numel(windowSizes));

for wi = 1:numel(windowSizes)
    wz = windowSizes(wi);
    for w = wz:nWeeks
        wStart = w - wz + 1;
        nConsistent = 0; nValid = 0;
        for p = 1:numel(participants)
            pIdx  = data.Participant == participants(p);
            signs = NaN(wz,1);
            for wk = 1:wz
                v = T.BAI2(pIdx & data.Week == (wStart + wk - 1));
                if ~isempty(v) && ~isnan(v)
                    signs(wk) = sign(v);
                end
            end
            validSigns = signs(~isnan(signs));
            if numel(validSigns) == wz
                nValid = nValid + 1;
                if all(validSigns == validSigns(1))
                    nConsistent = nConsistent + 1;
                end
            end
        end
        if nValid > 0
            rolling_consistency(w, wi) = nConsistent / nValid * 100;
        end
    end
end

fprintf('\n--- Rolling direction consistency ---\n');
for wi = 1:numel(windowSizes)
    fprintf('  %d-week window: %.1f%%\n', windowSizes(wi), ...
        mean(rolling_consistency(:,wi),'omitnan'));
end

figure('Name','Figure 6 - Rolling direction consistency','NumberTitle','off');
hold on;
for wi = 1:numel(windowSizes)
    plot(1:nWeeks, rolling_consistency(:,wi), ...
        [windowColors{wi} '-'], 'LineWidth', 1.8, ...
        'DisplayName', sprintf('%d-week window', windowSizes(wi)));
end
xlabel('Week','FontSize',11);
ylabel('Players with consistent direction (%)','FontSize',11);
title('Rolling direction consistency across the season','FontSize',12);
legend('Location','southwest','FontSize',10);
ylim([0 100]); xlim([1 nWeeks]); grid on; box on;

% -------------------------------------------------------------------------
%  10. SEASON-LONG SPAGHETTI PLOTS: BSA vs BAI-2 — FIGURE 7
% -------------------------------------------------------------------------
figure('Name','Figure 7 - Spaghetti plots','NumberTitle','off');

formulas_spaghetti = {'BSA','BAI2'};
titles_spaghetti   = {'(a) BSA (unsigned)','(b) BAI-2 (signed)'};

for sp = 1:2
    subplot(1,2,sp); hold on;
    fCode = formulas_spaghetti{sp};
    for p = 1:numel(participants)
        pIdx   = data.Participant == participants(p);
        pWeeks = data.Week(pIdx);
        pVals  = T.(fCode)(pIdx);
        [pWeeks, ord] = sort(pWeeks);
        plot(pWeeks, pVals(ord), 'Color', [0.6 0.6 0.6 0.5], 'LineWidth', 0.8);
    end
    yline(10,  'r--', 'LineWidth', 1.2, 'HandleVisibility','off');
    if sp == 2
        yline(-10, 'r--', 'LineWidth', 1.2, 'HandleVisibility','off');
        yline(0,   'k:',  'LineWidth', 1.0, 'HandleVisibility','off');
    end
    xlabel('Week','FontSize',11);
    ylabel('Asymmetry (%)','FontSize',11);
    title(titles_spaghetti{sp},'FontSize',12);
    xlim([1 nWeeks]); grid on; box on;
end

% -------------------------------------------------------------------------
%  11. WEEKLY ASYMMETRY BY FORMULA — FIGURE 1
% -------------------------------------------------------------------------
weeks     = 1:nWeeks;
weekMeans = NaN(nWeeks, nFormulas);
weekSDs   = NaN(nWeeks, nFormulas);

for w = 1:nWeeks
    wIdx = data.Week == w;
    for f = 1:nFormulas
        v = T.(formulaCodes{f})(wIdx);
        weekMeans(w,f) = mean(v,  'omitnan');
        weekSDs(w,f)   = std(v,   'omitnan');
    end
end

figure('Name','Figure 1 - Weekly mean asymmetry','NumberTitle','off');

% Panel a: signed formulas
subplot(2,1,1); hold on;
signedF = find(isSigned);
cmap_s  = lines(numel(signedF));
for fi = 1:numel(signedF)
    f = signedF(fi);
    plot(weeks, weekMeans(:,f), '-', 'Color', cmap_s(fi,:), ...
        'LineWidth', 1.5, 'DisplayName', formulaNames{f});
end
yline(0, 'k:', 'LineWidth', 1, 'HandleVisibility','off');
xlabel('Week','FontSize',10); ylabel('Mean asymmetry (%)','FontSize',10);
title('(a) Signed formulas','FontSize',11);
legend('Location','best','FontSize',9); grid on; box on; xlim([1 nWeeks]);

% Panel b: unsigned formulas
subplot(2,1,2); hold on;
unsignedF = find(~isSigned);
cmap_u    = lines(numel(unsignedF));
for fi = 1:numel(unsignedF)
    f = unsignedF(fi);
    plot(weeks, weekMeans(:,f), '-', 'Color', cmap_u(fi,:), ...
        'LineWidth', 1.5, 'DisplayName', formulaNames{f});
end
yline(10, 'k--', 'LineWidth', 1, 'HandleVisibility','off');
text(1.5, 10.5, '10% threshold', 'FontSize', 8);
xlabel('Week','FontSize',10); ylabel('Mean asymmetry (%)','FontSize',10);
title('(b) Unsigned formulas','FontSize',11);
legend('Location','best','FontSize',9); grid on; box on; xlim([1 nWeeks]);

% -------------------------------------------------------------------------
%  12. PER-PLAYER FLAGGING HEATMAP — FIGURE 3
% -------------------------------------------------------------------------
playerFlagRates = NaN(numel(participants), nFormulas);
for p = 1:numel(participants)
    pIdx = data.Participant == participants(p);
    for f = 1:nFormulas
        absV = abs(T.(formulaCodes{f})(pIdx));
        playerFlagRates(p,f) = mean(absV >= 10, 'omitnan') * 100;
    end
end

figure('Name','Figure 3 - Per-player flagging heatmap','NumberTitle','off');
imagesc(playerFlagRates);
colormap(flipud(hot));
clim([0 100]);
cb = colorbar; cb.Label.String = 'Flagging rate (%)';
set(gca, 'XTick',1:nFormulas, 'XTickLabel',formulaNames, ...
         'YTick',1:numel(participants), ...
         'YTickLabel',arrayfun(@(p) sprintf('P%02d',p), 1:numel(participants),'UniformOutput',false), ...
         'FontSize',10);
title('Per-player observations flagged above 10% (%)','FontSize',12);
xlabel('Formula','FontSize',11); ylabel('Player','FontSize',11);

% -------------------------------------------------------------------------
%  13. INDIVIDUAL vs POSITION-GROUP REPRESENTATIVENESS — FIGURES 8 & 9
% -------------------------------------------------------------------------
groups = [1 2 3 4];

% Use BAI2 for representativeness (main formula recommendation)
fCode = 'BAI2';

% Z-scores relative to position-group weekly mean
zScores   = NaN(height(data),1);
misFlagInd  = false(height(data),1);
misFlagGrp  = false(height(data),1);

for g = groups
    gMembers = unique(data.Participant(data.Group == g));
    if numel(gMembers) < 3, continue; end
    for w = 1:nWeeks
        wgIdx = data.Group == g & data.Week == w;
        vals  = T.(fCode)(wgIdx);
        if sum(~isnan(vals)) < 3, continue; end
        gMean = mean(vals,'omitnan');
        gSD   = std(vals, 'omitnan');
        if gSD == 0, continue; end
        rows  = find(wgIdx);
        for ri = 1:numel(rows)
            r = rows(ri);
            zScores(r) = (T.(fCode)(r) - gMean) / gSD;
            indFlag = abs(T.(fCode)(r)) >= 10;
            grpFlag = abs(gMean) >= 10;
            misFlagInd(r) = indFlag;
            misFlagGrp(r) = grpFlag;
        end
    end
end

validZ     = zScores(~isnan(zScores));
pct1SD     = mean(abs(validZ) > 1) * 100;
pct2SD     = mean(abs(validZ) > 2) * 100;
misClass   = mean(xor(misFlagInd(~isnan(zScores)), misFlagGrp(~isnan(zScores)))) * 100;

fprintf('\n--- Individual vs group representativeness ---\n');
fprintf('  Outside ±1 SD: %.1f%% (expected 31.7%%)\n', pct1SD);
fprintf('  Outside ±2 SD: %.1f%% (expected 4.6%%)\n',  pct2SD);
fprintf('  Mis-classification rate at 10%%: %.1f%%\n',  misClass);

% Figure 8: Z-score heatmap (player × week)
zMatrix = NaN(numel(participants), nWeeks);
for p = 1:numel(participants)
    for w = 1:nWeeks
        idx = data.Participant == participants(p) & data.Week == w;
        if any(idx)
            zMatrix(p,w) = zScores(find(idx,1));
        end
    end
end

figure('Name','Figure 8 - Z-score heatmap','NumberTitle','off');
imagesc(zMatrix, [-3 3]);
colormap(redblue_colormap());
colorbar;
set(gca, 'XTick',1:5:nWeeks, ...
         'YTick',1:numel(participants), ...
         'YTickLabel',arrayfun(@(p) sprintf('P%02d',p), 1:numel(participants),'UniformOutput',false), ...
         'FontSize',10);
xlabel('Week','FontSize',11); ylabel('Player','FontSize',11);
title('Z-scores relative to position-group mean (BAI-2)','FontSize',12);

% Figure 9: Individual vs group trajectories (4 panels)
groupTitles = {'(a) Goalkeepers','(b) Defenders','(c) Midfielders','(d) Forwards (n=1, excluded)'};
figure('Name','Figure 9 - Individual vs group trajectories','NumberTitle','off');

for g = 1:4
    subplot(2,2,g); hold on;
    gMembers = unique(data.Participant(data.Group == g));
    
    % Group mean and SD per week
    gMeanW = NaN(nWeeks,1);
    gSDW   = NaN(nWeeks,1);
    for w = 1:nWeeks
        v = T.(fCode)(data.Group == g & data.Week == w);
        if sum(~isnan(v)) >= 2
            gMeanW(w) = mean(v,'omitnan');
            gSDW(w)   = std(v, 'omitnan');
        end
    end
    
    % Shaded ±1 SD band
    validW = find(~isnan(gMeanW));
    if ~isempty(validW)
        xPatch = [validW; flipud(validW)];
        yPatch = [(gMeanW(validW)+gSDW(validW)); flipud(gMeanW(validW)-gSDW(validW))];
        patch(xPatch, yPatch, [0.6 0.8 1.0], 'FaceAlpha',0.4,'EdgeColor','none');
        plot(validW, gMeanW(validW), 'b-', 'LineWidth',2,'DisplayName','Group mean');
    end
    
    % Individual lines
    for p = 1:numel(gMembers)
        pIdx   = data.Participant == gMembers(p);
        pWeeks = data.Week(pIdx);
        pVals  = T.(fCode)(pIdx);
        [pWeeks,ord] = sort(pWeeks);
        plot(pWeeks, pVals(ord), 'Color',[0.5 0.5 0.5 0.6], ...
            'LineWidth',0.8,'HandleVisibility','off');
    end
    
    yline(10,  'r--','LineWidth',1.2,'HandleVisibility','off');
    yline(-10, 'r--','LineWidth',1.2,'HandleVisibility','off');
    title(groupTitles{g},'FontSize',11);
    xlabel('Week','FontSize',10); ylabel('BAI-2 (%)','FontSize',10);
    xlim([1 nWeeks]); grid on; box on;
end

% -------------------------------------------------------------------------
%  14. PER-PLAYER CUMULATIVE-MEAN STABILISATION — FIGURE 10
% -------------------------------------------------------------------------
% Tolerance: ±1 percentage point around final value
stabilTol   = 1.0;   % ±1 pp
stabilWks   = NaN(numel(participants), nFormulas);

for p = 1:numel(participants)
    pIdx   = data.Participant == participants(p);
    pData  = data(pIdx,:);
    [~,ord]= sort(pData.Week);
    pData  = pData(ord,:);
    
    for f = 1:nFormulas
        vals = abs(T.(formulaCodes{f})(pIdx));
        vals = vals(ord);
        validVals = vals(~isnan(vals));
        nObs = numel(validVals);
        if nObs < 5, continue; end
        
        finalMean = mean(validVals,'omitnan');
        stabIdx   = NaN;
        for n = 3:nObs
            cumMn = mean(validVals(1:n));
            if abs(cumMn - finalMean) <= stabilTol
                % Check it stays within tolerance for all subsequent obs
                staysIn = all(abs(arrayfun(@(k) mean(validVals(1:k)), n:nObs) - finalMean) <= stabilTol);
                if staysIn
                    stabIdx = n;
                    break;
                end
            end
        end
        stabilWks(p,f) = stabIdx;
    end
end

fprintf('\n--- Per-player stabilisation (median [range]) ---\n');
for f = 1:nFormulas
    sv = stabilWks(:,f);
    sv = sv(~isnan(sv));
    fprintf('  %-8s: median %d wk, range %d-%d wk\n', ...
        formulaNames{f}, round(median(sv)), min(sv), max(sv));
end

% Export stabilisation table
stabilTbl = array2table(stabilWks, 'VariableNames', strrep(formulaNames,'-','_'), ...
    'RowNames', arrayfun(@(p) sprintf('P%02d',p), 1:numel(participants),'UniformOutput',false));
writetable(stabilTbl, fullfile(scriptDir,'PaperB_Stabilisation.xlsx'), 'WriteRowNames',true);

% Figure 10: cumulative mean for BSA and BAI-2 (panels a & b)
figure('Name','Figure 10 - Cumulative mean stabilisation','NumberTitle','off');
stabilFormulas = {'BSA','BAI2'};
stabilTitles   = {'(a) BSA/%Diff','(b) BAI-2'};

for sp = 1:2
    subplot(1,2,sp); hold on;
    fCode2  = stabilFormulas{sp};
    fIdx    = strcmp(formulaCodes, fCode2);
    medStab = round(median(stabilWks(:,fIdx),'omitnan'));
    
    for p = 1:numel(participants)
        pIdx   = data.Participant == participants(p);
        vals   = abs(T.(fCode2)(pIdx));
        wks    = data.Week(pIdx);
        [wks,ord] = sort(wks);
        vals   = vals(ord);
        validV = vals(~isnan(vals));
        if numel(validV) < 5, continue; end
        
        cumMeans = arrayfun(@(n) mean(validV(1:n)), 1:numel(validV));
        plot(1:numel(validV), cumMeans, 'Color',[0.7 0.7 0.7],'LineWidth',0.8);
        
        sw = stabilWks(p,fIdx);
        if ~isnan(sw)
            plot(sw, cumMeans(sw), 'ro', 'MarkerSize',6,'MarkerFaceColor','r');
        end
    end
    
    xline(medStab, 'r-','LineWidth',1.5,'DisplayName',sprintf('Median = %d wk',medStab));
    xlabel('Observation number','FontSize',11);
    ylabel('Cumulative mean |asymmetry| (%)','FontSize',11);
    title(stabilTitles{sp},'FontSize',12);
    legend('Location','best','FontSize',9); grid on; box on;
end

% -------------------------------------------------------------------------
%  15. 95% CI HALF-WIDTH ANALYSIS
% -------------------------------------------------------------------------
targetObs  = [4 8 12 20];
MDC        = 6.6;   % Minimal detectable change from Light & Thorborg (2016)
ciResults  = NaN(numel(targetObs), nFormulas);

fprintf('\n--- 95%% CI half-widths (median across players) compared to MDC = %.1f%% ---\n', MDC);
fprintf('%-10s', 'n obs');
for f = 1:nFormulas, fprintf('%8s', formulaNames{f}); end; fprintf('\n');

for ti = 1:numel(targetObs)
    nObs = targetObs(ti);
    for f = 1:nFormulas
        halfWidths = NaN(numel(participants),1);
        for p = 1:numel(participants)
            pIdx  = data.Participant == participants(p);
            vals  = abs(T.(formulaCodes{f})(pIdx));
            wks   = data.Week(pIdx);
            [~,ord] = sort(wks);
            vals  = vals(ord);
            vals  = vals(~isnan(vals));
            if numel(vals) < nObs, continue; end
            subset = vals(1:nObs);
            se     = std(subset,'omitnan') / sqrt(nObs);
            df_    = nObs - 1;
            tCrit2 = tinv(0.975, df_);
            halfWidths(p) = tCrit2 * se;
        end
        ciResults(ti,f) = median(halfWidths,'omitnan');
    end
    fprintf('%-10d', nObs);
    for f = 1:nFormulas
        fprintf('%8.1f', ciResults(ti,f));
    end
    fprintf('\n');
end

% Export CI table
ciTbl = array2table(ciResults, 'VariableNames', strrep(formulaNames,'-','_'), ...
    'RowNames', arrayfun(@(n) sprintf('%d obs',n), targetObs,'UniformOutput',false));
writetable(ciTbl, fullfile(scriptDir,'PaperB_CI_HalfWidths.xlsx'), 'WriteRowNames',true);

% -------------------------------------------------------------------------
%  16. ABDUCTION AND RATIO SUPPLEMENTARY RESULTS
% -------------------------------------------------------------------------
% Abduction BSA
DL_abd  = data.Abd_Dom;
NDL_abd = data.Abd_NonDom;
BSA_abd = abs(DL_abd - NDL_abd) ./ max(DL_abd, NDL_abd) .* 100;
domStronger_abd = mean(DL_abd > NDL_abd,'omitnan') * 100;

% Adduction:Abduction ratio BSA
ratio_dom    = data.Add_Dom    ./ data.Abd_Dom;
ratio_nondom = data.Add_NonDom ./ data.Abd_NonDom;
BSA_ratio    = abs(ratio_dom - ratio_nondom) ./ max(ratio_dom, ratio_nondom) .* 100;
domStronger_ratio = mean(ratio_dom > ratio_nondom,'omitnan') * 100;

fprintf('\n--- Supplementary: abduction and ratio asymmetry (BSA) ---\n');
fprintf('  Abduction BSA:    %.1f ± %.1f%%, dominant stronger in %.1f%%\n', ...
    mean(BSA_abd,'omitnan'), std(BSA_abd,'omitnan'), domStronger_abd);
fprintf('  Add:Abd ratio BSA: %.1f ± %.1f%%, dominant stronger in %.1f%%\n', ...
    mean(BSA_ratio,'omitnan'), std(BSA_ratio,'omitnan'), domStronger_ratio);

fprintf('\n--- All analyses complete ---\n');
fprintf('Outputs written to: %s\n', scriptDir);

% =========================================================================
%  17. HELPER FUNCTIONS
% =========================================================================

function k = computeKappa(a, b)
    % Cohen's kappa for two binary vectors a and b
    a = double(a > 0);
    b = double(b > 0);
    n = numel(a);
    if n < 2, k = NaN; return; end
    po = sum(a == b) / n;
    pA = mean(a); pB = mean(b);
    pe = pA*pB + (1-pA)*(1-pB);
    if pe == 1, k = 1; return; end
    k = (po - pe) / (1 - pe);
end

function out = ternary(cond, a, b)
    % Inline conditional — MATLAB has no ternary operator
    if cond, out = a; else, out = b; end
end

function cmap = redblue_colormap(n)
    % Red-white-blue diverging colormap
    if nargin < 1, n = 256; end
    half = floor(n/2);
    red   = [linspace(0.2,1,half)', linspace(0.2,1,half)', ones(half,1)];
    blue  = [ones(n-half,1), linspace(1,0.2,n-half)', linspace(1,0.2,n-half)'];
    cmap  = [red; blue];
end