%% ======================================================================== % Group normative data and arbitrary asymmetry thresholds fail to % represent individual players in weekly hip strength monitoring % ======================================================================== % Callaway AJ, Dunne C, Williams JM. % % Data: 17 male Premier League players, weekly testing, 2019-20 season % Device: Kangatech KT360 fixed-frame dynamometer % Ethics: Bournemouth University ID 36906 % % This script reproduces all analyses reported in the manuscript. % Place this file in the same folder as Data_SeasonLong_Anonymous.xlsx. % No path editing is required. % % SECTION INDEX % PART 1: Data import & cleaning % PART 2: Derived variable computation % PART 3: Descriptive statistics (group, position, individual) % PART 4: Reference strategy comparison (baseline vs contralateral vs rolling) % PART 5: Threshold sensitivity analysis % PART 6: Rolling CoV stabilisation % PART 7: Match-day proximity analysis % PART 8: Week-to-week variability % PART 9: 3-time-point vs weekly monitoring comparison % PART 10: Individual vs group normative data (includes SWC computation) % PART 11: Linear mixed-effects models (longitudinal, MD, position, trajectories) % PART 12: Figures % PART 13: Export to Excel % % ── Measurement error (Dunne et al., 2022) ────────────────────────────── % Component CoV: Adduction = 7.35%, Abduction = 4.01% % Propagated ratio CoV (delta method: sqrt(CoVa^2 + CoVb^2)): % Add D:ND = 10.39%, Abd D:ND = 5.67%, Ipsilateral = 8.37% % % Smallest Worthwhile Change (SWC): % SWC = 0.2 × pooled within-player SD (raw Newtons) % Pooled within-player SD = sqrt(mean of per-player within-season variances) % This is computed in Part 10a and results in: % SWC Add ND = 12.7 N (within-player SD = 63.3 N) % SWC Add D = 13.0 N (within-player SD = 65.1 N) % SWC Abd ND = 10.3 N (within-player SD = 51.6 N) % SWC Abd D = 11.4 N (within-player SD = 56.8 N) % The corresponding between-player SDs (81.4, 88.0, 77.1, 72.7 N) are % consistently larger, confirming that the within-player SD is the % correct and more conservative denominator for a monitoring threshold. % Note: the LME residual SD (Part 11 / PaperA_ReviewerResponse.m) gives % near-identical values (12.2, 12.7, 10.2, 11.5 N); the small difference % arises because Part 10a includes the week-to-week trend in within-player % variance whereas the LME partials it out. Both are defensible; the % values reported in the manuscript and this header are from Part 10a. % % ── Version history ───────────────────────────────────────────────────── % v3 (June 2026): Corrected SWC documentation in header; added explicit % SWC computation block (Part 10a-ii) with between-player vs within-player % SD table printed to console. No changes to analysis logic. % v2 (prior): Original submitted version. % % Requirements: Statistics and Machine Learning Toolbox (for fitlme) % ======================================================================== clc; clear; close all; %% ======================================================================== % PART 1: DATA IMPORT % ======================================================================== % Locate data file relative to this script — no path editing required scriptDir = fileparts(mfilename('fullpath')); filePath = fullfile(scriptDir, 'Data_SeasonLong_Anonymous.xlsx'); if ~isfile(filePath) error(['Data file not found.\n' ... 'Expected: %s\n' ... 'Place Data_SeasonLong_Anonymous.xlsx in the same folder as this script.'], ... filePath); 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', ... 'AddND_raw','AddD_raw','AbdND_raw','AbdD_raw', ... 'Ratio_AddDvsND','Ratio_AbdDvsND','Ratio_DomAddAbd','Ratio_NDomAddAbd', ... 'Norm_AddND','Norm_AddD','Norm_AbdND','Norm_AbdD'}; 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',''); raw = readtable(filePath, opts, 'UseExcel', false); % Remove padding rows (missing Participant) and week 36 (single observation) raw = raw(~isnan(raw.Participant), :); raw = raw(raw.Week <= 35, :); fprintf('Imported: %d rows, %d players, weeks 1–%d\n', ... height(raw), numel(unique(raw.Participant(~isnan(raw.Participant)))), max(raw.Week)); %% ======================================================================== % PART 2: DERIVED VARIABLES (recomputed — not relying on Excel formulas) % ======================================================================== % Interlimb ratios raw.Ratio_AddDvsND = raw.AddD_raw ./ raw.AddND_raw; % Dom / Non-Dom Adduction raw.Ratio_AbdDvsND = raw.AbdD_raw ./ raw.AbdND_raw; % Dom / Non-Dom Abduction % Ipsilateral ratios raw.Ratio_DomAddAbd = raw.AddD_raw ./ raw.AbdD_raw; % Dom Add / Dom Abd raw.Ratio_NDomAddAbd = raw.AddND_raw ./ raw.AbdND_raw; % Non-Dom Add / Non-Dom Abd % Normalised to body mass (N/kg) raw.Norm_AddND = raw.AddND_raw ./ raw.Mass; raw.Norm_AddD = raw.AddD_raw ./ raw.Mass; raw.Norm_AbdND = raw.AbdND_raw ./ raw.Mass; raw.Norm_AbdD = raw.AbdD_raw ./ raw.Mass; % --- Constants --- TEST_COV_ADD = 7.35; % Dunne et al. (2022) — component measure TEST_COV_ABD = 4.01; % Propagated CoV thresholds for ratio variables (delta method: sqrt(CoVa^2 + CoVb^2)) PROP_COV_ADD_INTERLIMB = 10.39; % Add D:ND: sqrt(7.35^2 + 7.35^2) PROP_COV_ABD_INTERLIMB = 5.67; % Abd D:ND: sqrt(4.01^2 + 4.01^2) PROP_COV_IPSILATERAL = 8.37; % Dom/NDom Add:Abd: sqrt(7.35^2 + 4.01^2) % SWC constants — derived in Part 10a-ii below; verified values hardcoded here % for use in threshold comparisons throughout the script SWC_ADD_ND = 12.7; % N — 0.2 × within-player SD (63.3 N) SWC_ADD_D = 13.0; % N — 0.2 × within-player SD (65.1 N) SWC_ABD_ND = 10.3; % N — 0.2 × within-player SD (51.6 N) SWC_ABD_D = 11.4; % N — 0.2 × within-player SD (56.8 N) PRE_SEASON_WEEKS = 1:3; MAX_WEEK = 35; ROLLING_WINDOW = 4; % weeks % --- Convenience labels --- strengthVars = {'AddND_raw','AddD_raw','AbdND_raw','AbdD_raw'}; strengthLbl = {'Add Non-Dom','Add Dom','Abd Non-Dom','Abd Dom'}; ratioVars = {'Ratio_AddDvsND','Ratio_AbdDvsND','Ratio_DomAddAbd','Ratio_NDomAddAbd'}; ratioLbl = {'Add D:ND','Abd D:ND','Dom Add:Abd','NDom Add:Abd'}; normVars = {'Norm_AddND','Norm_AddD','Norm_AbdND','Norm_AbdD'}; normLbl = {'Norm Add ND','Norm Add D','Norm Abd ND','Norm Abd D'}; allVars = [strengthVars, ratioVars, normVars]; allLbl = [strengthLbl, ratioLbl, normLbl]; groupLabels = {'GK','Defender','Midfield','Forward'}; % Valid rows and unique players validIdx = ~isnan(raw.AddND_raw); players = unique(raw.Participant(validIdx)); nPlayers = numel(players); fprintf('Players with data: %d\n\n', nPlayers); %% ======================================================================== % PART 3: DESCRIPTIVE STATISTICS % ======================================================================== fprintf('============================================================\n'); fprintf(' PART 3: DESCRIPTIVE STATISTICS\n'); fprintf('============================================================\n'); % ---- 3a: Weekly group-level ---- weeklyMean = NaN(MAX_WEEK, numel(allVars)); weeklySD = NaN(MAX_WEEK, numel(allVars)); weeklyN = NaN(MAX_WEEK, numel(allVars)); weeklyCOV = NaN(MAX_WEEK, numel(allVars)); for w = 1:MAX_WEEK wk = raw(raw.Week == w, :); for v = 1:numel(allVars) vals = wk.(allVars{v}); vals(isnan(vals)) = []; weeklyN(w,v) = numel(vals); if numel(vals) >= 2 weeklyMean(w,v) = mean(vals); weeklySD(w,v) = std(vals); weeklyCOV(w,v) = (std(vals)/mean(vals))*100; elseif numel(vals) == 1 weeklyMean(w,v) = vals; end end end % ---- 3b: Individual player profiles ---- pStats = struct(); for p = 1:nPlayers pID = players(p); pData = raw(raw.Participant == pID & validIdx, :); pData = sortrows(pData, 'Week'); pStats(p).ID = pID; pStats(p).Group = pData.Group(1); pStats(p).Position = char(pData.Position(1)); pStats(p).nSess = height(pData); pStats(p).Mass = pData.Mass(1); for v = 1:numel(allVars) vals = pData.(allVars{v}); vals(isnan(vals)) = []; pStats(p).([allVars{v} '_mean']) = mean(vals); pStats(p).([allVars{v} '_sd']) = std(vals); pStats(p).([allVars{v} '_cov']) = (std(vals)/mean(vals))*100; pStats(p).([allVars{v} '_min']) = min(vals); pStats(p).([allVars{v} '_max']) = max(vals); pStats(p).([allVars{v} '_range']) = max(vals) - min(vals); end end % ---- 3c: Position-level season means ---- fprintf('\nPosition-level season averages:\n'); fprintf('%-10s %-18s %8s %8s %6s %5s\n','Position','Variable','Mean','SD','CoV%','n'); for g = 1:4 gData = raw(raw.Group == g & validIdx, :); if isempty(gData), continue; end for v = 1:4 vals = gData.(strengthVars{v}); vals(isnan(vals)) = []; fprintf('%-10s %-18s %8.1f %8.1f %5.1f%% %5d\n', ... groupLabels{g}, strengthLbl{v}, mean(vals), std(vals), ... (std(vals)/mean(vals))*100, numel(vals)); end fprintf('\n'); end %% ======================================================================== % PART 4: REFERENCE STRATEGY COMPARISON % ======================================================================== fprintf('============================================================\n'); fprintf(' PART 4: REFERENCE STRATEGY COMPARISON\n'); fprintf('============================================================\n'); maxRows = sum(validIdx); r_Participant = NaN(maxRows,1); r_Group = NaN(maxRows,1); r_Week = NaN(maxRows,1); r_Obs = NaN(maxRows,4); r_Base = NaN(maxRows,4); r_PctBase = NaN(maxRows,4); r_Ratio = NaN(maxRows,2); r_Roll = NaN(maxRows,4); r_PctRoll = NaN(maxRows,4); r_Flag_Base = NaN(maxRows,2); r_Flag_Contra = NaN(maxRows,2); r_Flag_Roll = NaN(maxRows,2); rc = 0; for p = 1:nPlayers pID = players(p); pAll = raw(raw.Participant == pID, :); pAll = sortrows(pAll, 'Week'); pre = pAll(ismember(pAll.Week, PRE_SEASON_WEEKS) & ~isnan(pAll.AddND_raw), :); baselines = [meanSafe(pre.AddND_raw), meanSafe(pre.AddD_raw), ... meanSafe(pre.AbdND_raw), meanSafe(pre.AbdD_raw)]; for w = 1:MAX_WEEK wRow = pAll(pAll.Week == w, :); if isempty(wRow) || isnan(wRow.AddND_raw(1)), continue; end rc = rc + 1; obs = [wRow.AddND_raw(1), wRow.AddD_raw(1), ... wRow.AbdND_raw(1), wRow.AbdD_raw(1)]; r_Participant(rc) = pID; r_Group(rc) = pAll.Group(1); r_Week(rc) = w; r_Obs(rc,:) = obs; r_Base(rc,:) = baselines; pctBase = ((obs - baselines) ./ baselines) * 100; r_PctBase(rc,:) = pctBase; ratio_add = obs(2) / obs(1); ratio_abd = obs(4) / obs(3); r_Ratio(rc,:) = [ratio_add, ratio_abd]; prevWks = max(1, w-ROLLING_WINDOW):(w-1); prevD = pAll(ismember(pAll.Week, prevWks) & ~isnan(pAll.AddND_raw), :); if height(prevD) >= 2 roll = [meanSafe(prevD.AddND_raw), meanSafe(prevD.AddD_raw), ... meanSafe(prevD.AbdND_raw), meanSafe(prevD.AbdD_raw)]; pctRoll = ((obs - roll) ./ roll) * 100; r_Roll(rc,:) = roll; r_PctRoll(rc,:) = pctRoll; end r_Flag_Base(rc,:) = [flagDrop(pctBase(1),-10), flagDrop(pctBase(2),-10)]; r_Flag_Contra(rc,:) = [flagAsym(ratio_add,0.10), flagAsym(ratio_abd,0.10)]; r_Flag_Roll(rc,:) = [flagDrop(r_PctRoll(rc,1),-10), flagDrop(r_PctRoll(rc,2),-10)]; end end r_Participant = r_Participant(1:rc); r_Group = r_Group(1:rc); r_Week = r_Week(1:rc); r_Obs = r_Obs(1:rc,:); r_Base = r_Base(1:rc,:); r_PctBase = r_PctBase(1:rc,:); r_Ratio = r_Ratio(1:rc,:); r_Roll = r_Roll(1:rc,:); r_PctRoll = r_PctRoll(1:rc,:); r_Flag_Base = r_Flag_Base(1:rc,:); r_Flag_Contra = r_Flag_Contra(1:rc,:); r_Flag_Roll = r_Flag_Roll(1:rc,:); ref = table(r_Participant, r_Group, r_Week, ... r_Obs(:,1), r_Obs(:,2), r_Obs(:,3), r_Obs(:,4), ... r_Base(:,1), r_Base(:,2), r_Base(:,3), r_Base(:,4), ... r_PctBase(:,1), r_PctBase(:,2), r_PctBase(:,3), r_PctBase(:,4), ... r_Ratio(:,1), r_Ratio(:,2), ... r_Roll(:,1), r_Roll(:,2), r_Roll(:,3), r_Roll(:,4), ... r_PctRoll(:,1), r_PctRoll(:,2), r_PctRoll(:,3), r_PctRoll(:,4), ... r_Flag_Base(:,1), r_Flag_Base(:,2), ... r_Flag_Contra(:,1), r_Flag_Contra(:,2), ... r_Flag_Roll(:,1), r_Flag_Roll(:,2), ... 'VariableNames', {'Participant','Group','Week', ... 'Obs_AddND','Obs_AddD','Obs_AbdND','Obs_AbdD', ... 'Base_AddND','Base_AddD','Base_AbdND','Base_AbdD', ... 'PctBase_AddND','PctBase_AddD','PctBase_AbdND','PctBase_AbdD', ... 'Ratio_Add','Ratio_Abd', ... 'Roll_AddND','Roll_AddD','Roll_AbdND','Roll_AbdD', ... 'PctRoll_AddND','PctRoll_AddD','PctRoll_AbdND','PctRoll_AbdD', ... 'Flag_Base_AddND','Flag_Base_AddD', ... 'Flag_Contra_Add','Flag_Contra_Abd', ... 'Flag_Roll_AddND','Flag_Roll_AddD'}); fprintf('Reference table: %d player-week observations\n', height(ref)); valid3 = ref(~isnan(ref.Flag_Base_AddND) & ~isnan(ref.Flag_Roll_AddND), :); nV = height(valid3); nB = sum(valid3.Flag_Base_AddND); nC = sum(valid3.Flag_Contra_Add); nR = sum(valid3.Flag_Roll_AddND); fprintf('\n--- Flagging rates: Adduction Non-Dom (>10%% decline) ---\n'); fprintf(' Observations with all 3 methods available: %d\n', nV); fprintf(' Baseline flags: %3d (%5.1f%%)\n', nB, (nB/nV)*100); fprintf(' Contralateral flags: %3d (%5.1f%%)\n', nC, (nC/nV)*100); fprintf(' Rolling flags: %3d (%5.1f%%)\n', nR, (nR/nV)*100); fB = valid3.Flag_Base_AddND; fC = valid3.Flag_Contra_Add; fR = valid3.Flag_Roll_AddND; allFlag = sum(fB==1 & fC==1 & fR==1); allClear = sum(fB==0 & fC==0 & fR==0); agree = allFlag + allClear; fprintf(' All agree (flag): %3d (%5.1f%%)\n', allFlag, (allFlag/nV)*100); fprintf(' All agree (clear): %3d (%5.1f%%)\n', allClear, (allClear/nV)*100); fprintf(' Total agreement: %3d (%5.1f%%)\n', agree, (agree/nV)*100); fprintf(' Disagreement: %3d (%5.1f%%)\n', nV-agree, ((nV-agree)/nV)*100); fprintf('\n--- Per-player flagging (Add Non-Dom) ---\n'); fprintf('%-6s %-10s %6s %8s %8s %8s\n','ID','Position','nSess','Base%','Contra%','Roll%'); for p = 1:nPlayers pID = players(p); pf = valid3(valid3.Participant == pID, :); if isempty(pf), continue; end n = height(pf); fprintf('%-6d %-10s %6d %7.1f%% %7.1f%% %7.1f%%\n', ... pID, pStats(p).Position, n, ... (sum(pf.Flag_Base_AddND)/n)*100, ... (sum(pf.Flag_Contra_Add)/n)*100, ... (sum(pf.Flag_Roll_AddND)/n)*100); end %% ======================================================================== % PART 5: THRESHOLD SENSITIVITY % ======================================================================== fprintf('\n============================================================\n'); fprintf(' PART 5: THRESHOLD SENSITIVITY\n'); fprintf('============================================================\n'); thresholds = [10, 15, PROP_COV_ADD_INTERLIMB]; thrLbl = {'10%','15%',sprintf('%.1f%% (Propagated CoV)',PROP_COV_ADD_INTERLIMB)}; pctB = ref.PctBase_AddND; pctR = ref.PctRoll_AddND; rAdd = ref.Ratio_Add; fprintf('\nAdduction Non-Dom — flagging by threshold:\n'); fprintf('%-22s %15s %15s %15s\n','Threshold','Baseline','Rolling','Contralateral'); for t = 1:numel(thresholds) th = thresholds(t); nb = sum(pctB < -th, 'omitnan'); tb = sum(~isnan(pctB)); nr = sum(pctR < -th, 'omitnan'); tr = sum(~isnan(pctR)); nc = sum(rAdd < (1-th/100) | rAdd > (1+th/100), 'omitnan'); tc = sum(~isnan(rAdd)); fprintf('%-22s %5d/%d (%4.1f%%) %5d/%d (%4.1f%%) %5d/%d (%4.1f%%)\n', ... thrLbl{t}, nb,tb,(nb/tb)*100, nr,tr,(nr/tr)*100, nc,tc,(nc/tc)*100); end %% ======================================================================== % PART 6: ROLLING CoV STABILISATION % ======================================================================== fprintf('\n============================================================\n'); fprintf(' PART 6: ROLLING CoV STABILISATION\n'); fprintf('============================================================\n'); maxSess = 25; cumCoV = NaN(nPlayers, maxSess, 4); for p = 1:nPlayers pData = raw(raw.Participant == players(p) & validIdx, :); pData = sortrows(pData, 'Week'); nS = height(pData); for s = 2:min(nS, maxSess) for v = 1:4 vals = pData.(strengthVars{v})(1:s); vals(isnan(vals)) = []; if numel(vals) >= 2 cumCoV(p,s,v) = (std(vals)/mean(vals))*100; end end end end fprintf('\nCumulative CoV (mean ± SD across players):\n'); fprintf('%-8s','Sess'); for v = 1:4, fprintf(' %18s',strengthLbl{v}); end fprintf('\n'); for s = 2:20 fprintf('%-8d', s); for v = 1:4 vals = cumCoV(:,s,v); vals(isnan(vals)) = []; if ~isempty(vals) fprintf(' %7.1f ± %5.1f ', mean(vals), std(vals)); else fprintf(' %18s','—'); end end fprintf('\n'); end %% ======================================================================== % PART 7: MATCH-DAY PROXIMITY % ======================================================================== fprintf('\n============================================================\n'); fprintf(' PART 7: MATCH-DAY PROXIMITY\n'); fprintf('============================================================\n'); mdGroups = [0 2 3 5]; mdLbl = {'Pre-Season','MD+2','MD+3','MD+5'}; for v = 1:numel(allVars) fprintf('\n%s:\n', allLbl{v}); fprintf(' %-12s %5s %8s %8s\n','MD','n','Mean','SD'); for m = 1:numel(mdGroups) vals = raw.(allVars{v})(raw.MD_Group == mdGroups(m)); vals(isnan(vals)) = []; fprintf(' %-12s %5d %8.2f %8.2f\n', mdLbl{m}, numel(vals), mean(vals), std(vals)); end end fprintf('\nEffect sizes (Cohen''s d): MD+2 vs MD+3 vs MD+5\n'); fprintf('%-22s %8s %8s %8s\n','Variable','2v3','2v5','3v5'); for v = 1:numel(allVars) v2 = raw.(allVars{v})(raw.MD_Group==2); v2(isnan(v2))=[]; v3 = raw.(allVars{v})(raw.MD_Group==3); v3(isnan(v3))=[]; v5 = raw.(allVars{v})(raw.MD_Group==5); v5(isnan(v5))=[]; fprintf('%-22s %8.2f %8.2f %8.2f\n', allLbl{v}, ... cohensD(v2,v3), cohensD(v2,v5), cohensD(v3,v5)); end %% ======================================================================== % PART 8: WEEK-TO-WEEK VARIABILITY % ======================================================================== fprintf('\n============================================================\n'); fprintf(' PART 8: WEEK-TO-WEEK VARIABILITY\n'); fprintf('============================================================\n'); wk2wk = []; for p = 1:nPlayers pData = raw(raw.Participant == players(p) & validIdx, :); pData = sortrows(pData, 'Week'); if height(pData) < 2, continue; end for s = 2:height(pData) for v = 1:4 prev = pData.(strengthVars{v})(s-1); curr = pData.(strengthVars{v})(s); if ~isnan(prev) && ~isnan(curr) && prev > 0 wk2wk(end+1,:) = [players(p), pData.Week(s), v, ... ((curr-prev)/prev)*100, pData.Week(s)-pData.Week(s-1)]; end end end end fprintf('\nConsecutive-week %% change (gap = 1 week):\n'); fprintf('%-18s %5s %8s %8s %8s %8s\n','Variable','n','Mean%','SD%','Min%','Max%'); for v = 1:4 d = wk2wk(wk2wk(:,3)==v & wk2wk(:,5)==1, 4); fprintf('%-18s %5d %7.1f%% %7.1f%% %7.1f%% %7.1f%%\n', ... strengthLbl{v}, numel(d), mean(d), std(d), min(d), max(d)); end fprintf('\n%% sessions exceeding thresholds (|change|):\n'); fprintf('%-18s %8s %8s %8s\n','Variable','>10%','>15%','>CoV'); for v = 1:4 d = abs(wk2wk(wk2wk(:,3)==v & wk2wk(:,5)==1, 4)); n = numel(d); cv = ifelse(v<=2, TEST_COV_ADD, TEST_COV_ABD); fprintf('%-18s %7.1f%% %7.1f%% %7.1f%%\n', strengthLbl{v}, ... (sum(d>10)/n)*100, (sum(d>15)/n)*100, (sum(d>cv)/n)*100); end %% ======================================================================== % PART 9: THREE-TIME-POINT vs WEEKLY COMPARISON % ======================================================================== fprintf('\n============================================================\n'); fprintf(' PART 9: 3-TIME-POINT vs WEEKLY MONITORING\n'); fprintf('============================================================\n'); TP = [3, 18, 34]; tpLbl = {'TP1 (Wk3)','TP2 (Wk18)','TP3 (Wk34)'}; for v = 1:numel(allVars) tp1vals = raw.(allVars{v})(raw.Week == TP(1)); tp1vals(isnan(tp1vals)) = []; tp1m = mean(tp1vals); flagWks = []; for w = 4:MAX_WEEK wvals = raw.(allVars{v})(raw.Week == w); wvals(isnan(wvals)) = []; if numel(wvals) >= 3 pctChg = ((mean(wvals) - tp1m) / tp1m) * 100; if abs(pctChg) > 10 flagWks(end+1) = w; end end end if isempty(flagWks), wkStr = 'None'; else, wkStr = strjoin(string(flagWks),', '); end fprintf('%-22s: %2d weeks deviate >10%% from TP1 [%s]\n', ... allLbl{v}, numel(flagWks), wkStr); end fprintf('\n--- 3-time-point group means ---\n'); for v = 1:numel(allVars) fprintf('\n%s:\n', allLbl{v}); for t = 1:3 vals = raw.(allVars{v})(raw.Week == TP(t)); vals(isnan(vals)) = []; fprintf(' %-12s Mean=%8.2f SD=%8.2f n=%d\n', tpLbl{t}, mean(vals), std(vals), numel(vals)); end end %% ======================================================================== % PART 10: INDIVIDUAL vs GROUP NORMATIVE DATA % ======================================================================== fprintf('\n============================================================\n'); fprintf(' PART 10: INDIVIDUAL vs GROUP NORMATIVE DATA\n'); fprintf('============================================================\n'); % ---- 10a-i: Between-player vs within-player variance (ICC-style) -------- % % NOTE: This ICC(1,1)-style computation is retained for reference and for % deriving the SWC (see Part 10a-ii below). For the manuscript Table 3 % variance partition coefficients (VPC), use the LME-derived values from % PaperA_ReviewerResponse.m, which are coherent with the fixed-effects % structure of the longitudinal models (Week effect partialled out). % % The Part 10a var_within includes the week-to-week seasonal trend within % each player's data. The LME residual does not (it partials out Week). % This makes Part 10a var_within slightly larger than the LME residual % variance, with the largest difference for abduction variables where the % seasonal trend is most pronounced. Both are defensible for the SWC % denominator; the SWC values in this manuscript use Part 10a (see 10a-ii). fprintf('\n--- 10a-i: Between-player vs Within-player Variance ---\n'); fprintf('%-22s %10s %10s %10s %8s\n', ... 'Variable','Var_Between','Var_Within','Ratio B/W','ICC'); icc_results = NaN(numel(allVars), 4); for v = 1:numel(allVars) playerMeans = NaN(nPlayers,1); playerVars = NaN(nPlayers,1); for p = 1:nPlayers pData = raw(raw.Participant == players(p) & validIdx, :); vals = pData.(allVars{v}); vals(isnan(vals)) = []; if numel(vals) >= 2 playerMeans(p) = mean(vals); playerVars(p) = var(vals); end end playerMeans(isnan(playerMeans)) = []; playerVars(isnan(playerVars)) = []; var_between = var(playerMeans); var_within = mean(playerVars); ratio_bw = var_between / var_within; icc = var_between / (var_between + var_within); icc_results(v,:) = [var_between, var_within, ratio_bw, icc]; fprintf('%-22s %10.2f %10.2f %10.2f %8.3f\n', ... allLbl{v}, var_between, var_within, ratio_bw, icc); end fprintf('\n Interpretation: ICC > 0.75 = most variance is between players.\n'); fprintf(' Ratio B/W > 1 = between-player variance exceeds within-player variance.\n'); fprintf(' See PaperA_ReviewerResponse.m for LME-derived VPC values (Table 3).\n'); % ---- 10a-ii: Explicit SWC computation (Reviewer 2 transparency request) - % % SWC = 0.2 × pooled within-player SD % Pooled within-player SD = sqrt(mean(per-player within-season variance)) % This is equivalent to sqrt(var_within) from the ICC computation above, % but applied to raw strength in Newtons (not normalised N/kg). % % The between-player SD is reported alongside so readers can verify that % between-SD > within-SD (consistent with the variance partition in Table 3) % and that the SWC uses the correct (smaller) denominator. fprintf('\n--- 10a-ii: SWC computation — between-player vs within-player SD ---\n'); fprintf('%-12s %16s %15s %14s %10s\n', ... 'Variable','Between-SD (N)','Within-SD (N)','SD Ratio','SWC (N)'); fprintf('%s\n', repmat('-',1,70)); swcVars = {'AddND_raw','AddD_raw','AbdND_raw','AbdD_raw'}; swcLbl = {'Add ND','Add D','Abd ND','Abd D'}; swc_between_SD = NaN(4,1); swc_within_SD = NaN(4,1); swc_values = NaN(4,1); for v = 1:4 col = swcVars{v}; p_means_v = NaN(nPlayers,1); p_vars_v = NaN(nPlayers,1); for p = 1:nPlayers pData = raw(raw.Participant == players(p) & validIdx, :); vals = pData.(col); vals(isnan(vals)) = []; if numel(vals) >= 2 p_means_v(p) = mean(vals); p_vars_v(p) = var(vals); end end p_means_v(isnan(p_means_v)) = []; p_vars_v(isnan(p_vars_v)) = []; sd_between = std(p_means_v); % SD of per-player season means sd_within = sqrt(mean(p_vars_v)); % pooled within-player SD sd_ratio = sd_between / sd_within; swc = 0.2 * sd_within; swc_between_SD(v) = sd_between; swc_within_SD(v) = sd_within; swc_values(v) = swc; fprintf('%-12s %16.1f %15.1f %14.2f %10.1f\n', ... swcLbl{v}, sd_between, sd_within, sd_ratio, swc); end fprintf('\n Between-SD > Within-SD for all variables (ratio 1.14–1.70×).\n'); fprintf(' SWC = 0.2 × Within-SD: uses the smaller, correct denominator.\n'); fprintf(' A between-SD-based SWC would give larger thresholds (16–18 N),\n'); fprintf(' not the values reported in the manuscript.\n'); % ---- 10b: % sessions outside group norm bands --------------------------- fprintf('\n--- 10b: %% of Sessions Outside Group Normative Bands ---\n'); analyseVars = [strengthVars, ratioVars]; analyseLbl = [strengthLbl, ratioLbl]; outside1SD = NaN(nPlayers, numel(analyseVars)); outside2SD = NaN(nPlayers, numel(analyseVars)); nSessValid = NaN(nPlayers, numel(analyseVars)); for v = 1:numel(analyseVars) vColIdx = find(strcmp(allVars, analyseVars{v})); for p = 1:nPlayers pData = raw(raw.Participant == players(p) & validIdx, :); pData = sortrows(pData, 'Week'); count1 = 0; count2 = 0; nValid = 0; for s = 1:height(pData) w = pData.Week(s); val = pData.(analyseVars{v})(s); if isnan(val) || w > MAX_WEEK, continue; end gMean = weeklyMean(w, vColIdx); gSD = weeklySD(w, vColIdx); if isnan(gMean) || isnan(gSD), continue; end nValid = nValid + 1; if val < (gMean - gSD) || val > (gMean + gSD) count1 = count1 + 1; end if val < (gMean - 2*gSD) || val > (gMean + 2*gSD) count2 = count2 + 1; end end if nValid > 0 outside1SD(p,v) = (count1 / nValid) * 100; outside2SD(p,v) = (count2 / nValid) * 100; nSessValid(p,v) = nValid; end end end fprintf('\nGroup-level summary (mean %% of player-sessions outside band):\n'); fprintf('%-22s %12s %12s\n', 'Variable', 'Outside ±1SD', 'Outside ±2SD'); for v = 1:numel(analyseVars) fprintf('%-22s %10.1f%% %10.1f%%\n', analyseLbl{v}, ... nanmean(outside1SD(:,v)), nanmean(outside2SD(:,v))); end fprintf('\nPer-player detail (Add Non-Dom strength):\n'); fprintf('%-6s %-10s %6s %12s %12s\n','ID','Position','nSess','Outside ±1SD','Outside ±2SD'); for p = 1:nPlayers fprintf('%-6d %-10s %6.0f %10.1f%% %10.1f%%\n', ... players(p), pStats(p).Position, nSessValid(p,1), ... outside1SD(p,1), outside2SD(p,1)); end fprintf('\nPer-player detail (Add Dom:Non-Dom ratio):\n'); fprintf('%-6s %-10s %6s %12s %12s\n','ID','Position','nSess','Outside ±1SD','Outside ±2SD'); vIdx_ratio = find(strcmp(analyseVars, 'Ratio_AddDvsND')); for p = 1:nPlayers fprintf('%-6d %-10s %6.0f %10.1f%% %10.1f%%\n', ... players(p), pStats(p).Position, nSessValid(p,vIdx_ratio), ... outside1SD(p,vIdx_ratio), outside2SD(p,vIdx_ratio)); end % ---- 10c: Asymmetry threshold exceedance (injury-free) ------------------ fprintf('\n--- 10c: Asymmetry Threshold Exceedance (All Players Injury-Free) ---\n'); asymVars = {'Ratio_AddDvsND','Ratio_AbdDvsND','Ratio_DomAddAbd','Ratio_NDomAddAbd'}; asymLbl = {'Add D:ND','Abd D:ND','Dom Add:Abd','NDom Add:Abd'}; fprintf('\nAdd Dom:Non-Dom ratio — per-player threshold exceedance:\n'); fprintf('%-6s %-10s %6s %10s %10s %10s %10s\n', ... 'ID','Position','nSess','Mean Asym%','SD Asym%','>10% Sess','>15% Sess'); for p = 1:nPlayers pData = raw(raw.Participant == players(p) & validIdx, :); ratios = pData.Ratio_AddDvsND; ratios(isnan(ratios)) = []; if isempty(ratios), continue; end asymPct = abs(ratios - 1) * 100; nS = numel(asymPct); fprintf('%-6d %-10s %6d %9.1f%% %9.1f%% %9.1f%% %9.1f%%\n', ... players(p), pStats(p).Position, nS, ... mean(asymPct), std(asymPct), ... (sum(asymPct > 10)/nS)*100, ... (sum(asymPct > 15)/nS)*100); end fprintf('\nGroup summary — %% of ALL sessions exceeding thresholds (injury-free squad):\n'); fprintf('%-22s %6s %10s %10s %10s %10s %10s\n', ... 'Ratio','n','Mean Asym%','SD','>10%','>15%','>CoV'); for v = 1:numel(asymVars) vals = raw.(asymVars{v})(validIdx); vals(isnan(vals)) = []; asymPct = abs(vals - 1) * 100; nA = numel(asymPct); if v == 1, covThr = PROP_COV_ADD_INTERLIMB; elseif v == 2, covThr = PROP_COV_ABD_INTERLIMB; else, covThr = PROP_COV_IPSILATERAL; end fprintf('%-22s %6d %9.1f%% %9.1f%% %9.1f%% %9.1f%% %9.1f%%\n', ... asymLbl{v}, nA, mean(asymPct), std(asymPct), ... (sum(asymPct > 10)/nA)*100, ... (sum(asymPct > 15)/nA)*100, ... (sum(asymPct > covThr)/nA)*100); end % ---- 10d: Individual vs group normative band width ---------------------- fprintf('\n--- 10d: Individual vs Group Normative Band Width ---\n'); fprintf('%-22s %14s %14s %10s\n', 'Variable', 'Group SD (avg)', 'Indiv SD (avg)', 'Ratio G/I'); for v = 1:numel(analyseVars) vColIdx = find(strcmp(allVars, analyseVars{v})); avgGroupSD = nanmean(weeklySD(:, vColIdx)); indivSDs = NaN(nPlayers,1); for p = 1:nPlayers pData = raw(raw.Participant == players(p) & validIdx, :); vals = pData.(analyseVars{v}); vals(isnan(vals)) = []; if numel(vals) >= 2 indivSDs(p) = std(vals); end end avgIndivSD = nanmean(indivSDs); ratio_gi = avgGroupSD / avgIndivSD; fprintf('%-22s %14.2f %14.2f %10.2f\n', ... analyseLbl{v}, avgGroupSD, avgIndivSD, ratio_gi); end fprintf('\n Ratio > 1: group band is wider than individual band.\n'); fprintf(' Group norms are too imprecise for individual monitoring.\n'); %% ======================================================================== % PART 11: LINEAR MIXED-EFFECTS MODELS % ======================================================================== fprintf('\n============================================================\n'); fprintf(' PART 11: LINEAR MIXED-EFFECTS MODELS\n'); fprintf('============================================================\n'); mdl_tbl = raw(validIdx, :); mdl_tbl.Participant = categorical(mdl_tbl.Participant); mdl_tbl.Group = categorical(mdl_tbl.Group); mdl_tbl.Player_Group = categorical(mdl_tbl.Player_Group); mdl_tbl.MD_Group_nom = categorical(mdl_tbl.MD_Group); modVars = [normVars, ratioVars]; modLbl = [normLbl, ratioLbl]; % ---------- Model 1: Longitudinal change over season ---------- fprintf('\n--- MODEL 1: Longitudinal (Week fixed; Participant + MD random) ---\n'); fprintf('%-22s %10s %10s %8s %8s %5s\n', ... 'Variable','Estimate','SE','t','p','Sig'); m1_results = cell(numel(modVars), 6); for v = 1:numel(modVars) dv = modVars{v}; try lme = fitlme(mdl_tbl, ... sprintf('%s ~ Week + (1|Participant) + (1|MD_Group_nom)', dv)); [~, ~, feStats] = fixedEffects(lme, 'DFMethod','satterthwaite'); est = feStats.Estimate(2); se = feStats.SE(2); tval = feStats.tStat(2); pval = feStats.pValue(2); sig = sigStr(pval); fprintf('%-22s %10.6f %10.6f %8.3f %8s %5s\n', ... modLbl{v}, est, se, tval, pStr(pval), sig); m1_results(v,:) = {modLbl{v}, est, se, tval, pval, sig}; catch ME fprintf('%-22s FAILED: %s\n', modLbl{v}, ME.message); m1_results(v,:) = {modLbl{v}, NaN, NaN, NaN, NaN, ''}; end end % ---------- Model 2: Match-day proximity ---------- fprintf('\n--- MODEL 2: Match-day proximity (MD_Group fixed; Participant + Week random) ---\n'); inseason = mdl_tbl(mdl_tbl.MD_Group > 0, :); inseason.MD_Group_nom = categorical(inseason.MD_Group); fprintf('In-season observations: %d\n', height(inseason)); for v = 1:numel(modVars) dv = modVars{v}; try lme = fitlme(inseason, ... sprintf('%s ~ MD_Group_nom + (1|Participant) + (1|Week)', dv)); [~,~,feStats] = fixedEffects(lme, 'DFMethod','satterthwaite'); fprintf('\n %s:\n', modLbl{v}); fprintf(' %-15s %10s %10s %8s %8s\n','Effect','Estimate','SE','t','p'); for r = 1:height(feStats) fprintf(' %-15s %10.4f %10.4f %8.3f %8s\n', ... char(feStats.Name(r)), feStats.Estimate(r), feStats.SE(r), ... feStats.tStat(r), pStr(feStats.pValue(r))); end catch ME fprintf('\n %s: FAILED — %s\n', modLbl{v}, ME.message); end end % ---------- Model 3: Random slopes ---------- fprintf('\n--- MODEL 3: Random slopes — do players change at different rates? ---\n'); fprintf('%-22s %10s %10s %10s %8s\n', 'Variable','AIC_RI','AIC_RS','LRT_p','Sig'); for v = 1:numel(modVars) dv = modVars{v}; try lme_ri = fitlme(mdl_tbl, sprintf('%s ~ Week + (1|Participant)', dv)); lme_rs = fitlme(mdl_tbl, sprintf('%s ~ Week + (Week|Participant)', dv)); comp = compare(lme_ri, lme_rs); pval = comp.pValue(2); sig = sigStr(pval); fprintf('%-22s %10.1f %10.1f %10s %8s\n', ... modLbl{v}, lme_ri.ModelCriterion.AIC, lme_rs.ModelCriterion.AIC, ... pStr(pval), sig); if pval < 0.05 re = randomEffects(lme_rs); slopes = re(2:2:end) + lme_rs.Coefficients.Estimate(2); fprintf(' → Individual slopes range: %.5f to %.5f\n', min(slopes), max(slopes)); end catch ME fprintf('%-22s FAILED: %s\n', modLbl{v}, ME.message); end end % ---------- Model 4: Positional differences ---------- fprintf('\n--- MODEL 4: Positional differences ---\n'); pos_tbl = mdl_tbl(mdl_tbl.Group ~= '4', :); pos_tbl.Player_Group = removecats(pos_tbl.Player_Group); fprintf('Observations (excl. Forward): %d\n', height(pos_tbl)); for v = 1:numel(modVars) dv = modVars{v}; try lme = fitlme(pos_tbl, ... sprintf('%s ~ Player_Group + (1|Participant) + (1|Week)', dv)); [~,~,feStats] = fixedEffects(lme, 'DFMethod','satterthwaite'); fprintf('\n %s:\n', modLbl{v}); fprintf(' %-20s %10s %10s %8s %8s\n','Effect','Estimate','SE','t','p'); for r = 1:height(feStats) fprintf(' %-20s %10.4f %10.4f %8.3f %8s\n', ... char(feStats.Name(r)), feStats.Estimate(r), feStats.SE(r), ... feStats.tStat(r), pStr(feStats.pValue(r))); end catch ME fprintf('\n %s: FAILED — %s\n', modLbl{v}, ME.message); end end % ---------- Model 5: Week × Position interaction ---------- fprintf('\n--- MODEL 5: Week × Position interaction ---\n'); fprintf('%-22s %10s %10s %10s %8s\n','Variable','AIC_main','AIC_int','LRT_p','Sig'); for v = 1:numel(modVars) dv = modVars{v}; try lme_main = fitlme(pos_tbl, ... sprintf('%s ~ Week + Player_Group + (1|Participant)', dv)); lme_int = fitlme(pos_tbl, ... sprintf('%s ~ Week * Player_Group + (1|Participant)', dv)); comp = compare(lme_main, lme_int); pval = comp.pValue(2); fprintf('%-22s %10.1f %10.1f %10s %8s\n', ... modLbl{v}, lme_main.ModelCriterion.AIC, lme_int.ModelCriterion.AIC, ... pStr(pval), sigStr(pval)); catch ME fprintf('%-22s FAILED: %s\n', modLbl{v}, ME.message); end end %% ======================================================================== % PART 12: FIGURES % ======================================================================== fprintf('\n============================================================\n'); fprintf(' PART 12: FIGURES\n'); fprintf('============================================================\n'); colGrey = [0.80 0.80 0.80]; colBlack = [0 0 0]; %#ok colRed = [0.85 0.15 0.15]; %#ok LW_THRESH = 2.0; LW_MEAN = 1.5; LW_INDIV = 0.9; MS_MEAN = 3; figDir = fullfile(scriptDir, 'Figures'); if ~exist(figDir, 'dir'), mkdir(figDir); end % Figure 1: Weekly mean ± 95% CI fig1 = figure('Name','Fig1_WeeklyMeans','Position',[50 50 1800 800],'Color','w'); tiledlayout(2, 4, 'TileSpacing','compact','Padding','compact'); normIdxMap = [9 10 11 12]; normYlbl = {'Norm Add ND (N/kg)','Norm Add D (N/kg)','Norm Abd ND (N/kg)','Norm Abd D (N/kg)'}; panelLbl = {'(a)','(b)','(c)','(d)','(e)','(f)','(g)','(h)'}; for v = 1:4 nexttile; vIdx = normIdxMap(v); m = weeklyMean(:,vIdx); s = weeklySD(:,vIdx); n = weeklyN(:,vIdx); sem = s ./ sqrt(n); x = (1:MAX_WEEK)'; ok = ~isnan(m); fill([x(ok); flipud(x(ok))], [m(ok)+1.96*sem(ok); flipud(m(ok)-1.96*sem(ok))], ... colGrey, 'EdgeColor','none','FaceAlpha',0.5); hold on; plot(x(ok), m(ok), 'k-o','MarkerSize',MS_MEAN,'LineWidth',LW_MEAN,'MarkerFaceColor','k'); xline(3.5,'--r','LineWidth',LW_THRESH); xlabel('Week','FontSize',9); ylabel(normYlbl{v},'FontSize',9); text(0.02,0.97,panelLbl{v},'Units','normalized','FontSize',10,'FontWeight','bold','VerticalAlignment','top'); xlim([0 36]); set(gca,'Color','w','FontSize',9); hold off; end ratioYlbl = {'Add D:ND','Abd D:ND','Dom Add:Abd','NDom Add:Abd'}; for v = 1:4 nexttile; vIdx = v + 4; m = weeklyMean(:,vIdx); s = weeklySD(:,vIdx); n = weeklyN(:,vIdx); sem = s ./ sqrt(n); x = (1:MAX_WEEK)'; ok = ~isnan(m); fill([x(ok); flipud(x(ok))], [m(ok)+1.96*sem(ok); flipud(m(ok)-1.96*sem(ok))], ... colGrey, 'EdgeColor','none','FaceAlpha',0.5); hold on; plot(x(ok), m(ok), 'k-o','MarkerSize',MS_MEAN,'LineWidth',LW_MEAN,'MarkerFaceColor','k'); yline(1.0,':k','LineWidth',LW_THRESH); yline(1.1,':r','LineWidth',LW_THRESH); yline(0.9,':r','LineWidth',LW_THRESH); xline(3.5,'--r','LineWidth',LW_THRESH); xlabel('Week','FontSize',9); ylabel(ratioYlbl{v},'FontSize',9); text(0.02,0.97,panelLbl{v+4},'Units','normalized','FontSize',10,'FontWeight','bold','VerticalAlignment','top'); xlim([0 36]); set(gca,'Color','w','FontSize',9); hold off; end saveFig(fig1, figDir, 'Fig1_WeeklyMeans'); % Figure 2: Individual trajectories fig2 = figure('Name','Fig2_IndivTrajectories','Position',[50 50 1200 500],'Color','w'); cmap2 = lines(nPlayers); tiledlayout(1, 2, 'TileSpacing','compact','Padding','compact'); nexttile; hold on; for p = 1:nPlayers pD = sortrows(raw(raw.Participant == players(p) & validIdx, :),'Week'); plot(pD.Week, pD.AddND_raw, '-', 'Color',[cmap2(p,:) 0.5],'LineWidth',LW_INDIV); end xlabel('Week'); ylabel('Force (N)'); text(0.02,0.97,'(a)','Units','normalized','FontSize',10,'FontWeight','bold','VerticalAlignment','top'); xlim([0 36]); set(gca,'Color','w'); hold off; nexttile; hold on; for p = 1:nPlayers pD = sortrows(raw(raw.Participant == players(p) & validIdx, :),'Week'); r = pD.Ratio_AddDvsND; wk = pD.Week; ok = ~isnan(r); plot(wk(ok), r(ok), '-', 'Color',[cmap2(p,:) 0.5],'LineWidth',LW_INDIV); end yline(1.0,':k','LineWidth',LW_THRESH); yline(1.1,':r','LineWidth',LW_THRESH); yline(0.9,':r','LineWidth',LW_THRESH); xlabel('Week'); ylabel('Ratio'); text(0.02,0.97,'(b)','Units','normalized','FontSize',10,'FontWeight','bold','VerticalAlignment','top'); xlim([0 36]); set(gca,'Color','w'); hold off; saveFig(fig2, figDir, 'Fig2_IndivTrajectories'); % Figure 3: Individual trajectories vs group normative bands fig3 = figure('Name','Fig3_NormBands','Position',[50 50 1400 900],'Color','w'); plotVarsF3 = {'AddND_raw','AddD_raw','Ratio_AddDvsND','Ratio_DomAddAbd'}; plotLblF3 = {'Add Non-Dom (N)','Add Dom (N)','Add Dom:Non-Dom Ratio','Dom Add:Abd Ratio'}; panelLblF3 = {'(a)','(b)','(c)','(d)'}; tiledlayout(2, 2, 'TileSpacing','compact','Padding','compact'); cmap3 = lines(nPlayers); for v = 1:4 nexttile; vColIdx = find(strcmp(allVars, plotVarsF3{v})); m = weeklyMean(:,vColIdx); s = weeklySD(:,vColIdx); x = (1:MAX_WEEK)'; ok = ~isnan(m) & ~isnan(s); fill([x(ok); flipud(x(ok))], [m(ok)+2*s(ok); flipud(m(ok)-2*s(ok))], ... [0.90 0.90 0.90], 'EdgeColor','none','FaceAlpha',0.6); hold on; fill([x(ok); flipud(x(ok))], [m(ok)+s(ok); flipud(m(ok)-s(ok))], ... [0.75 0.75 0.75], 'EdgeColor','none','FaceAlpha',0.6); plot(x(ok), m(ok), 'k-', 'LineWidth', 2); for p = 1:nPlayers pD = sortrows(raw(raw.Participant == players(p) & validIdx, :),'Week'); vals = pD.(plotVarsF3{v}); wks = pD.Week; okp = ~isnan(vals); plot(wks(okp), vals(okp), '-', 'Color',[cmap3(p,:) 0.4], 'LineWidth', LW_INDIV); end if v >= 3 yline(1.0,':k','LineWidth',LW_THRESH); yline(1.1,':r','LineWidth',LW_THRESH); yline(0.9,':r','LineWidth',LW_THRESH); end xlabel('Week'); ylabel(plotLblF3{v}); text(0.02,0.97,panelLblF3{v},'Units','normalized','FontSize',10,'FontWeight','bold','VerticalAlignment','top'); xlim([0 36]); set(gca,'Color','w'); hold off; end saveFig(fig3, figDir, 'Fig3_NormBands'); % Figure 4: Asymmetry exceedance heatmap fig4 = figure('Name','Fig4_AsymHeatmap','Position',[50 50 900 700],'Color','w'); asymHeat = NaN(nPlayers, numel(asymVars)); playerLabels9 = cell(nPlayers, 1); for p = 1:nPlayers playerLabels9{p} = sprintf('P%d (%s)', players(p), pStats(p).Position); for v = 1:numel(asymVars) pData = raw(raw.Participant == players(p) & validIdx, :); vals = pData.(asymVars{v}); vals(isnan(vals)) = []; if ~isempty(vals) asymHeat(p,v) = (sum(abs(vals-1)*100 > 10) / numel(vals)) * 100; end end end imagesc(asymHeat); colormap(flipud(hot)); colorbar; caxis([0 100]); set(gca,'YTick',1:nPlayers,'YTickLabel',playerLabels9, ... 'XTick',1:numel(asymVars),'XTickLabel',asymLbl,'Color','w','FontSize',8); xlabel('Ratio Type'); ylabel('Player'); saveFig(fig4, figDir, 'Fig4_AsymHeatmap'); % Figure 5: Per-player flagging fig5 = figure('Name','Fig5_Flagging','Position',[50 50 1300 500],'Color','w'); flagPcts = []; pIDs = []; for p = 1:nPlayers pID = players(p); pf = valid3(valid3.Participant == pID, :); if height(pf) < 5, continue; end n = height(pf); flagPcts(end+1,:) = [ ... (sum(pf.Flag_Base_AddND)/n)*100, ... (sum(pf.Flag_Contra_Add)/n)*100, ... (sum(pf.Flag_Roll_AddND)/n)*100]; pIDs(end+1) = pID; end bar(flagPcts); set(gca,'XTickLabel',string(pIDs),'Color','w'); xlabel('Player ID'); ylabel('% Sessions Flagged'); legend({'Baseline','Contralateral','Rolling'},'Location','northeastoutside'); ylim([0 100]); saveFig(fig5, figDir, 'Fig5_FlaggingComparison'); % Figure 6: Baseline vs Rolling scatter fig6 = figure('Name','Fig6_Scatter','Position',[50 50 600 500],'Color','w'); pB2 = ref.PctBase_AddND; pR2 = ref.PctRoll_AddND; wk2 = ref.Week; ok2 = ~isnan(pB2) & ~isnan(pR2); scatter(pB2(ok2), pR2(ok2), 20, wk2(ok2), 'filled','MarkerFaceAlpha',0.6); colormap(parula); cb = colorbar; ylabel(cb,'Week'); hold on; plot([-60 60],[-60 60],'--k','LineWidth',LW_THRESH); xline(-10,':r','LineWidth',LW_THRESH); yline(-10,':r','LineWidth',LW_THRESH); xlabel('% Change from Baseline'); ylabel('% Change from Rolling Avg'); axis equal; axis([-60 60 -60 60]); set(gca,'Color','w'); hold off; saveFig(fig6, figDir, 'Fig6_BaselineVsRolling'); fprintf('\nAll figures saved to: %s/\n', figDir); %% ======================================================================== % PART 13: EXPORT TO EXCEL % ======================================================================== fprintf('\n============================================================\n'); fprintf(' PART 13: EXPORT\n'); fprintf('============================================================\n'); outFile = fullfile(scriptDir, 'PaperA_Results.xlsx'); wkTbl = array2table([(1:MAX_WEEK)', weeklyN(:,1:4), weeklyMean(:,1:4), weeklySD(:,1:4), ... weeklyMean(:,5:8), weeklySD(:,5:8), weeklyMean(:,9:12), weeklySD(:,9:12)], ... 'VariableNames', ['Week', ... strcat(strengthVars,'_n'), strcat(strengthVars,'_mean'), strcat(strengthVars,'_sd'), ... strcat(ratioVars,'_mean'), strcat(ratioVars,'_sd'), ... strcat(normVars,'_mean'), strcat(normVars,'_sd')]); writetable(wkTbl, outFile, 'Sheet','WeeklyDescriptives'); pTbl = struct2table(pStats); writetable(pTbl, outFile, 'Sheet','PlayerProfiles'); writetable(ref, outFile, 'Sheet','ReferenceComparison'); m1Tbl = cell2table(m1_results, ... 'VariableNames',{'Variable','Estimate','SE','t','p','Sig'}); writetable(m1Tbl, outFile, 'Sheet','Model1_Longitudinal'); % SWC summary sheet (new in v3 — for Reviewer 2 transparency) swcTbl = table(swcLbl', swc_between_SD, swc_within_SD, ... swc_between_SD ./ swc_within_SD, swc_values, ... 'VariableNames',{'Variable','Between_SD_N','Within_SD_N','SD_Ratio','SWC_N'}); writetable(swcTbl, outFile, 'Sheet','SWC_Computation'); fprintf('Exported: %s\n', outFile); fprintf('\n=== ALL ANALYSES COMPLETE ===\n'); %% ======================================================================== % HELPER FUNCTIONS % ======================================================================== function saveFig(figHandle, folder, name) set(figHandle, 'Color','w','InvertHardcopy','off'); pngFile = fullfile(folder, [name '.png']); tiffFile = fullfile(folder, [name '.tiff']); if exist('exportgraphics','file') exportgraphics(figHandle, pngFile, 'Resolution',300); exportgraphics(figHandle, tiffFile, 'Resolution',300); else print(figHandle, pngFile, '-dpng', '-r300'); print(figHandle, tiffFile, '-dtiff', '-r300'); end fprintf(' Saved: %s.png & .tiff\n', name); end function m = meanSafe(x) x(isnan(x)) = []; if isempty(x), m = NaN; else, m = mean(x); end end function f = flagDrop(pctChange, threshold) if isnan(pctChange), f = NaN; else, f = double(pctChange < threshold); end end function f = flagAsym(ratio, tolerance) if isnan(ratio), f = NaN; else, f = double(ratio < (1-tolerance) | ratio > (1+tolerance)); end end function d = cohensD(x1, x2) n1 = numel(x1); n2 = numel(x2); poolSD = sqrt(((n1-1)*std(x1)^2 + (n2-1)*std(x2)^2) / (n1+n2-2)); if poolSD == 0, d = 0; else, d = (mean(x1)-mean(x2)) / poolSD; end end function s = pStr(p) if p < 0.001, s = '<0.001'; else, s = sprintf('%.4f', p); end end function s = sigStr(p) if p < 0.001, s = '***'; elseif p < 0.01, s = '**'; elseif p < 0.05, s = '*'; else, s = ''; end end function r = ifelse(cond, a, b) if cond, r = a; else, r = b; end end