%% === Count number of animals per Sex × Exposure group (High level, Adult only) ===

T = UpdatedFinalTable;  % for convenience

% 1) Convert key variables to string
T.SubjectID = string(T.SubjectID);
T.Sex       = string(T.Sex);
T.Exposure  = string(T.Exposure);
T.Level     = string(T.Level);
T.Age       = string(T.Age);

% 2) Subset to Level = "High" and Age = "Adult"
idx = (T.Level == "High") & (T.Age == "Adult");
Tsub = T(idx, :);

% 3) Define experimental group as Sex_Exposure
Tsub.ExperimentalGroup = Tsub.Sex + "_" + Tsub.Exposure;

% 4) Unique SubjectID × ExperimentalGroup
varsForUniq = ["SubjectID", "ExperimentalGroup"];
T_unique = unique(Tsub(:, varsForUniq), "rows");

% 5) Count unique animals per Sex × Exposure group
animalCounts = groupcounts(T_unique, "ExperimentalGroup");

% 6) Rename only the count column
animalCounts.Properties.VariableNames{end} = 'N_Animals';

disp(animalCounts);

%% === Unique neurons per ExperimentalGroup × DIVISION (High, Adult only) ===

T = UpdatedFinalTable;

% 1) Ensure appropriate types
T.SubjectID = string(T.SubjectID);
T.Sex       = string(T.Sex);
T.Exposure  = string(T.Exposure);
T.Level     = string(T.Level);
T.Age       = string(T.Age);
T.DIVISION  = string(T.DIVISION);
T.UNIT      = string(T.UNIT);

% 2) Filter to Level = "High" and Age = "Adult"
idx  = (T.Level == "High") & (T.Age == "Adult");
Tsub = T(idx, :);

% 3) Define experimental group as Sex_Exposure
Tsub.ExperimentalGroup = Tsub.Sex + "_" + Tsub.Exposure;

% 4) Each neuron is uniquely defined by:
%    SubjectID × UNIT × DIVISION × ExperimentalGroup
varsForNeuron   = ["SubjectID", "UNIT", "DIVISION", "ExperimentalGroup"];
T_uniqueNeurons = unique(Tsub(:, varsForNeuron), "rows");

% 5) Count neurons per ExperimentalGroup × DIVISION
neuronCounts = groupcounts(T_uniqueNeurons, ["ExperimentalGroup", "DIVISION"]);

% 6) Rename the last column (GroupCount) to N_Neurons
neuronCounts.Properties.VariableNames{end} = 'N_Neurons';

disp(neuronCounts);

%% === Median‐based iMM, iRS, iPE per ExperimentalGroup × DIVISION ===
% Uses frequency-pair entries (NO neuron merging)

T = UpdatedFinalTable;

% 1) Ensure appropriate types
T.SubjectID = string(T.SubjectID);
T.Sex       = string(T.Sex);
T.Exposure  = string(T.Exposure);
T.Level     = string(T.Level);
T.Age       = string(T.Age);
T.DIVISION  = string(T.DIVISION);

% 2) Filter to Level = "High" and Age = "Adult"
idx  = (T.Level == "High") & (T.Age == "Adult");
Tsub = T(idx, :);

% If you want only URETHANE:
% Tsub.PREPARATION = string(Tsub.PREPARATION);
% Tsub = Tsub(Tsub.PREPARATION == "URETHANE", :);

% 3) Define experimental group as Sex_Exposure
Tsub.ExperimentalGroup = Tsub.Sex + "_" + Tsub.Exposure;

% 4) Compute medians of NOR_DEV, NOR_CTR, NOR_STD per ExperimentalGroup × Division
varsToMedian = ["NOR_DEV", "NOR_CTR", "NOR_STD"];

medianTbl = groupsummary( ...
    Tsub, ...
    ["ExperimentalGroup", "DIVISION"], ...   % << NOW CORRECTLY SEGREGATED BY DIVISION
    "median", ...
    varsToMedian);

% groupsummary creates:
% GroupCount, median_NOR_DEV, median_NOR_CTR, median_NOR_STD

% 5) Derive median iMM, iRS, iPE
medianTbl.Median_iMM = medianTbl.median_NOR_DEV - medianTbl.median_NOR_STD;  % DEV - STD
medianTbl.Median_iRS = medianTbl.median_NOR_CTR - medianTbl.median_NOR_STD;  % CTR - STD
medianTbl.Median_iPE = medianTbl.median_NOR_DEV - medianTbl.median_NOR_CTR;  % DEV - CTR

% 6) Optional: Clean table (use cell array of char vectors)
median_i_components = medianTbl(:, { ...
    'ExperimentalGroup', 'DIVISION', ...
    'median_NOR_DEV', 'median_NOR_CTR', 'median_NOR_STD', ...
    'Median_iMM', 'Median_iRS', 'Median_iPE'});

disp(median_i_components);


%% Merge high level tables from GRIN2B animals

% Combine all manually collected tables
finalTable = [T1; T2; T3; T4];  % List all your T-tables

%% %% Check experimental groups

% Assuming finalTable is already loaded

% Extract all unique values from the 'PREPARATION' column
experimentalGroups = unique(finalTable.PREPARATION);

% Display the unique experimental groups
disp('Unique Experimental Groups:');
disp(experimentalGroups);

%% Generate separate columns for each Between-Subject Variables (Grouping Factors of sex, age and exposure)

% Assuming finalTable is already loaded

% Initialize new columns for Sex, Age, and Exposure
Sex = cell(height(finalTable), 1);
Age = cell(height(finalTable), 1);
Exposure = cell(height(finalTable), 1);

% Loop through each row of the finalTable to extract the information from 'PREPARATION'
for i = 1:height(finalTable)
    % Get the preparation name for the current row
    preparationName = finalTable.PREPARATION{i};
    
    % Handle the special case for "URETHANE" in the control condition
    if strcmp(preparationName, 'URETHANE')
        Sex{i} = 'Female';
        Age{i} = 'Adult';
        Exposure{i} = 'Control';
    else
        % Determine Sex based on the preparation name
        if contains(preparationName, 'Female')
            Sex{i} = 'Female';
        elseif contains(preparationName, 'Male')
            Sex{i} = 'Male';
        end
        
        % Determine Age based on the preparation name
        if contains(preparationName, 'Adolescente')
            Age{i} = 'Prepubertal';
        elseif contains(preparationName, 'Adult')
            Age{i} = 'Adult';
        end
        
        % Determine Exposure based on the preparation name
        if contains(preparationName, 'Control')
            Exposure{i} = 'Control';
        else
            % If the name doesn't contain 'Control', assume GRIN2B 
            Exposure{i} = 'GRIN2B';
        end
    end
end

% Add the new columns to the finalTable
finalTable.Sex = Sex;
finalTable.Age = Age;
finalTable.Exposure = Exposure;

% Display the updated table with new columns
disp('Updated Final Table with Sex, Age, and Exposure:');
disp(finalTable);

% Save the updated table as a .mat file
save('UpdatedfinalTable.mat', 'finalTable');

%% Check if finalTable and UpdatedfinalTable are equal

% Load finalTable and UpdatedfinalTable.mat files
load('finalTable.mat', 'finalTable'); 
load('UpdatedfinalTable.mat', 'UpdatedfinalTable'); 

% Step 1: Extract 'ANIMAL' columns from both tables
finalTableAnimal = finalTable.ANIMAL;
UpdatedfinalTableAnimal = UpdatedfinalTable.ANIMAL;

% Step 2: Count the number of unique 'ANIMAL' entries in both tables
numUniquefinalTableAnimals = numel(unique(finalTableAnimal));
numUniqueUpdatedfinalTableAnimals = numel(unique(UpdatedfinalTableAnimal));

% Step 3: Compare the number of unique 'ANIMAL' entries
if numUniquefinalTableAnimals == numUniqueUpdatedfinalTableAnimals
    disp('The number of unique ANIMAL entries are equal in both tables.');
else
    disp('The number of unique ANIMAL entries are NOT equal in both tables.');
    disp(['finalTable has ', num2str(numUniquefinalTableAnimals), ' unique ANIMAL entries.']);
    disp(['UpdatedfinalTable has ', num2str(numUniqueUpdatedfinalTableAnimals), ' unique ANIMAL entries.']);
end

%% Check for, e.g., nSFR similarity between finalTable and UpdatedfinalTable
 
% Load finalTable and UpdatedfinalTable.mat files
load('finalTable.mat', 'finalTable'); 
load('UpdatedfinalTable.mat', 'UpdatedfinalTable'); 

% Step 1: Extract 'ANIMAL' and 'nSFR' columns from both tables
finalTableAnimal = finalTable.ANIMAL;
UpdatedfinalTableAnimal = UpdatedfinalTable.ANIMAL;
finalTableNSFR = finalTable.nSFR;
UpdatedfinalTableNSFR = UpdatedfinalTable.nSFR;

% Step 2: Sort both tables by 'ANIMAL'
[~, idxFinal] = sort(finalTableAnimal);
[~, idxUpdated] = sort(UpdatedfinalTableAnimal);

finalTableSorted = finalTable(idxFinal, :);
UpdatedfinalTableSorted = UpdatedfinalTable(idxUpdated, :);

% Step 3: Compare the 'nSFR' values for each corresponding 'ANIMAL'
if isequal(finalTableSorted.nSFR, UpdatedfinalTableSorted.nSFR)
    disp('The nSFR values are equal for each ANIMAL in both tables.');
else
    disp('The nSFR values are NOT equal for each ANIMAL in both tables.');
    
    % Step 4: Identify differences in nSFR values
    diffRows = find(finalTableSorted.nSFR ~= UpdatedfinalTableSorted.nSFR);
    
    % Step 5: Display the differences
    disp('Rows with different nSFR values:');
    disp(diffRows);
end

%% Change "ANIMAL" header to "SubjectID"

% Load UpdatedfinalTable.mat file
load('UpdatedfinalTable.mat', 'UpdatedfinalTable'); 

% Step 1: Rename the 'ANIMAL' column to 'SubjectID'
UpdatedfinalTable.Properties.VariableNames{'ANIMAL'} = 'SubjectID';

% Step 2: Save the updated table back to a .mat file (optional)
save('UpdatedfinalTable.mat', 'UpdatedfinalTable');

% Confirm the change
disp('Column name "ANIMAL" has been changed to "SubjectID" in UpdatedfinalTable.');

%% Add 'High label'

UpdatedfinalTable.Level = repmat({'High'}, height(UpdatedfinalTable), 1);

%% Merge tables to the definitive one

% Ensure both tables have the same column order
T1 = UpdatedFinalTableOriginal;          % e.g., from 'Low' level
T2 = UpdatedfinalTable;   % e.g., from 'High' level

% Step 1: Add missing 'Group' column to T2 if needed
if ~ismember('Group', T2.Properties.VariableNames)
    T2.Group = repmat({'Unknown'}, height(T2), 1);  % or use actual label if known
end

% Step 2: Reorder columns of T2 to match T1
T2 = T2(:, T1.Properties.VariableNames);

% Step 3: Concatenate the two tables
MergedTable = [T1; T2];

%% We have distincts ages and stimulation levels, we must address only adults under high stimulation 


%% ---------------- GLMM for nSFR (Normal Distribution) with Aggregation ------------------ %%
% Revised for current experiment: Sex × Exposure (Genotype) × Division

%% Step 1: Load Data
data = load('UpdatedFinalTable.mat');

% If the .mat file contains a struct with the table inside:
if isstruct(data) && isfield(data, 'UpdatedFinalTable')
    tbl = data.UpdatedFinalTable;
elseif istable(data)
    tbl = data;
else
    error('Cannot find variable "UpdatedFinalTable" in the loaded .mat file.');
end

%% Step 2: Convert Categorical Variables
tbl.Sex = categorical(tbl.Sex);              
tbl.Exposure = categorical(tbl.Exposure);    
tbl.DIVISION = categorical(tbl.DIVISION);    

if iscell(tbl.SubjectID)
    tbl.SubjectID = categorical(tbl.SubjectID);
end

%% Step 2.5: Filter to include only Adult animals under High stimulation
% Check if required columns exist
if ~all(ismember({'Age', 'Level'}, tbl.Properties.VariableNames))
    error('Required columns "Age" and "Level" not found in the table.');
end

tbl = tbl(strcmp(tbl.Age, 'Adult') & strcmp(tbl.Level, 'High'), :);

%% Step 3: Aggregate by Subject × Division
[G, subjectList, divisionList] = findgroups(tbl.SubjectID, tbl.DIVISION);
median_nSFR = splitapply(@median, tbl.nSFR, G);

T_agg = table(subjectList, divisionList, median_nSFR, ...
    'VariableNames', {'SubjectID', 'DIVISION', 'nSFR'});


%% Step 5: Add Subject-level Metadata (Sex, Exposure)
T_subject = unique(tbl(:, {'SubjectID', 'Sex', 'Exposure'}));  % One row per subject
T_agg = outerjoin(T_agg, T_subject, 'Keys', 'SubjectID', 'MergeKeys', true);

%% Step 6: Fit GLMM (Normal, Identity Link, Laplace Fit)
glme = fitglme(T_agg, ...
    'nSFR ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');

%% Step 7: Display Model Results
disp(glme);

%% Get 95% profile likelihood confidence intervals
ci = coefCI(glme, 'Alpha', 0.05);  % 95% CI = 1 - 0.05
% Display labeled profile likelihood CIs
fprintf('\nProfile Likelihood 95%% Confidence Intervals:\n');
for i = 1:length(glme.CoefficientNames)
    fprintf('%s: [%.3f, %.3f]\n', glme.CoefficientNames{i}, ci(i, 1), ci(i, 2));
end

%% Step 8: AIC, BIC, LogLikelihood
AIC = glme.ModelCriterion.AIC;
BIC = glme.ModelCriterion.BIC;
LogLikelihood = glme.LogLikelihood;

disp(['AIC: ', num2str(AIC)]);
disp(['BIC: ', num2str(BIC)]);
disp(['Log Likelihood: ', num2str(LogLikelihood)]);

%% Step 9: Compute ICC
[psi, mse] = covarianceParameters(glme);
var_intercept = psi{1}(1,1);  % Random intercept variance (SubjectID)
var_residual = mse;           % Residual variance
ICC = var_intercept / (var_intercept + var_residual);
fprintf('Intraclass Correlation (ICC): %.3f\n', ICC);

%% Step 10: Residual Diagnostics
residuals = glme.Residuals.Raw;

figure;
histogram(residuals, 30);
xlabel('Residuals');
ylabel('Frequency');
title('Histogram of Residuals');
grid on;

figure;
qqplot(residuals);
title('QQ Plot of Residuals');
grid on;

%% Step 11: K-Fold Cross-Validation Setup for MedAE
K = 10;  % Specify number of folds for cross-validation
cv = cvpartition(size(T_agg, 1), 'KFold', K);  % Create K-fold partition object

medae_kfold = zeros(K, 1);  % Store MedAE for each fold

%% Step 12: Loop for K-Fold Cross-Validation (Modified for MedAE)
for k = 1:K
    % Split data into training and test sets for the current fold
    trainData = T_agg(training(cv, k), :);
    testData = T_agg(test(cv, k), :);
    
    % Fit the GLMM model on training data
    glme = fitglme(T_agg, ...
    'nSFR ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');
    
    % Predict on the test data
    cv_predictions = predict(glme, testData);
    
    % Actual values (test data)
    cv_actual = testData.nSFR;
    
    % Calculate absolute errors
    cv_error = abs(cv_predictions - cv_actual);  % Absolute errors
    
    % Calculate MedAE for this fold
    medae_kfold(k) = median(cv_error);      % Median Absolute Error
end

% Calculate the average MedAE across all folds
avg_medae = mean(medae_kfold);

% Display the results
fprintf('Median Absolute Error (MedAE) (K-Fold CV): %.3f\n', avg_medae);

%% Step 12: GLMM post-hoc pairwise contrasts (nSFR) -- robust version
fprintf('\n=== Step 12: GLMM post-hoc pairwise contrasts (nSFR) ===\n');

beta_hat   = fixedEffects(glme);
Sigma_beta = glme.CoefficientCovariance;
coefNames  = glme.CoefficientNames;

divLevels = categories(T_agg.DIVISION);
sexLevels = categories(T_agg.Sex);
expLevels = categories(T_agg.Exposure);

% Wald function: take a contrast vector L, return estimate and SE
wald1 = @(L) deal(L * beta_hat, sqrt(max(L * Sigma_beta * L', 0)));

rows = {};
idx  = 1;

%% (A) Sex contrast: Female vs Male within each Exposure × DIVISION
for d = 1:numel(divLevels)
    for e = 1:numel(expLevels)
        div = char(divLevels{d});
        exp = char(expLevels{e});

        Xf = makeXrow_divSexExp(coefNames, div, 'Female', exp);
        Xm = makeXrow_divSexExp(coefNames, div, 'Male',   exp);

        if isempty(Xf) || isempty(Xm)
            continue;
        end

        L = Xf - Xm;     % Female − Male

        [est, se] = wald1(L);
        if se == 0 || isnan(se)
            continue;
        end

        z     = est / se;
        p     = 2 * normcdf(-abs(z));
        ci_lo = est - 1.96 * se;
        ci_hi = est + 1.96 * se;

        fprintf('DIV=%s | Exposure=%s | Female−Male: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
            div, exp, est, se, z, p);

        rows(idx,:) = { ...
            div, 'Sex_Female_vs_Male', ...   % ContrastLabel
            'Sex', ...                       % ContrastType
            exp, ...                         % ConditioningLevel (Exposure)
            'Female','Male', ...             % Level1, Level2
            est,se,z,p,ci_lo,ci_hi};
        idx = idx + 1;
    end
end

%% (B) Exposure contrasts within each Sex × DIVISION
pairsExp = {'Control','GRIN2B'; 'Control','VPA'; 'GRIN2B','VPA'};

for d = 1:numel(divLevels)
    for s = 1:numel(sexLevels)
        div = char(divLevels{d});
        sex = char(sexLevels{s});

        for k = 1:size(pairsExp,1)
            E1 = pairsExp{k,1};  % numerator
            E2 = pairsExp{k,2};  % denominator

            X1 = makeXrow_divSexExp(coefNames, div, sex, E1);
            X2 = makeXrow_divSexExp(coefNames, div, sex, E2);

            if isempty(X1) || isempty(X2)
                continue;
            end

            L = X1 - X2;    % E1 − E2

            [est, se] = wald1(L);
            if se == 0 || isnan(se)
                continue;
            end

            z     = est / se;
            p     = 2 * normcdf(-abs(z));
            ci_lo = est - 1.96 * se;
            ci_hi = est + 1.96 * se;

            fprintf('DIV=%s | Sex=%s | %s−%s: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
                div, sex, E1, E2, est, se, z, p);

            rows(idx,:) = { ...
                div, sprintf('Exposure_%s_minus_%s',E1,E2), ... % ContrastLabel
                'Exposure', ...         % ContrastType
                sex, ...                % ConditioningLevel (Sex)
                E1, E2, ...             % Level1, Level2
                est,se,z,p,ci_lo,ci_hi};
            idx = idx + 1;
        end
    end
end

% Build summary table
ContrastTable = cell2table(rows, 'VariableNames', { ...
    'DIVISION', ...
    'ContrastLabel', ...      % e.g. 'Sex_Female_vs_Male' or 'Exposure_Control_minus_VPA'
    'ContrastType', ...       % 'Sex' or 'Exposure'
    'ConditioningLevel', ...  % Exposure (for Sex contrasts) or Sex (for Exposure contrasts)
    'Level1', 'Level2', ...   % Level1 − Level2
    'Estimate', 'SE', 'zValue', 'pValue', ...
    'CI_Lower', 'CI_Upper'});

disp('=== GLMM-based contrasts (nSFR) ===');
disp(ContrastTable);


%% Step 14: Visualization of GLMM input estimates using violinplotwhole
% Here we visualize the subject-level median nSFR (T_agg),
% grouped by Sex × Exposure (FaC, MaC, FaV, MaV, FaG, MaG),
% separately for CNIC and RCIC, using your violinplotwhole function.

% 1) Build a plotting table in the format expected by violinplotwhole:
%    variables needed: DIVISION, Level, Group, iDE (here, iDE = nSFR)
T_plot = T_agg(:, {'SubjectID','DIVISION','Sex','Exposure','nSFR'});

% Level: all are Adult / High by construction of tbl
T_plot.Level = repmat("High", height(T_plot), 1);   % string array is fine

% Group code: FaC, MaC, FaV, MaV, FaG, MaG
T_plot.Group = strings(height(T_plot),1);
for i = 1:height(T_plot)
    sex = char(T_plot.Sex(i));
    exp = char(T_plot.Exposure(i));

    if strcmp(sex,'Female') && strcmp(exp,'Control')
        T_plot.Group(i) = 'FaC';
    elseif strcmp(sex,'Male') && strcmp(exp,'Control')
        T_plot.Group(i) = 'MaC';
    elseif strcmp(sex,'Female') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'FaV';
    elseif strcmp(sex,'Male') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'MaV';
    elseif strcmp(sex,'Female') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'FaG';
    elseif strcmp(sex,'Male') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'MaG';
    else
        T_plot.Group(i) = "NA";
    end
end

T_plot.Group = categorical(T_plot.Group);

% The function expects the y-variable to be named iDE → alias nSFR as iDE
T_plot.iDE = T_plot.nSFR;

% Make sure DIVISION is categorical (should already be, but just in case)
T_plot.DIVISION = categorical(T_plot.DIVISION);

% 2) Call your violinplotwhole for each division at High level
levels = "High";   % we only have High in this analysis

% Lemniscal division (CNIC)
if any(T_plot.DIVISION == "CNIC")
    violinplotwhole(T_plot, "CNIC", levels, glme);
end

% Non-lemniscal division (RCIC)
if any(T_plot.DIVISION == "RCIC")
    violinplotwhole(T_plot, "RCIC", levels, glme);
end

%% Now, generalized linear mixed effects model to assess contextual elicited activity and predictive components

%% %% ------------- GLMM for NOR_DEV (Normal distribution) ----------- %%

%% Step 1: Load Data (Assuming already done in previous steps)
data = load('UpdatedFinalTable.mat');
tbl = data.UpdatedFinalTable;

%% Step 2: Convert Categorical Variables
tbl.Sex = categorical(tbl.Sex);              
tbl.Exposure = categorical(tbl.Exposure);    
tbl.DIVISION = categorical(tbl.DIVISION);    

if iscell(tbl.SubjectID)
    tbl.SubjectID = categorical(tbl.SubjectID);
end

%% Step 2.5: Filter to include only Adult animals under High stimulation
% Check if required columns exist
if ~all(ismember({'Age', 'Level'}, tbl.Properties.VariableNames))
    error('Required columns "Age" and "Level" not found in the table.');
end

tbl = tbl(strcmp(tbl.Age, 'Adult') & strcmp(tbl.Level, 'High'), :);


%% Step 3: Aggregate by Subject × Division
[G, subjectList, divisionList] = findgroups(tbl.SubjectID, tbl.DIVISION);
median_NOR_DEV = splitapply(@median, tbl.NOR_DEV, G);

T_agg = table(subjectList, divisionList, median_NOR_DEV, ...
    'VariableNames', {'SubjectID', 'DIVISION', 'NOR_DEV'});


%% Step 5: Add Subject-level Metadata (Sex, Exposure)
T_subject = unique(tbl(:, {'SubjectID', 'Sex', 'Exposure'}));  % One row per subject
T_agg = outerjoin(T_agg, T_subject, 'Keys', 'SubjectID', 'MergeKeys', true);

%% Step 5: Correlation Check
tbl_numeric = tbl(:, {'DIVISION', 'Sex', 'Exposure'}); 
tbl_numeric.DIVISION = double(tbl_numeric.DIVISION);  % Convert to numeric
tbl_numeric.Sex = double(tbl_numeric.Sex);  % Convert to numeric
tbl_numeric.Exposure = double(tbl_numeric.Exposure);  % Convert to numeric

% Calculate the correlation matrix
corr_matrix = corr(tbl_numeric{:, :});  % Compute correlation between numeric variables
disp('Correlation Matrix:');
disp(corr_matrix);

%% Step 6: Fit GLMM (Normal, Identity Link, Laplace Fit)
glme = fitglme(T_agg, ...
    'NOR_DEV ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');

%% Step 7: Display Model Results
disp(glme);

%% Step 8: Model Fit Statistics (AIC and BIC)
AIC = glme.ModelCriterion.AIC;
BIC = glme.ModelCriterion.BIC;
LogLikelihood = glme.LogLikelihood;

disp(['AIC: ', num2str(AIC)]);
disp(['BIC: ', num2str(BIC)]);
disp(['Log Likelihood: ', num2str(LogLikelihood)]);

%% Step 9: Compute Intraclass Correlation (ICC)
[psi, mse] = covarianceParameters(glme);

% Random intercept variance (SubjectID)
var_intercept = psi{1}(1,1);

% Residual variance
var_residual = mse;

% Compute ICC
ICC = var_intercept / (var_intercept + var_residual);
fprintf('Intraclass Correlation (ICC): %.3f\n', ICC);

%% Step 10: Visualize Residuals
residuals = glme.Residuals.Raw;

% Histogram of Residuals
figure;
histogram(residuals, 30);
xlabel('Residuals');
ylabel('Frequency');
title('Histogram of Residuals');
grid on;

% QQ Plot of Residuals
figure;
qqplot(residuals);
title('QQ Plot of Residuals');
grid on;

%% Step 11: K-Fold Cross-Validation Setup for MedAE
K = 10;  % Specify number of folds for cross-validation
cv = cvpartition(size(T_agg, 1), 'KFold', K);  % Create K-fold partition object

medae_kfold = zeros(K, 1);  % Store MedAE for each fold

%% Step 12: Loop for K-Fold Cross-Validation (Modified for MedAE)
for k = 1:K
    % Split data into training and test sets for the current fold
    trainData = T_agg(training(cv, k), :);
    testData = T_agg(test(cv, k), :);
    
    % Fit the GLMM model on training data
    glme = fitglme(T_agg, ...
    'NOR_DEV ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');
    
    % Predict on the test data
    cv_predictions = predict(glme, testData);
    
    % Actual values (test data)
    cv_actual = testData.NOR_DEV;
    
    % Calculate absolute errors
    cv_error = abs(cv_predictions - cv_actual);  % Absolute errors
    
    % Calculate MedAE for this fold
    medae_kfold(k) = median(cv_error);      % Median Absolute Error
end

% Calculate the average MedAE across all folds
avg_medae = mean(medae_kfold);

% Display the results
fprintf('Median Absolute Error (MedAE) (K-Fold CV): %.3f\n', avg_medae);

%% Step 12: GLMM post-hoc pairwise contrasts (DEV) -- robust version
fprintf('\n=== Step 12: GLMM post-hoc pairwise contrasts (DEV) ===\n');

beta_hat   = fixedEffects(glme);
Sigma_beta = glme.CoefficientCovariance;
coefNames  = glme.CoefficientNames;

divLevels = categories(T_agg.DIVISION);
sexLevels = categories(T_agg.Sex);
expLevels = categories(T_agg.Exposure);

% Wald function: take a contrast vector L, return estimate and SE
wald1 = @(L) deal(L * beta_hat, sqrt(max(L * Sigma_beta * L', 0)));

rows = {};
idx  = 1;

%% (A) Sex contrast: Female vs Male within each Exposure × DIVISION
for d = 1:numel(divLevels)
    for e = 1:numel(expLevels)
        div = char(divLevels{d});
        exp = char(expLevels{e});

        Xf = makeXrow_divSexExp(coefNames, div, 'Female', exp);
        Xm = makeXrow_divSexExp(coefNames, div, 'Male',   exp);

        if isempty(Xf) || isempty(Xm)
            continue;
        end

        L = Xf - Xm;     % Female − Male

        [est, se] = wald1(L);
        if se == 0 || isnan(se)
            continue;
        end

        z     = est / se;
        p     = 2 * normcdf(-abs(z));
        ci_lo = est - 1.96 * se;
        ci_hi = est + 1.96 * se;

        fprintf('DIV=%s | Exposure=%s | Female−Male: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
            div, exp, est, se, z, p);

        rows(idx,:) = { ...
            div, 'Sex_Female_vs_Male', ...   % ContrastLabel
            'Sex', ...                       % ContrastType
            exp, ...                         % ConditioningLevel (Exposure)
            'Female','Male', ...             % Level1, Level2
            est,se,z,p,ci_lo,ci_hi};
        idx = idx + 1;
    end
end

%% (B) Exposure contrasts within each Sex × DIVISION
pairsExp = {'Control','GRIN2B'; 'Control','VPA'; 'GRIN2B','VPA'};

for d = 1:numel(divLevels)
    for s = 1:numel(sexLevels)
        div = char(divLevels{d});
        sex = char(sexLevels{s});

        for k = 1:size(pairsExp,1)
            E1 = pairsExp{k,1};  % numerator
            E2 = pairsExp{k,2};  % denominator

            X1 = makeXrow_divSexExp(coefNames, div, sex, E1);
            X2 = makeXrow_divSexExp(coefNames, div, sex, E2);

            if isempty(X1) || isempty(X2)
                continue;
            end

            L = X1 - X2;    % E1 − E2

            [est, se] = wald1(L);
            if se == 0 || isnan(se)
                continue;
            end

            z     = est / se;
            p     = 2 * normcdf(-abs(z));
            ci_lo = est - 1.96 * se;
            ci_hi = est + 1.96 * se;

            fprintf('DIV=%s | Sex=%s | %s−%s: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
                div, sex, E1, E2, est, se, z, p);

            rows(idx,:) = { ...
                div, sprintf('Exposure_%s_minus_%s',E1,E2), ... % ContrastLabel
                'Exposure', ...         % ContrastType
                sex, ...                % ConditioningLevel (Sex)
                E1, E2, ...             % Level1, Level2
                est,se,z,p,ci_lo,ci_hi};
            idx = idx + 1;
        end
    end
end

% Build summary table
ContrastTable = cell2table(rows, 'VariableNames', { ...
    'DIVISION', ...
    'ContrastLabel', ...      % e.g. 'Sex_Female_vs_Male' or 'Exposure_Control_minus_VPA'
    'ContrastType', ...       % 'Sex' or 'Exposure'
    'ConditioningLevel', ...  % Exposure (for Sex contrasts) or Sex (for Exposure contrasts)
    'Level1', 'Level2', ...   % Level1 − Level2
    'Estimate', 'SE', 'zValue', 'pValue', ...
    'CI_Lower', 'CI_Upper'});

disp('=== GLMM-based contrasts (DEV) ===');
disp(ContrastTable);


%% Step 14: Visualization of GLMM input estimates using violinplotwhole
% Here we visualize the subject-level median nSFR (T_agg),
% grouped by Sex × Exposure (FaC, MaC, FaV, MaV, FaG, MaG),
% separately for CNIC and RCIC, using your violinplotwhole function.

% 1) Build a plotting table in the format expected by violinplotwhole:
%    variables needed: DIVISION, Level, Group, iDE (here, iDE = nSFR)
T_plot = T_agg(:, {'SubjectID','DIVISION','Sex','Exposure','NOR_DEV'});

% Level: all are Adult / High by construction of tbl
T_plot.Level = repmat("High", height(T_plot), 1);   % string array is fine

% Group code: FaC, MaC, FaV, MaV, FaG, MaG
T_plot.Group = strings(height(T_plot),1);
for i = 1:height(T_plot)
    sex = char(T_plot.Sex(i));
    exp = char(T_plot.Exposure(i));

    if strcmp(sex,'Female') && strcmp(exp,'Control')
        T_plot.Group(i) = 'FaC';
    elseif strcmp(sex,'Male') && strcmp(exp,'Control')
        T_plot.Group(i) = 'MaC';
    elseif strcmp(sex,'Female') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'FaV';
    elseif strcmp(sex,'Male') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'MaV';
    elseif strcmp(sex,'Female') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'FaG';
    elseif strcmp(sex,'Male') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'MaG';
    else
        T_plot.Group(i) = "NA";
    end
end

T_plot.Group = categorical(T_plot.Group);

% The function expects the y-variable to be named iDE → alias nSFR as iDE
T_plot.iDE = T_plot.NOR_DEV;

% Make sure DIVISION is categorical (should already be, but just in case)
T_plot.DIVISION = categorical(T_plot.DIVISION);

% 2) Call your violinplotwhole for each division at High level
levels = "High";   % we only have High in this analysis

% Lemniscal division (CNIC)
if any(T_plot.DIVISION == "CNIC")
    violinplotwhole(T_plot, "CNIC", levels, glme);
end

% Non-lemniscal division (RCIC)
if any(T_plot.DIVISION == "RCIC")
    violinplotwhole(T_plot, "RCIC", levels, glme);
end

%% %% ------------- GLMM for NOR_CAS (Normal distribution) ----------- %%

%% Step 1: Load Data (Assuming already done in previous steps)
data = load('UpdatedFinalTable.mat');
tbl = data.UpdatedFinalTable;

%% Step 2: Convert Categorical Variables
tbl.Sex = categorical(tbl.Sex);              
tbl.Exposure = categorical(tbl.Exposure);    
tbl.DIVISION = categorical(tbl.DIVISION);    

if iscell(tbl.SubjectID)
    tbl.SubjectID = categorical(tbl.SubjectID);
end

%% Step 2.5: Filter to include only Adult animals under High stimulation
% Check if required columns exist
if ~all(ismember({'Age', 'Level'}, tbl.Properties.VariableNames))
    error('Required columns "Age" and "Level" not found in the table.');
end

tbl = tbl(strcmp(tbl.Age, 'Adult') & strcmp(tbl.Level, 'High'), :);


%% Step 3: Aggregate by Subject × Division
[G, subjectList, divisionList] = findgroups(tbl.SubjectID, tbl.DIVISION);
median_NOR_CTR = splitapply(@median, tbl.NOR_CTR, G);

T_agg = table(subjectList, divisionList, median_NOR_CTR, ...
    'VariableNames', {'SubjectID', 'DIVISION', 'NOR_CTR'});


%% Step 5: Add Subject-level Metadata (Sex, Exposure)
T_subject = unique(tbl(:, {'SubjectID', 'Sex', 'Exposure'}));  % One row per subject
T_agg = outerjoin(T_agg, T_subject, 'Keys', 'SubjectID', 'MergeKeys', true);

%% Step 5: Correlation Check
tbl_numeric = tbl(:, {'DIVISION', 'Sex', 'Exposure'}); 
tbl_numeric.DIVISION = double(tbl_numeric.DIVISION);  % Convert to numeric
tbl_numeric.Sex = double(tbl_numeric.Sex);  % Convert to numeric
tbl_numeric.Exposure = double(tbl_numeric.Exposure);  % Convert to numeric

% Calculate the correlation matrix
corr_matrix = corr(tbl_numeric{:, :});  % Compute correlation between numeric variables
disp('Correlation Matrix:');
disp(corr_matrix);

%% Step 6: Fit GLMM (Normal, Identity Link, Laplace Fit)
glme = fitglme(T_agg, ...
    'NOR_CTR ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');

%% Step 7: Display Model Results
disp(glme);

%% Step 8: Model Fit Statistics (AIC and BIC)
AIC = glme.ModelCriterion.AIC;
BIC = glme.ModelCriterion.BIC;
LogLikelihood = glme.LogLikelihood;

disp(['AIC: ', num2str(AIC)]);
disp(['BIC: ', num2str(BIC)]);
disp(['Log Likelihood: ', num2str(LogLikelihood)]);

%% Step 9: Compute Intraclass Correlation (ICC)
[psi, mse] = covarianceParameters(glme);

% Random intercept variance (SubjectID)
var_intercept = psi{1}(1,1);

% Residual variance
var_residual = mse;

% Compute ICC
ICC = var_intercept / (var_intercept + var_residual);
fprintf('Intraclass Correlation (ICC): %.3f\n', ICC);

%% Step 10: Visualize Residuals
residuals = glme.Residuals.Raw;

% Histogram of Residuals
figure;
histogram(residuals, 30);
xlabel('Residuals');
ylabel('Frequency');
title('Histogram of Residuals');
grid on;

% QQ Plot of Residuals
figure;
qqplot(residuals);
title('QQ Plot of Residuals');
grid on;

%% Step 11: K-Fold Cross-Validation Setup for MedAE
K = 10;  % Specify number of folds for cross-validation
cv = cvpartition(size(T_agg, 1), 'KFold', K);  % Create K-fold partition object

medae_kfold = zeros(K, 1);  % Store MedAE for each fold

%% Step 12: Loop for K-Fold Cross-Validation (Modified for MedAE)
for k = 1:K
    % Split data into training and test sets for the current fold
    trainData = T_agg(training(cv, k), :);
    testData = T_agg(test(cv, k), :);
    
    % Fit the GLMM model on training data
    glme = fitglme(T_agg, ...
    'NOR_CTR ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');
    
    % Predict on the test data
    cv_predictions = predict(glme, testData);
    
    % Actual values (test data)
    cv_actual = testData.NOR_CTR;
    
    % Calculate absolute errors
    cv_error = abs(cv_predictions - cv_actual);  % Absolute errors
    
    % Calculate MedAE for this fold
    medae_kfold(k) = median(cv_error);      % Median Absolute Error
end

% Calculate the average MedAE across all folds
avg_medae = mean(medae_kfold);

% Display the results
fprintf('Median Absolute Error (MedAE) (K-Fold CV): %.3f\n', avg_medae);

%% Step 12: GLMM post-hoc pairwise contrasts (CAS) -- robust version
fprintf('\n=== Step 12: GLMM post-hoc pairwise contrasts (CAS) ===\n');

beta_hat   = fixedEffects(glme);
Sigma_beta = glme.CoefficientCovariance;
coefNames  = glme.CoefficientNames;

divLevels = categories(T_agg.DIVISION);
sexLevels = categories(T_agg.Sex);
expLevels = categories(T_agg.Exposure);

% Wald function: take a contrast vector L, return estimate and SE
wald1 = @(L) deal(L * beta_hat, sqrt(max(L * Sigma_beta * L', 0)));

rows = {};
idx  = 1;

%% (A) Sex contrast: Female vs Male within each Exposure × DIVISION
for d = 1:numel(divLevels)
    for e = 1:numel(expLevels)
        div = char(divLevels{d});
        exp = char(expLevels{e});

        Xf = makeXrow_divSexExp(coefNames, div, 'Female', exp);
        Xm = makeXrow_divSexExp(coefNames, div, 'Male',   exp);

        if isempty(Xf) || isempty(Xm)
            continue;
        end

        L = Xf - Xm;     % Female − Male

        [est, se] = wald1(L);
        if se == 0 || isnan(se)
            continue;
        end

        z     = est / se;
        p     = 2 * normcdf(-abs(z));
        ci_lo = est - 1.96 * se;
        ci_hi = est + 1.96 * se;

        fprintf('DIV=%s | Exposure=%s | Female−Male: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
            div, exp, est, se, z, p);

        rows(idx,:) = { ...
            div, 'Sex_Female_vs_Male', ...   % ContrastLabel
            'Sex', ...                       % ContrastType
            exp, ...                         % ConditioningLevel (Exposure)
            'Female','Male', ...             % Level1, Level2
            est,se,z,p,ci_lo,ci_hi};
        idx = idx + 1;
    end
end

%% (B) Exposure contrasts within each Sex × DIVISION
pairsExp = {'Control','GRIN2B'; 'Control','VPA'; 'GRIN2B','VPA'};

for d = 1:numel(divLevels)
    for s = 1:numel(sexLevels)
        div = char(divLevels{d});
        sex = char(sexLevels{s});

        for k = 1:size(pairsExp,1)
            E1 = pairsExp{k,1};  % numerator
            E2 = pairsExp{k,2};  % denominator

            X1 = makeXrow_divSexExp(coefNames, div, sex, E1);
            X2 = makeXrow_divSexExp(coefNames, div, sex, E2);

            if isempty(X1) || isempty(X2)
                continue;
            end

            L = X1 - X2;    % E1 − E2

            [est, se] = wald1(L);
            if se == 0 || isnan(se)
                continue;
            end

            z     = est / se;
            p     = 2 * normcdf(-abs(z));
            ci_lo = est - 1.96 * se;
            ci_hi = est + 1.96 * se;

            fprintf('DIV=%s | Sex=%s | %s−%s: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
                div, sex, E1, E2, est, se, z, p);

            rows(idx,:) = { ...
                div, sprintf('Exposure_%s_minus_%s',E1,E2), ... % ContrastLabel
                'Exposure', ...         % ContrastType
                sex, ...                % ConditioningLevel (Sex)
                E1, E2, ...             % Level1, Level2
                est,se,z,p,ci_lo,ci_hi};
            idx = idx + 1;
        end
    end
end

% Build summary table
ContrastTable = cell2table(rows, 'VariableNames', { ...
    'DIVISION', ...
    'ContrastLabel', ...      % e.g. 'Sex_Female_vs_Male' or 'Exposure_Control_minus_VPA'
    'ContrastType', ...       % 'Sex' or 'Exposure'
    'ConditioningLevel', ...  % Exposure (for Sex contrasts) or Sex (for Exposure contrasts)
    'Level1', 'Level2', ...   % Level1 − Level2
    'Estimate', 'SE', 'zValue', 'pValue', ...
    'CI_Lower', 'CI_Upper'});

disp('=== GLMM-based contrasts (CAS) ===');
disp(ContrastTable);


%% Step 14: Visualization of GLMM input estimates using violinplotwhole
% Here we visualize the subject-level median nSFR (T_agg),
% grouped by Sex × Exposure (FaC, MaC, FaV, MaV, FaG, MaG),
% separately for CNIC and RCIC, using your violinplotwhole function.

% 1) Build a plotting table in the format expected by violinplotwhole:
%    variables needed: DIVISION, Level, Group, iDE (here, iDE = nSFR)
T_plot = T_agg(:, {'SubjectID','DIVISION','Sex','Exposure','NOR_CTR'});

% Level: all are Adult / High by construction of tbl
T_plot.Level = repmat("High", height(T_plot), 1);   % string array is fine

% Group code: FaC, MaC, FaV, MaV, FaG, MaG
T_plot.Group = strings(height(T_plot),1);
for i = 1:height(T_plot)
    sex = char(T_plot.Sex(i));
    exp = char(T_plot.Exposure(i));

    if strcmp(sex,'Female') && strcmp(exp,'Control')
        T_plot.Group(i) = 'FaC';
    elseif strcmp(sex,'Male') && strcmp(exp,'Control')
        T_plot.Group(i) = 'MaC';
    elseif strcmp(sex,'Female') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'FaV';
    elseif strcmp(sex,'Male') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'MaV';
    elseif strcmp(sex,'Female') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'FaG';
    elseif strcmp(sex,'Male') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'MaG';
    else
        T_plot.Group(i) = "NA";
    end
end

T_plot.Group = categorical(T_plot.Group);

% The function expects the y-variable to be named iDE → alias nSFR as iDE
T_plot.iDE = T_plot.NOR_CTR;

% Make sure DIVISION is categorical (should already be, but just in case)
T_plot.DIVISION = categorical(T_plot.DIVISION);

% 2) Call your violinplotwhole for each division at High level
levels = "High";   % we only have High in this analysis

% Lemniscal division (CNIC)
if any(T_plot.DIVISION == "CNIC")
    violinplotwhole(T_plot, "CNIC", levels, glme);
end

% Non-lemniscal division (RCIC)
if any(T_plot.DIVISION == "RCIC")
    violinplotwhole(T_plot, "RCIC", levels, glme);
end


%% %% ------------- GLMM for NOR_STD (Normal distribution) ----------- %%

%% Step 1: Load Data (Assuming already done in previous steps)
data = load('UpdatedFinalTable.mat');
tbl = data.UpdatedFinalTable;

%% Step 2: Convert Categorical Variables
tbl.Sex = categorical(tbl.Sex);              
tbl.Exposure = categorical(tbl.Exposure);    
tbl.DIVISION = categorical(tbl.DIVISION);    

if iscell(tbl.SubjectID)
    tbl.SubjectID = categorical(tbl.SubjectID);
end

%% Step 2.5: Filter to include only Adult animals under High stimulation
% Check if required columns exist
if ~all(ismember({'Age', 'Level'}, tbl.Properties.VariableNames))
    error('Required columns "Age" and "Level" not found in the table.');
end

tbl = tbl(strcmp(tbl.Age, 'Adult') & strcmp(tbl.Level, 'High'), :);


%% Step 3: Aggregate by Subject × Division
[G, subjectList, divisionList] = findgroups(tbl.SubjectID, tbl.DIVISION);
median_NOR_STD = splitapply(@median, tbl.NOR_STD, G);

T_agg = table(subjectList, divisionList, median_NOR_STD, ...
    'VariableNames', {'SubjectID', 'DIVISION', 'NOR_STD'});


%% Step 5: Add Subject-level Metadata (Sex, Exposure)
T_subject = unique(tbl(:, {'SubjectID', 'Sex', 'Exposure'}));  % One row per subject
T_agg = outerjoin(T_agg, T_subject, 'Keys', 'SubjectID', 'MergeKeys', true);

%% Step 5: Correlation Check
tbl_numeric = tbl(:, {'DIVISION', 'Sex', 'Exposure'}); 
tbl_numeric.DIVISION = double(tbl_numeric.DIVISION);  % Convert to numeric
tbl_numeric.Sex = double(tbl_numeric.Sex);  % Convert to numeric
tbl_numeric.Exposure = double(tbl_numeric.Exposure);  % Convert to numeric

% Calculate the correlation matrix
corr_matrix = corr(tbl_numeric{:, :});  % Compute correlation between numeric variables
disp('Correlation Matrix:');
disp(corr_matrix);

%% Step 6: Fit GLMM (Normal, Identity Link, Laplace Fit)
glme = fitglme(T_agg, ...
    'NOR_STD ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');

%% Step 7: Display Model Results
disp(glme);

%% Step 8: Model Fit Statistics (AIC and BIC)
AIC = glme.ModelCriterion.AIC;
BIC = glme.ModelCriterion.BIC;
LogLikelihood = glme.LogLikelihood;

disp(['AIC: ', num2str(AIC)]);
disp(['BIC: ', num2str(BIC)]);
disp(['Log Likelihood: ', num2str(LogLikelihood)]);

%% Step 9: Compute Intraclass Correlation (ICC)
[psi, mse] = covarianceParameters(glme);

% Random intercept variance (SubjectID)
var_intercept = psi{1}(1,1);

% Residual variance
var_residual = mse;

% Compute ICC
ICC = var_intercept / (var_intercept + var_residual);
fprintf('Intraclass Correlation (ICC): %.3f\n', ICC);

%% Step 10: Visualize Residuals
residuals = glme.Residuals.Raw;

% Histogram of Residuals
figure;
histogram(residuals, 30);
xlabel('Residuals');
ylabel('Frequency');
title('Histogram of Residuals');
grid on;

% QQ Plot of Residuals
figure;
qqplot(residuals);
title('QQ Plot of Residuals');
grid on;

%% Step 11: K-Fold Cross-Validation Setup for MedAE
K = 10;  % Specify number of folds for cross-validation
cv = cvpartition(size(T_agg, 1), 'KFold', K);  % Create K-fold partition object

medae_kfold = zeros(K, 1);  % Store MedAE for each fold

%% Step 12: Loop for K-Fold Cross-Validation (Modified for MedAE)
for k = 1:K
    % Split data into training and test sets for the current fold
    trainData = T_agg(training(cv, k), :);
    testData = T_agg(test(cv, k), :);
    
    % Fit the GLMM model on training data
    glme = fitglme(T_agg, ...
    'NOR_STD ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');
    
    % Predict on the test data
    cv_predictions = predict(glme, testData);
    
    % Actual values (test data)
    cv_actual = testData.NOR_STD;
    
    % Calculate absolute errors
    cv_error = abs(cv_predictions - cv_actual);  % Absolute errors
    
    % Calculate MedAE for this fold
    medae_kfold(k) = median(cv_error);      % Median Absolute Error
end

% Calculate the average MedAE across all folds
avg_medae = mean(medae_kfold);

% Display the results
fprintf('Median Absolute Error (MedAE) (K-Fold CV): %.3f\n', avg_medae);

%% Step 12: GLMM post-hoc pairwise contrasts (STD) -- robust version
fprintf('\n=== Step 12: GLMM post-hoc pairwise contrasts (STD) ===\n');

beta_hat   = fixedEffects(glme);
Sigma_beta = glme.CoefficientCovariance;
coefNames  = glme.CoefficientNames;

divLevels = categories(T_agg.DIVISION);
sexLevels = categories(T_agg.Sex);
expLevels = categories(T_agg.Exposure);

% Wald function: take a contrast vector L, return estimate and SE
wald1 = @(L) deal(L * beta_hat, sqrt(max(L * Sigma_beta * L', 0)));

rows = {};
idx  = 1;

%% (A) Sex contrast: Female vs Male within each Exposure × DIVISION
for d = 1:numel(divLevels)
    for e = 1:numel(expLevels)
        div = char(divLevels{d});
        exp = char(expLevels{e});

        Xf = makeXrow_divSexExp(coefNames, div, 'Female', exp);
        Xm = makeXrow_divSexExp(coefNames, div, 'Male',   exp);

        if isempty(Xf) || isempty(Xm)
            continue;
        end

        L = Xf - Xm;     % Female − Male

        [est, se] = wald1(L);
        if se == 0 || isnan(se)
            continue;
        end

        z     = est / se;
        p     = 2 * normcdf(-abs(z));
        ci_lo = est - 1.96 * se;
        ci_hi = est + 1.96 * se;

        fprintf('DIV=%s | Exposure=%s | Female−Male: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
            div, exp, est, se, z, p);

        rows(idx,:) = { ...
            div, 'Sex_Female_vs_Male', ...   % ContrastLabel
            'Sex', ...                       % ContrastType
            exp, ...                         % ConditioningLevel (Exposure)
            'Female','Male', ...             % Level1, Level2
            est,se,z,p,ci_lo,ci_hi};
        idx = idx + 1;
    end
end

%% (B) Exposure contrasts within each Sex × DIVISION
pairsExp = {'Control','GRIN2B'; 'Control','VPA'; 'GRIN2B','VPA'};

for d = 1:numel(divLevels)
    for s = 1:numel(sexLevels)
        div = char(divLevels{d});
        sex = char(sexLevels{s});

        for k = 1:size(pairsExp,1)
            E1 = pairsExp{k,1};  % numerator
            E2 = pairsExp{k,2};  % denominator

            X1 = makeXrow_divSexExp(coefNames, div, sex, E1);
            X2 = makeXrow_divSexExp(coefNames, div, sex, E2);

            if isempty(X1) || isempty(X2)
                continue;
            end

            L = X1 - X2;    % E1 − E2

            [est, se] = wald1(L);
            if se == 0 || isnan(se)
                continue;
            end

            z     = est / se;
            p     = 2 * normcdf(-abs(z));
            ci_lo = est - 1.96 * se;
            ci_hi = est + 1.96 * se;

            fprintf('DIV=%s | Sex=%s | %s−%s: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
                div, sex, E1, E2, est, se, z, p);

            rows(idx,:) = { ...
                div, sprintf('Exposure_%s_minus_%s',E1,E2), ... % ContrastLabel
                'Exposure', ...         % ContrastType
                sex, ...                % ConditioningLevel (Sex)
                E1, E2, ...             % Level1, Level2
                est,se,z,p,ci_lo,ci_hi};
            idx = idx + 1;
        end
    end
end

% Build summary table
ContrastTable = cell2table(rows, 'VariableNames', { ...
    'DIVISION', ...
    'ContrastLabel', ...      % e.g. 'Sex_Female_vs_Male' or 'Exposure_Control_minus_VPA'
    'ContrastType', ...       % 'Sex' or 'Exposure'
    'ConditioningLevel', ...  % Exposure (for Sex contrasts) or Sex (for Exposure contrasts)
    'Level1', 'Level2', ...   % Level1 − Level2
    'Estimate', 'SE', 'zValue', 'pValue', ...
    'CI_Lower', 'CI_Upper'});

disp('=== GLMM-based contrasts (STD) ===');
disp(ContrastTable);


%% Step 14: Visualization of GLMM input estimates using violinplotwhole
% Here we visualize the subject-level median nSFR (T_agg),
% grouped by Sex × Exposure (FaC, MaC, FaV, MaV, FaG, MaG),
% separately for CNIC and RCIC, using your violinplotwhole function.

% 1) Build a plotting table in the format expected by violinplotwhole:
%    variables needed: DIVISION, Level, Group, iDE (here, iDE = nSFR)
T_plot = T_agg(:, {'SubjectID','DIVISION','Sex','Exposure','NOR_STD'});

% Level: all are Adult / High by construction of tbl
T_plot.Level = repmat("High", height(T_plot), 1);   % string array is fine

% Group code: FaC, MaC, FaV, MaV, FaG, MaG
T_plot.Group = strings(height(T_plot),1);
for i = 1:height(T_plot)
    sex = char(T_plot.Sex(i));
    exp = char(T_plot.Exposure(i));

    if strcmp(sex,'Female') && strcmp(exp,'Control')
        T_plot.Group(i) = 'FaC';
    elseif strcmp(sex,'Male') && strcmp(exp,'Control')
        T_plot.Group(i) = 'MaC';
    elseif strcmp(sex,'Female') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'FaV';
    elseif strcmp(sex,'Male') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'MaV';
    elseif strcmp(sex,'Female') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'FaG';
    elseif strcmp(sex,'Male') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'MaG';
    else
        T_plot.Group(i) = "NA";
    end
end

T_plot.Group = categorical(T_plot.Group);

% The function expects the y-variable to be named iDE → alias nSFR as iDE
T_plot.iDE = T_plot.NOR_STD;

% Make sure DIVISION is categorical (should already be, but just in case)
T_plot.DIVISION = categorical(T_plot.DIVISION);

% 2) Call your violinplotwhole for each division at High level
levels = "High";   % we only have High in this analysis

% Lemniscal division (CNIC)
if any(T_plot.DIVISION == "CNIC")
    violinplotwhole(T_plot, "CNIC", levels, glme);
end

% Non-lemniscal division (RCIC)
if any(T_plot.DIVISION == "RCIC")
    violinplotwhole(T_plot, "RCIC", levels, glme);
end


%% ------------ Predcitive components of mismatch processing ------ %%
% With GLMMs as they work fine and we use the same Distribution lapalce
% with quasinewton optimizer

%% %% ------------- GLMM for iMM (Normal distribution) ----------- %%

%% Step 1: Load Data (Assuming already done in previous steps)
data = load('UpdatedFinalTable.mat');
tbl = data.UpdatedFinalTable;

%% Step 2: Convert Categorical Variables
tbl.Sex = categorical(tbl.Sex);              
tbl.Exposure = categorical(tbl.Exposure);    
tbl.DIVISION = categorical(tbl.DIVISION);    

if iscell(tbl.SubjectID)
    tbl.SubjectID = categorical(tbl.SubjectID);
end

%% Step 2.5: Filter to include only Adult animals under High stimulation
% Check if required columns exist
if ~all(ismember({'Age', 'Level'}, tbl.Properties.VariableNames))
    error('Required columns "Age" and "Level" not found in the table.');
end

tbl = tbl(strcmp(tbl.Age, 'Adult') & strcmp(tbl.Level, 'High'), :);


%% Step 3: Aggregate by Subject × Division
[G, subjectList, divisionList] = findgroups(tbl.SubjectID, tbl.DIVISION);
median_iSSA = splitapply(@median, tbl.iSSA, G);

T_agg = table(subjectList, divisionList, median_iSSA, ...
    'VariableNames', {'SubjectID', 'DIVISION', 'iSSA'});


%% Step 5: Add Subject-level Metadata (Sex, Exposure)
T_subject = unique(tbl(:, {'SubjectID', 'Sex', 'Exposure'}));  % One row per subject
T_agg = outerjoin(T_agg, T_subject, 'Keys', 'SubjectID', 'MergeKeys', true);

%% Step 5: Correlation Check
tbl_numeric = tbl(:, {'DIVISION', 'Sex', 'Exposure'}); 
tbl_numeric.DIVISION = double(tbl_numeric.DIVISION);  % Convert to numeric
tbl_numeric.Sex = double(tbl_numeric.Sex);  % Convert to numeric
tbl_numeric.Exposure = double(tbl_numeric.Exposure);  % Convert to numeric

% Calculate the correlation matrix
corr_matrix = corr(tbl_numeric{:, :});  % Compute correlation between numeric variables
disp('Correlation Matrix:');
disp(corr_matrix);

%% Step 6: Fit GLMM (Normal, Identity Link, Laplace Fit)
glme = fitglme(T_agg, ...
    'iSSA ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');

%% Step 7: Display Model Results
disp(glme);

%% Step 8: Model Fit Statistics (AIC and BIC)
AIC = glme.ModelCriterion.AIC;
BIC = glme.ModelCriterion.BIC;
LogLikelihood = glme.LogLikelihood;

disp(['AIC: ', num2str(AIC)]);
disp(['BIC: ', num2str(BIC)]);
disp(['Log Likelihood: ', num2str(LogLikelihood)]);

%% Step 9: Compute Intraclass Correlation (ICC)
[psi, mse] = covarianceParameters(glme);

% Random intercept variance (SubjectID)
var_intercept = psi{1}(1,1);

% Residual variance
var_residual = mse;

% Compute ICC
ICC = var_intercept / (var_intercept + var_residual);
fprintf('Intraclass Correlation (ICC): %.3f\n', ICC);

%% Step 10: Visualize Residuals
residuals = glme.Residuals.Raw;

% Histogram of Residuals
figure;
histogram(residuals, 30);
xlabel('Residuals');
ylabel('Frequency');
title('Histogram of Residuals');
grid on;

% QQ Plot of Residuals
figure;
qqplot(residuals);
title('QQ Plot of Residuals');
grid on;

%% Step 11: K-Fold Cross-Validation Setup for MedAE
K = 10;  % Specify number of folds for cross-validation
cv = cvpartition(size(T_agg, 1), 'KFold', K);  % Create K-fold partition object

medae_kfold = zeros(K, 1);  % Store MedAE for each fold

%% Step 12: Loop for K-Fold Cross-Validation (Modified for MedAE)
for k = 1:K
    % Split data into training and test sets for the current fold
    trainData = T_agg(training(cv, k), :);
    testData = T_agg(test(cv, k), :);
    
    % Fit the GLMM model on training data
    glme = fitglme(T_agg, ...
    'iSSA ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');
    
    % Predict on the test data
    cv_predictions = predict(glme, testData);
    
    % Actual values (test data)
    cv_actual = testData.iSSA;
    
    % Calculate absolute errors
    cv_error = abs(cv_predictions - cv_actual);  % Absolute errors
    
    % Calculate MedAE for this fold
    medae_kfold(k) = median(cv_error);      % Median Absolute Error
end

% Calculate the average MedAE across all folds
avg_medae = mean(medae_kfold);

% Display the results
fprintf('Median Absolute Error (MedAE) (K-Fold CV): %.3f\n', avg_medae);

%% Step 12: GLMM post-hoc pairwise contrasts (iMM) -- robust version
fprintf('\n=== Step 12: GLMM post-hoc pairwise contrasts (iMM) ===\n');

beta_hat   = fixedEffects(glme);
Sigma_beta = glme.CoefficientCovariance;
coefNames  = glme.CoefficientNames;

divLevels = categories(T_agg.DIVISION);
sexLevels = categories(T_agg.Sex);
expLevels = categories(T_agg.Exposure);

% Wald function: take a contrast vector L, return estimate and SE
wald1 = @(L) deal(L * beta_hat, sqrt(max(L * Sigma_beta * L', 0)));

rows = {};
idx  = 1;

%% (A) Sex contrast: Female vs Male within each Exposure × DIVISION
for d = 1:numel(divLevels)
    for e = 1:numel(expLevels)
        div = char(divLevels{d});
        exp = char(expLevels{e});

        Xf = makeXrow_divSexExp(coefNames, div, 'Female', exp);
        Xm = makeXrow_divSexExp(coefNames, div, 'Male',   exp);

        if isempty(Xf) || isempty(Xm)
            continue;
        end

        L = Xf - Xm;     % Female − Male

        [est, se] = wald1(L);
        if se == 0 || isnan(se)
            continue;
        end

        z     = est / se;
        p     = 2 * normcdf(-abs(z));
        ci_lo = est - 1.96 * se;
        ci_hi = est + 1.96 * se;

        fprintf('DIV=%s | Exposure=%s | Female−Male: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
            div, exp, est, se, z, p);

        rows(idx,:) = { ...
            div, 'Sex_Female_vs_Male', ...   % ContrastLabel
            'Sex', ...                       % ContrastType
            exp, ...                         % ConditioningLevel (Exposure)
            'Female','Male', ...             % Level1, Level2
            est,se,z,p,ci_lo,ci_hi};
        idx = idx + 1;
    end
end

%% (B) Exposure contrasts within each Sex × DIVISION
pairsExp = {'Control','GRIN2B'; 'Control','VPA'; 'GRIN2B','VPA'};

for d = 1:numel(divLevels)
    for s = 1:numel(sexLevels)
        div = char(divLevels{d});
        sex = char(sexLevels{s});

        for k = 1:size(pairsExp,1)
            E1 = pairsExp{k,1};  % numerator
            E2 = pairsExp{k,2};  % denominator

            X1 = makeXrow_divSexExp(coefNames, div, sex, E1);
            X2 = makeXrow_divSexExp(coefNames, div, sex, E2);

            if isempty(X1) || isempty(X2)
                continue;
            end

            L = X1 - X2;    % E1 − E2

            [est, se] = wald1(L);
            if se == 0 || isnan(se)
                continue;
            end

            z     = est / se;
            p     = 2 * normcdf(-abs(z));
            ci_lo = est - 1.96 * se;
            ci_hi = est + 1.96 * se;

            fprintf('DIV=%s | Sex=%s | %s−%s: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
                div, sex, E1, E2, est, se, z, p);

            rows(idx,:) = { ...
                div, sprintf('Exposure_%s_minus_%s',E1,E2), ... % ContrastLabel
                'Exposure', ...         % ContrastType
                sex, ...                % ConditioningLevel (Sex)
                E1, E2, ...             % Level1, Level2
                est,se,z,p,ci_lo,ci_hi};
            idx = idx + 1;
        end
    end
end

% Build summary table
ContrastTable = cell2table(rows, 'VariableNames', { ...
    'DIVISION', ...
    'ContrastLabel', ...      % e.g. 'Sex_Female_vs_Male' or 'Exposure_Control_minus_VPA'
    'ContrastType', ...       % 'Sex' or 'Exposure'
    'ConditioningLevel', ...  % Exposure (for Sex contrasts) or Sex (for Exposure contrasts)
    'Level1', 'Level2', ...   % Level1 − Level2
    'Estimate', 'SE', 'zValue', 'pValue', ...
    'CI_Lower', 'CI_Upper'});

disp('=== GLMM-based contrasts (iMM) ===');
disp(ContrastTable);


%% Step 14: Visualization of GLMM input estimates using violinplotwhole
% Here we visualize the subject-level median nSFR (T_agg),
% grouped by Sex × Exposure (FaC, MaC, FaV, MaV, FaG, MaG),
% separately for CNIC and RCIC, using your violinplotwhole function.

% 1) Build a plotting table in the format expected by violinplotwhole:
%    variables needed: DIVISION, Level, Group, iDE (here, iDE = nSFR)
T_plot = T_agg(:, {'SubjectID','DIVISION','Sex','Exposure','iSSA'});

% Level: all are Adult / High by construction of tbl
T_plot.Level = repmat("High", height(T_plot), 1);   % string array is fine

% Group code: FaC, MaC, FaV, MaV, FaG, MaG
T_plot.Group = strings(height(T_plot),1);
for i = 1:height(T_plot)
    sex = char(T_plot.Sex(i));
    exp = char(T_plot.Exposure(i));

    if strcmp(sex,'Female') && strcmp(exp,'Control')
        T_plot.Group(i) = 'FaC';
    elseif strcmp(sex,'Male') && strcmp(exp,'Control')
        T_plot.Group(i) = 'MaC';
    elseif strcmp(sex,'Female') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'FaV';
    elseif strcmp(sex,'Male') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'MaV';
    elseif strcmp(sex,'Female') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'FaG';
    elseif strcmp(sex,'Male') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'MaG';
    else
        T_plot.Group(i) = "NA";
    end
end

T_plot.Group = categorical(T_plot.Group);

% The function expects the y-variable to be named iDE → alias nSFR as iDE
T_plot.iDE = T_plot.iSSA;

% Make sure DIVISION is categorical (should already be, but just in case)
T_plot.DIVISION = categorical(T_plot.DIVISION);

% 2) Call your violinplotwhole for each division at High level
levels = "High";   % we only have High in this analysis

% Lemniscal division (CNIC)
if any(T_plot.DIVISION == "CNIC")
    violinplotwhole(T_plot, "CNIC", levels, glme);
end

% Non-lemniscal division (RCIC)
if any(T_plot.DIVISION == "RCIC")
    violinplotwhole(T_plot, "RCIC", levels, glme);
end

%% %% ------------- GLMM for iRS (Normal distribution) ----------- %%

%% Step 1: Load Data (Assuming already done in previous steps)
data = load('UpdatedFinalTable.mat');
tbl = data.UpdatedFinalTable;

%% Step 2: Convert Categorical Variables
tbl.Sex = categorical(tbl.Sex);              
tbl.Exposure = categorical(tbl.Exposure);    
tbl.DIVISION = categorical(tbl.DIVISION);    

if iscell(tbl.SubjectID)
    tbl.SubjectID = categorical(tbl.SubjectID);
end

%% Step 2.5: Filter to include only Adult animals under High stimulation
% Check if required columns exist
if ~all(ismember({'Age', 'Level'}, tbl.Properties.VariableNames))
    error('Required columns "Age" and "Level" not found in the table.');
end

tbl = tbl(strcmp(tbl.Age, 'Adult') & strcmp(tbl.Level, 'High'), :);


%% Step 3: Aggregate by Subject × Division
[G, subjectList, divisionList] = findgroups(tbl.SubjectID, tbl.DIVISION);
median_iRA = splitapply(@median, tbl.iRA, G);

T_agg = table(subjectList, divisionList, median_iRA, ...
    'VariableNames', {'SubjectID', 'DIVISION', 'iRA'});


%% Step 5: Add Subject-level Metadata (Sex, Exposure)
T_subject = unique(tbl(:, {'SubjectID', 'Sex', 'Exposure'}));  % One row per subject
T_agg = outerjoin(T_agg, T_subject, 'Keys', 'SubjectID', 'MergeKeys', true);

%% Step 5: Correlation Check
tbl_numeric = tbl(:, {'DIVISION', 'Sex', 'Exposure'}); 
tbl_numeric.DIVISION = double(tbl_numeric.DIVISION);  % Convert to numeric
tbl_numeric.Sex = double(tbl_numeric.Sex);  % Convert to numeric
tbl_numeric.Exposure = double(tbl_numeric.Exposure);  % Convert to numeric

% Calculate the correlation matrix
corr_matrix = corr(tbl_numeric{:, :});  % Compute correlation between numeric variables
disp('Correlation Matrix:');
disp(corr_matrix);

%% Step 6: Fit GLMM (Normal, Identity Link, Laplace Fit)
glme = fitglme(T_agg, ...
    'iRA ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');

%% Step 7: Display Model Results
disp(glme);

%% Step 8: Model Fit Statistics (AIC and BIC)
AIC = glme.ModelCriterion.AIC;
BIC = glme.ModelCriterion.BIC;
LogLikelihood = glme.LogLikelihood;

disp(['AIC: ', num2str(AIC)]);
disp(['BIC: ', num2str(BIC)]);
disp(['Log Likelihood: ', num2str(LogLikelihood)]);

%% Step 9: Compute Intraclass Correlation (ICC)
[psi, mse] = covarianceParameters(glme);

% Random intercept variance (SubjectID)
var_intercept = psi{1}(1,1);

% Residual variance
var_residual = mse;

% Compute ICC
ICC = var_intercept / (var_intercept + var_residual);
fprintf('Intraclass Correlation (ICC): %.3f\n', ICC);

%% Step 10: Visualize Residuals
residuals = glme.Residuals.Raw;

% Histogram of Residuals
figure;
histogram(residuals, 30);
xlabel('Residuals');
ylabel('Frequency');
title('Histogram of Residuals');
grid on;

% QQ Plot of Residuals
figure;
qqplot(residuals);
title('QQ Plot of Residuals');
grid on;

%% Step 11: K-Fold Cross-Validation Setup for MedAE
K = 10;  % Specify number of folds for cross-validation
cv = cvpartition(size(T_agg, 1), 'KFold', K);  % Create K-fold partition object

medae_kfold = zeros(K, 1);  % Store MedAE for each fold

%% Step 12: Loop for K-Fold Cross-Validation (Modified for MedAE)
for k = 1:K
    % Split data into training and test sets for the current fold
    trainData = T_agg(training(cv, k), :);
    testData = T_agg(test(cv, k), :);
    
    % Fit the GLMM model on training data
    glme = fitglme(T_agg, ...
    'iRA ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');
    
    % Predict on the test data
    cv_predictions = predict(glme, testData);
    
    % Actual values (test data)
    cv_actual = testData.iRA;
    
    % Calculate absolute errors
    cv_error = abs(cv_predictions - cv_actual);  % Absolute errors
    
    % Calculate MedAE for this fold
    medae_kfold(k) = median(cv_error);      % Median Absolute Error
end

% Calculate the average MedAE across all folds
avg_medae = mean(medae_kfold);

% Display the results
fprintf('Median Absolute Error (MedAE) (K-Fold CV): %.3f\n', avg_medae);

%% Step 12: GLMM post-hoc pairwise contrasts (iRS) -- robust version
fprintf('\n=== Step 12: GLMM post-hoc pairwise contrasts (iRS) ===\n');

beta_hat   = fixedEffects(glme);
Sigma_beta = glme.CoefficientCovariance;
coefNames  = glme.CoefficientNames;

divLevels = categories(T_agg.DIVISION);
sexLevels = categories(T_agg.Sex);
expLevels = categories(T_agg.Exposure);

% Wald function: take a contrast vector L, return estimate and SE
wald1 = @(L) deal(L * beta_hat, sqrt(max(L * Sigma_beta * L', 0)));

rows = {};
idx  = 1;

%% (A) Sex contrast: Female vs Male within each Exposure × DIVISION
for d = 1:numel(divLevels)
    for e = 1:numel(expLevels)
        div = char(divLevels{d});
        exp = char(expLevels{e});

        Xf = makeXrow_divSexExp(coefNames, div, 'Female', exp);
        Xm = makeXrow_divSexExp(coefNames, div, 'Male',   exp);

        if isempty(Xf) || isempty(Xm)
            continue;
        end

        L = Xf - Xm;     % Female − Male

        [est, se] = wald1(L);
        if se == 0 || isnan(se)
            continue;
        end

        z     = est / se;
        p     = 2 * normcdf(-abs(z));
        ci_lo = est - 1.96 * se;
        ci_hi = est + 1.96 * se;

        fprintf('DIV=%s | Exposure=%s | Female−Male: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
            div, exp, est, se, z, p);

        rows(idx,:) = { ...
            div, 'Sex_Female_vs_Male', ...   % ContrastLabel
            'Sex', ...                       % ContrastType
            exp, ...                         % ConditioningLevel (Exposure)
            'Female','Male', ...             % Level1, Level2
            est,se,z,p,ci_lo,ci_hi};
        idx = idx + 1;
    end
end

%% (B) Exposure contrasts within each Sex × DIVISION
pairsExp = {'Control','GRIN2B'; 'Control','VPA'; 'GRIN2B','VPA'};

for d = 1:numel(divLevels)
    for s = 1:numel(sexLevels)
        div = char(divLevels{d});
        sex = char(sexLevels{s});

        for k = 1:size(pairsExp,1)
            E1 = pairsExp{k,1};  % numerator
            E2 = pairsExp{k,2};  % denominator

            X1 = makeXrow_divSexExp(coefNames, div, sex, E1);
            X2 = makeXrow_divSexExp(coefNames, div, sex, E2);

            if isempty(X1) || isempty(X2)
                continue;
            end

            L = X1 - X2;    % E1 − E2

            [est, se] = wald1(L);
            if se == 0 || isnan(se)
                continue;
            end

            z     = est / se;
            p     = 2 * normcdf(-abs(z));
            ci_lo = est - 1.96 * se;
            ci_hi = est + 1.96 * se;

            fprintf('DIV=%s | Sex=%s | %s−%s: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
                div, sex, E1, E2, est, se, z, p);

            rows(idx,:) = { ...
                div, sprintf('Exposure_%s_minus_%s',E1,E2), ... % ContrastLabel
                'Exposure', ...         % ContrastType
                sex, ...                % ConditioningLevel (Sex)
                E1, E2, ...             % Level1, Level2
                est,se,z,p,ci_lo,ci_hi};
            idx = idx + 1;
        end
    end
end

% Build summary table
ContrastTable = cell2table(rows, 'VariableNames', { ...
    'DIVISION', ...
    'ContrastLabel', ...      % e.g. 'Sex_Female_vs_Male' or 'Exposure_Control_minus_VPA'
    'ContrastType', ...       % 'Sex' or 'Exposure'
    'ConditioningLevel', ...  % Exposure (for Sex contrasts) or Sex (for Exposure contrasts)
    'Level1', 'Level2', ...   % Level1 − Level2
    'Estimate', 'SE', 'zValue', 'pValue', ...
    'CI_Lower', 'CI_Upper'});

disp('=== GLMM-based contrasts (iRS) ===');
disp(ContrastTable);


%% Step 14: Visualization of GLMM input estimates using violinplotwhole
% Here we visualize the subject-level median nSFR (T_agg),
% grouped by Sex × Exposure (FaC, MaC, FaV, MaV, FaG, MaG),
% separately for CNIC and RCIC, using your violinplotwhole function.

% 1) Build a plotting table in the format expected by violinplotwhole:
%    variables needed: DIVISION, Level, Group, iDE (here, iDE = nSFR)
T_plot = T_agg(:, {'SubjectID','DIVISION','Sex','Exposure','iRA'});

% Level: all are Adult / High by construction of tbl
T_plot.Level = repmat("High", height(T_plot), 1);   % string array is fine

% Group code: FaC, MaC, FaV, MaV, FaG, MaG
T_plot.Group = strings(height(T_plot),1);
for i = 1:height(T_plot)
    sex = char(T_plot.Sex(i));
    exp = char(T_plot.Exposure(i));

    if strcmp(sex,'Female') && strcmp(exp,'Control')
        T_plot.Group(i) = 'FaC';
    elseif strcmp(sex,'Male') && strcmp(exp,'Control')
        T_plot.Group(i) = 'MaC';
    elseif strcmp(sex,'Female') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'FaV';
    elseif strcmp(sex,'Male') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'MaV';
    elseif strcmp(sex,'Female') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'FaG';
    elseif strcmp(sex,'Male') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'MaG';
    else
        T_plot.Group(i) = "NA";
    end
end

T_plot.Group = categorical(T_plot.Group);

% The function expects the y-variable to be named iDE → alias nSFR as iDE
T_plot.iDE = T_plot.iRA;

% Make sure DIVISION is categorical (should already be, but just in case)
T_plot.DIVISION = categorical(T_plot.DIVISION);

% 2) Call your violinplotwhole for each division at High level
levels = "High";   % we only have High in this analysis

% Lemniscal division (CNIC)
if any(T_plot.DIVISION == "CNIC")
    violinplotwhole(T_plot, "CNIC", levels, glme);
end

% Non-lemniscal division (RCIC)
if any(T_plot.DIVISION == "RCIC")
    violinplotwhole(T_plot, "RCIC", levels, glme);
end


%% %% ------------- GLMM for iPE (Normal distribution) ----------- %%

%% Step 1: Load Data (Assuming already done in previous steps)
data = load('UpdatedFinalTable.mat');
tbl = data.UpdatedFinalTable;

%% Step 2: Convert Categorical Variables
tbl.Sex = categorical(tbl.Sex);              
tbl.Exposure = categorical(tbl.Exposure);    
tbl.DIVISION = categorical(tbl.DIVISION);    

if iscell(tbl.SubjectID)
    tbl.SubjectID = categorical(tbl.SubjectID);
end

%% Step 2.5: Filter to include only Adult animals under High stimulation
% Check if required columns exist
if ~all(ismember({'Age', 'Level'}, tbl.Properties.VariableNames))
    error('Required columns "Age" and "Level" not found in the table.');
end

tbl = tbl(strcmp(tbl.Age, 'Adult') & strcmp(tbl.Level, 'High'), :);


%% Step 3: Aggregate by Subject × Division
[G, subjectList, divisionList] = findgroups(tbl.SubjectID, tbl.DIVISION);
median_iDE = splitapply(@median, tbl.iDE, G);

T_agg = table(subjectList, divisionList, median_iDE, ...
    'VariableNames', {'SubjectID', 'DIVISION', 'iDE'});


%% Step 5: Add Subject-level Metadata (Sex, Exposure)
T_subject = unique(tbl(:, {'SubjectID', 'Sex', 'Exposure'}));  % One row per subject
T_agg = outerjoin(T_agg, T_subject, 'Keys', 'SubjectID', 'MergeKeys', true);

%% Step 5: Correlation Check
tbl_numeric = tbl(:, {'DIVISION', 'Sex', 'Exposure'}); 
tbl_numeric.DIVISION = double(tbl_numeric.DIVISION);  % Convert to numeric
tbl_numeric.Sex = double(tbl_numeric.Sex);  % Convert to numeric
tbl_numeric.Exposure = double(tbl_numeric.Exposure);  % Convert to numeric

% Calculate the correlation matrix
corr_matrix = corr(tbl_numeric{:, :});  % Compute correlation between numeric variables
disp('Correlation Matrix:');
disp(corr_matrix);

%% Step 6: Fit GLMM (Normal, Identity Link, Laplace Fit)
glme = fitglme(T_agg, ...
    'iDE ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');

%% Step 7: Display Model Results
disp(glme);

%% Step 8: Model Fit Statistics (AIC and BIC)
AIC = glme.ModelCriterion.AIC;
BIC = glme.ModelCriterion.BIC;
LogLikelihood = glme.LogLikelihood;

disp(['AIC: ', num2str(AIC)]);
disp(['BIC: ', num2str(BIC)]);
disp(['Log Likelihood: ', num2str(LogLikelihood)]);

%% Step 9: Compute Intraclass Correlation (ICC)
[psi, mse] = covarianceParameters(glme);

% Random intercept variance (SubjectID)
var_intercept = psi{1}(1,1);

% Residual variance
var_residual = mse;

% Compute ICC
ICC = var_intercept / (var_intercept + var_residual);
fprintf('Intraclass Correlation (ICC): %.3f\n', ICC);

%% Step 10: Visualize Residuals
residuals = glme.Residuals.Raw;

% Histogram of Residuals
figure;
histogram(residuals, 30);
xlabel('Residuals');
ylabel('Frequency');
title('Histogram of Residuals');
grid on;

% QQ Plot of Residuals
figure;
qqplot(residuals);
title('QQ Plot of Residuals');
grid on;

%% Step 11: K-Fold Cross-Validation Setup for MedAE
K = 10;  % Specify number of folds for cross-validation
cv = cvpartition(size(T_agg, 1), 'KFold', K);  % Create K-fold partition object

medae_kfold = zeros(K, 1);  % Store MedAE for each fold

%% Step 12: Loop for K-Fold Cross-Validation (Modified for MedAE)
for k = 1:K
    % Split data into training and test sets for the current fold
    trainData = T_agg(training(cv, k), :);
    testData = T_agg(test(cv, k), :);
    
    % Fit the GLMM model on training data
    glme = fitglme(T_agg, ...
    'iDE ~ Sex*Exposure*DIVISION + (1 + DIVISION | SubjectID)', ...
    'Distribution', 'Normal', 'Link', 'identity', ...
    'FitMethod', 'Laplace', 'Optimizer', 'quasinewton');
    
    % Predict on the test data
    cv_predictions = predict(glme, testData);
    
    % Actual values (test data)
    cv_actual = testData.iDE;
    
    % Calculate absolute errors
    cv_error = abs(cv_predictions - cv_actual);  % Absolute errors
    
    % Calculate MedAE for this fold
    medae_kfold(k) = median(cv_error);      % Median Absolute Error
end

% Calculate the average MedAE across all folds
avg_medae = mean(medae_kfold);

% Display the results
fprintf('Median Absolute Error (MedAE) (K-Fold CV): %.3f\n', avg_medae);

%% Step 12: GLMM post-hoc pairwise contrasts (iRS) -- robust version
fprintf('\n=== Step 12: GLMM post-hoc pairwise contrasts (iPE) ===\n');

beta_hat   = fixedEffects(glme);
Sigma_beta = glme.CoefficientCovariance;
coefNames  = glme.CoefficientNames;

divLevels = categories(T_agg.DIVISION);
sexLevels = categories(T_agg.Sex);
expLevels = categories(T_agg.Exposure);

% Wald function: take a contrast vector L, return estimate and SE
wald1 = @(L) deal(L * beta_hat, sqrt(max(L * Sigma_beta * L', 0)));

rows = {};
idx  = 1;

%% (A) Sex contrast: Female vs Male within each Exposure × DIVISION
for d = 1:numel(divLevels)
    for e = 1:numel(expLevels)
        div = char(divLevels{d});
        exp = char(expLevels{e});

        Xf = makeXrow_divSexExp(coefNames, div, 'Female', exp);
        Xm = makeXrow_divSexExp(coefNames, div, 'Male',   exp);

        if isempty(Xf) || isempty(Xm)
            continue;
        end

        L = Xf - Xm;     % Female − Male

        [est, se] = wald1(L);
        if se == 0 || isnan(se)
            continue;
        end

        z     = est / se;
        p     = 2 * normcdf(-abs(z));
        ci_lo = est - 1.96 * se;
        ci_hi = est + 1.96 * se;

        fprintf('DIV=%s | Exposure=%s | Female−Male: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
            div, exp, est, se, z, p);

        rows(idx,:) = { ...
            div, 'Sex_Female_vs_Male', ...   % ContrastLabel
            'Sex', ...                       % ContrastType
            exp, ...                         % ConditioningLevel (Exposure)
            'Female','Male', ...             % Level1, Level2
            est,se,z,p,ci_lo,ci_hi};
        idx = idx + 1;
    end
end

%% (B) Exposure contrasts within each Sex × DIVISION
pairsExp = {'Control','GRIN2B'; 'Control','VPA'; 'GRIN2B','VPA'};

for d = 1:numel(divLevels)
    for s = 1:numel(sexLevels)
        div = char(divLevels{d});
        sex = char(sexLevels{s});

        for k = 1:size(pairsExp,1)
            E1 = pairsExp{k,1};  % numerator
            E2 = pairsExp{k,2};  % denominator

            X1 = makeXrow_divSexExp(coefNames, div, sex, E1);
            X2 = makeXrow_divSexExp(coefNames, div, sex, E2);

            if isempty(X1) || isempty(X2)
                continue;
            end

            L = X1 - X2;    % E1 − E2

            [est, se] = wald1(L);
            if se == 0 || isnan(se)
                continue;
            end

            z     = est / se;
            p     = 2 * normcdf(-abs(z));
            ci_lo = est - 1.96 * se;
            ci_hi = est + 1.96 * se;

            fprintf('DIV=%s | Sex=%s | %s−%s: est=%.3f, SE=%.3f, z=%.2f, p=%.4f\n', ...
                div, sex, E1, E2, est, se, z, p);

            rows(idx,:) = { ...
                div, sprintf('Exposure_%s_minus_%s',E1,E2), ... % ContrastLabel
                'Exposure', ...         % ContrastType
                sex, ...                % ConditioningLevel (Sex)
                E1, E2, ...             % Level1, Level2
                est,se,z,p,ci_lo,ci_hi};
            idx = idx + 1;
        end
    end
end

% Build summary table
ContrastTable = cell2table(rows, 'VariableNames', { ...
    'DIVISION', ...
    'ContrastLabel', ...      % e.g. 'Sex_Female_vs_Male' or 'Exposure_Control_minus_VPA'
    'ContrastType', ...       % 'Sex' or 'Exposure'
    'ConditioningLevel', ...  % Exposure (for Sex contrasts) or Sex (for Exposure contrasts)
    'Level1', 'Level2', ...   % Level1 − Level2
    'Estimate', 'SE', 'zValue', 'pValue', ...
    'CI_Lower', 'CI_Upper'});

disp('=== GLMM-based contrasts (iPE) ===');
disp(ContrastTable);


%% Step 14: Visualization of GLMM input estimates using violinplotwhole
% Here we visualize the subject-level median nSFR (T_agg),
% grouped by Sex × Exposure (FaC, MaC, FaV, MaV, FaG, MaG),
% separately for CNIC and RCIC, using your violinplotwhole function.

% 1) Build a plotting table in the format expected by violinplotwhole:
%    variables needed: DIVISION, Level, Group, iDE (here, iDE = nSFR)
T_plot = T_agg(:, {'SubjectID','DIVISION','Sex','Exposure','iDE'});

% Level: all are Adult / High by construction of tbl
T_plot.Level = repmat("High", height(T_plot), 1);   % string array is fine

% Group code: FaC, MaC, FaV, MaV, FaG, MaG
T_plot.Group = strings(height(T_plot),1);
for i = 1:height(T_plot)
    sex = char(T_plot.Sex(i));
    exp = char(T_plot.Exposure(i));

    if strcmp(sex,'Female') && strcmp(exp,'Control')
        T_plot.Group(i) = 'FaC';
    elseif strcmp(sex,'Male') && strcmp(exp,'Control')
        T_plot.Group(i) = 'MaC';
    elseif strcmp(sex,'Female') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'FaV';
    elseif strcmp(sex,'Male') && strcmp(exp,'VPA')
        T_plot.Group(i) = 'MaV';
    elseif strcmp(sex,'Female') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'FaG';
    elseif strcmp(sex,'Male') && strcmp(exp,'GRIN2B')
        T_plot.Group(i) = 'MaG';
    else
        T_plot.Group(i) = "NA";
    end
end

T_plot.Group = categorical(T_plot.Group);

% The function expects the y-variable to be named iDE → alias nSFR as iDE
T_plot.iDE = T_plot.iDE;

% Make sure DIVISION is categorical (should already be, but just in case)
T_plot.DIVISION = categorical(T_plot.DIVISION);

% 2) Call your violinplotwhole for each division at High level
levels = "High";   % we only have High in this analysis

% Lemniscal division (CNIC)
if any(T_plot.DIVISION == "CNIC")
    violinplotwhole(T_plot, "CNIC", levels, glme);
end

% Non-lemniscal division (RCIC)
if any(T_plot.DIVISION == "RCIC")
    violinplotwhole(T_plot, "RCIC", levels, glme);
end

%% --- Shor-term (first 10 STD trials) computations of adaptation magnitudes (use of b coefficients) --- %%

%% Merged table generation

% List of all group-specific TargetTables and their metadata
files = {
    'TargetTable_Control_Female.mat',   'Control',  'Female';
    'TargetTable_Control_Male.mat',     'Control',  'Male';
    'TargetTable_GRIN2B_Female.mat',    'GRIN2B',   'Female';
    'TargetTable_GRIN2B_Male.mat',      'GRIN2B',   'Male';
    'TargetTable_VPA_Female.mat',       'VPA',      'Female';
    'TargetTable_VPA_Male.mat',         'VPA',      'Male'
};

% Initialize combined .Sheet cell array
MergedSheet = {};

for i = 1:size(files,1)
    filename  = files{i,1};
    condition = files{i,2};
    sex       = files{i,3};
    age       = 'Adult';  % All are Adult by design

    % Load each file
    S = load(filename);
    varname = fieldnames(S);
    TargetTable = S.(varname{1});

    % Check .Sheet
    if ~isprop(TargetTable, 'Sheet') || ~iscell(TargetTable.Sheet)
        error('Invalid or missing .Sheet in %s', filename);
    end

    % Identify the header row
    sheet = TargetTable.Sheet;
    header_idx = find(cellfun(@(x) ischar(x) && strcmpi(x, 'ANIMAL'), sheet(:,1)), 1);
    if isempty(header_idx)
        error('No header found in %s', filename);
    end
    header = sheet(header_idx, :);
    data   = sheet(header_idx+1:end, :);

    % Create new columns
    nRows = size(data, 1);
    ConditionCol = repmat({condition}, nRows, 1);
    SexCol       = repmat({sex}, nRows, 1);
    AgeCol       = repmat({age}, nRows, 1);

    % Append new columns to header and data
    newHeader = [header, {'Condition', 'Sex', 'Age'}];
    newData   = [data, ConditionCol, SexCol, AgeCol];

    % Recombine and append
    sheet_with_meta = [newHeader; newData];
    MergedSheet = [MergedSheet; sheet_with_meta];
end

% Final merging steps
header_rows = cellfun(@(x) ischar(x) && strcmpi(x, 'ANIMAL'), MergedSheet(:,1));
header = MergedSheet(find(header_rows, 1), :);
MergedSheetClean = MergedSheet(~header_rows, :);
cleanVarNames = matlab.lang.makeValidName(header);
MergedSheetTable = cell2table(MergedSheetClean, 'VariableNames', cleanVarNames);
% Create combined group label (e.g., "GRIN2B_Adult_Female")
MergedSheetTable.ExpGroup = strcat( ...
    string(MergedSheetTable.Condition), '_', ...
    string(MergedSheetTable.Age), '_', ...
    string(MergedSheetTable.Sex));

%% Short-term (10 first STD trials) adaptation indices

load('MergedSheetTable') % from each experimentla group out of the six animal groups
STA_table = compute_STA_metrics(MergedSheetTable);

% Filter out invalid PercentDrop (NaN or Inf)
%valid_idx   = isfinite(STA_table.PercentDrop);
%clean_table = STA_table(valid_idx, :);


%% Plot the b coefficient for the first 10 STD trials

% Unique divisions and groups
divisions = unique(STA_table.Division);
groups = unique(STA_table.ExpGroup, 'stable');  % preserve order

% Generate distinct colors per group
groupColors = lines(numel(groups));  % Or define manually if you prefer

figure;
for d = 1:numel(divisions)
    subplot(1, numel(divisions), d);
    hold on;

    divData = STA_table(strcmp(STA_table.Division, divisions{d}), :);

    medians = zeros(numel(groups),1);
    sems = zeros(numel(groups),1);

    for g = 1:numel(groups)
        grpData = divData(strcmp(divData.ExpGroup, groups{g}), :);
        bvals = grpData.b;

        medians(g) = median(bvals, 'omitnan');

         % Print median to command window
        fprintf('Division: %s | Group: %s | Median b: %.4f\n', ...
            divisions{d}, groups{g}, medians(g));

        % Bootstrap SE of median
        nBoot = 1000;
        bootMeds = bootstrp(nBoot, @median, bvals);
        sems(g) = std(bootMeds);

        % Scatter individual points (fully opaque)
        x = repmat(g, length(bvals), 1) + 0.1*randn(length(bvals), 1);
        scatter(x, bvals, 20, groupColors(g,:), 'filled', 'MarkerFaceAlpha', 1);
    end

    % Plot median points with error bars, median marker filled with group color
    for g = 1:numel(groups)
        errorbar(g, medians(g), sems(g), 'o', ...
            'MarkerSize', 6, ...            % smaller markers
            'MarkerFaceColor', groupColors(g,:), ...
            'MarkerEdgeColor', 'k', ...
            'LineWidth', 0.8, ...           % thinner error bars
            'CapSize', 8, ...
            'Color', 'k');                  % error bars black
    end

    xlim([0.5, numel(groups)+0.5]);
    xticks(1:numel(groups));
    xticklabels(groups);
    xtickangle(45);
    ylabel('b coefficient (median ± SE)');
    title(divisions{d});
    grid on;
end
sgtitle('Median b coefficient ± SE of median per group by IC division');

%% Plot the b coefficient for the first 10 STD trials

% Unique divisions and groups
divisions = unique(STA_table.Division);
groups = unique(STA_table.ExpGroup, 'stable');  % preserve order

% Generate distinct colors per group
groupColors = lines(numel(groups));  % Or define manually if you prefer

figure;
for d = 1:numel(divisions)
    subplot(1, numel(divisions), d);
    hold on;

    divData = STA_table(strcmp(STA_table.Division, divisions{d}), :);

    medians = zeros(numel(groups),1);
    sems = zeros(numel(groups),1);

    for g = 1:numel(groups)
        grpData = divData(strcmp(divData.ExpGroup, groups{g}), :);
        bvals = grpData.b;

        medians(g) = median(bvals, 'omitnan');

        % Print median to command window
        fprintf('Division: %s | Group: %s | Median b: %.4f\n', ...
            divisions{d}, groups{g}, medians(g));

        % Bootstrap SE of median
        nBoot = 1000;
        bootMeds = bootstrp(nBoot, @median, bvals);
        sems(g) = std(bootMeds);

        % Scatter individual points (fully opaque)
        x = repmat(g, length(bvals), 1) + 0.1*randn(length(bvals), 1);
        scatter(x, bvals, 20, groupColors(g,:), 'filled', 'MarkerFaceAlpha', 1);
    end

    % Plot median points with error bars, median marker filled with group color
    for g = 1:numel(groups)
        errorbar(g, medians(g), sems(g), 'o', ...
            'MarkerSize', 6, ...
            'MarkerFaceColor', groupColors(g,:), ...
            'MarkerEdgeColor', 'k', ...
            'LineWidth', 0.8, ...
            'CapSize', 8, ...
            'Color', 'k');
    end

    xlim([0.5, numel(groups)+0.5]);
    xticks(1:numel(groups));
    xticklabels(groups);
    xtickangle(45);
    ylabel('b coefficient (median ± SE)');
    title(divisions{d});
    grid on;
end
sgtitle('Median b coefficient ± SE of median per group by IC division');

%% Match column names and merge STA_table with MergedSheetTable

% Rename uppercase headers in MergedSheetTable to match STA_table
MergedSheetTable.Properties.VariableNames{strcmp(MergedSheetTable.Properties.VariableNames, 'ANIMAL')}   = 'Animal';
MergedSheetTable.Properties.VariableNames{strcmp(MergedSheetTable.Properties.VariableNames, 'UNIT')}     = 'Unit';
MergedSheetTable.Properties.VariableNames{strcmp(MergedSheetTable.Properties.VariableNames, 'DIVISION')} = 'Division';

% Merge tables
combinedTable = innerjoin(STA_table, MergedSheetTable, ...
    'Keys', {'Animal', 'Unit', 'Division'});

%% Spearman correlation between b and iRA by division and group

divisions = unique(combinedTable.Division, 'stable');
groups = unique(combinedTable.ExpGroup_MergedSheetTable, 'stable');

results = table('Size',[0 5], ...
    'VariableTypes', {'string','string','double','double','double'}, ...
    'VariableNames', {'Division','ExpGroup','SpearmanRho','pValue','R_squared'});

% Matrix to store rho and n for later plotting
rhoMat = NaN(numel(groups), numel(divisions));
nMat   = NaN(numel(groups), numel(divisions));

for d = 1:numel(divisions)
    for g = 1:numel(groups)
        idx = strcmp(combinedTable.Division, divisions{d}) & ...
              strcmp(combinedTable.ExpGroup_MergedSheetTable, groups{g});

        b_sub = combinedTable.b(idx);
        iRA_sub = combinedTable.iRA(idx);

        valid = ~isnan(b_sub) & ~isnan(iRA_sub);
        n = sum(valid);
        if n < 5
            continue;
        end

        [rho, pval] = corr(b_sub(valid), iRA_sub(valid), 'Type', 'Spearman');

        % Store in results table
        newRow = {divisions{d}, groups{g}, rho, pval, rho^2};
        results = [results; newRow];

        % Store in matrix
        rhoMat(g,d) = rho;
        nMat(g,d)   = n;
    end
end

disp(results);
% writetable(results, 'SpearmanCorrelationResults.csv');
% FDR correction
[~, ~, ~, adj_p] = fdr_bh(results.pValue);

% Add to table
results.FDR_pValue = adj_p;

% Optional: flag significance
results.FDR_significant = adj_p < 0.05;

% Display updated results
disp(results);

%% Plot grouped bar chart with Fisher-z confidence intervals

figure;
b = bar(rhoMat, 'grouped');
set(gca, 'XTickLabel', groups, 'XTickLabelRotation', 45);
ylabel('Spearman \rho');
ylim([-1 1]);
legend(divisions, 'Location','best');
title('Spearman \rho between b and iRA by Group & Division');
grid on;
hold on;

for d = 1:numel(divisions)
    xCenters = b(d).XEndPoints;
    for g = 1:numel(groups)
        r = rhoMat(g,d);
        n = nMat(g,d);
        if isnan(r) || n < 4
            continue;
        end

        % Fisher-z transform
        z = atanh(r);
        SEz = 1 / sqrt(n - 3);
        CIz = 1.96 * SEz;
        lo = tanh(z - CIz);
        hi = tanh(z + CIz);
        errLow = r - lo;
        errHigh = hi - r;

        % Plot error bar
        errorbar(xCenters(g), r, errLow, errHigh, ...
            'LineWidth', 0.5, ...
            'CapSize', 6, ...
            'Color', [0.4 0.4 0.4]);
    end
end

hold off;


%% ------- Merged heatmaps for pairwise comparisons ------ %%

%% Merged heatmaps for predictive components (iSSA, iRA, iDE) across CNIC and RCIC

% Load and rename variables from .mat files
load('resultsTable_iSSA.mat'); resultsTableiSSA = resultsTable_iSSA;
load('resultsTable_iRA.mat');  resultsTableiRA  = resultsTable_iRA;
load('resultsTable_iDE.mat');  resultsTableiDE  = resultsTable_iDE;

% Combine all group names from the three result tables
allGroups = unique([resultsTableiSSA.Group1; resultsTableiSSA.Group2; ...
                    resultsTableiRA.Group1;  resultsTableiRA.Group2; ...
                    resultsTableiDE.Group1;  resultsTableiDE.Group2]);

% Set response labels and loaded tables
responses = {'iSSA', 'iRA', 'iDE'};
resultsTables = {resultsTableiSSA, resultsTableiRA, resultsTableiDE};
divisions = {'CNIC', 'RCIC'};

% Determine global min/max p-values across all divisions for color scaling
minP = inf; maxP = -inf;
for r = 1:numel(responses)
    tbl = resultsTables{r};
    for d = 1:numel(divisions)
        rows = strcmp(tbl.Division, divisions{d});
        pVals = tbl.Corrected_pValue(rows);
        minP = min([minP; pVals]);
        maxP = max([maxP; pVals]);
    end
end

% Create heatmaps for each division
for d = 1:numel(divisions)
    division = divisions{d};
    
    % Initialize 3D matrix for the heatmap
    heatMatrix = NaN(numel(allGroups), numel(allGroups), numel(responses));
    
    for r = 1:numel(responses)
        tbl = resultsTables{r};
        subTbl = tbl(strcmp(tbl.Division, division), :);
        
        for k = 1:height(subTbl)
            g1 = subTbl.Group1{k};
            g2 = subTbl.Group2{k};
            i = find(strcmp(allGroups, g1));
            j = find(strcmp(allGroups, g2));
            if ~isempty(i) && ~isempty(j)
                heatMatrix(i, j, r) = subTbl.Corrected_pValue(k);
            end
        end
    end
    
    % Plot merged heatmap
    figure('Name', [division, ' - Predictive Components'], 'NumberTitle', 'off');
    imagesc(mean(heatMatrix, 3, 'omitnan'));
    colormap(parula); colorbar; caxis([minP, maxP]);
    title([division, ' - Predictive Components']);
    xticks(1:numel(allGroups)); yticks(1:numel(allGroups));
    xticklabels(allGroups); yticklabels(allGroups);
    xtickangle(45); axis square;
    
    % Annotate with individual corrected p-values, color-coded per component
    for i = 1:numel(allGroups)
        for j = 1:numel(allGroups)
            for r = 1:numel(responses)
                val = heatMatrix(i,j,r);
                if ~isnan(val)
                    yOffset = 0.25 * (r - 2);  % -0.25 for iSSA, 0 for iRA, +0.25 for iDE
                    switch responses{r}
                        case 'iDE', color = 'r';
                        case 'iSSA', color = 'b';
                        case 'iRA', color = 'g';
                    end
                    text(j, i + yOffset, sprintf('%.3f', val), ...
                         'HorizontalAlignment', 'center', 'FontSize', 7, 'Color', color);
                end
            end
        end
    end
end

%% Merged heatmaps for evoked responses (DEV, CTR, STD) across CNIC and RCIC

% Load and rename variables from .mat files
load('resultsTable_DEV.mat'); resultsTableDEV = resultsTable_DEV;
load('resultsTable_CTR.mat'); resultsTableCTR = resultsTable_CTR;
load('resultsTable_STD.mat'); resultsTableSTD = resultsTable_STD;

% Combine all group names from the three result tables
allGroups = unique([resultsTableDEV.Group1; resultsTableDEV.Group2; ...
                    resultsTableCTR.Group1; resultsTableCTR.Group2; ...
                    resultsTableSTD.Group1; resultsTableSTD.Group2]);

% Set response labels and loaded tables
responses = {'DEV', 'CTR', 'STD'};
resultsTables = {resultsTableDEV, resultsTableCTR, resultsTableSTD};
divisions = {'CNIC', 'RCIC'};

% Determine global min/max p-values across all divisions for color scaling
minP = inf; maxP = -inf;
for r = 1:numel(responses)
    tbl = resultsTables{r};
    for d = 1:numel(divisions)
        rows = strcmp(tbl.Division, divisions{d});
        pVals = tbl.Corrected_pValue(rows);
        minP = min([minP; pVals]);
        maxP = max([maxP; pVals]);
    end
end

% Create heatmaps for each division
for d = 1:numel(divisions)
    division = divisions{d};
    
    % Initialize 3D matrix for the heatmap
    heatMatrix = NaN(numel(allGroups), numel(allGroups), numel(responses));
    
    for r = 1:numel(responses)
        tbl = resultsTables{r};
        subTbl = tbl(strcmp(tbl.Division, division), :);
        
        for k = 1:height(subTbl)
            g1 = subTbl.Group1{k};
            g2 = subTbl.Group2{k};
            i = find(strcmp(allGroups, g1));
            j = find(strcmp(allGroups, g2));
            if ~isempty(i) && ~isempty(j)
                heatMatrix(i, j, r) = subTbl.Corrected_pValue(k);
            end
        end
    end
    
    % Plot merged heatmap
    figure('Name', [division, ' - Evoked Responses'], 'NumberTitle', 'off');
    imagesc(mean(heatMatrix, 3, 'omitnan'));
    colormap(parula); colorbar; caxis([minP, maxP]);
    title([division, ' - Evoked Responses']);
    xticks(1:numel(allGroups)); yticks(1:numel(allGroups));
    xticklabels(allGroups); yticklabels(allGroups);
    xtickangle(45); axis square;
    
    % Annotate with individual corrected p-values, color-coded per response
    for i = 1:numel(allGroups)
        for j = 1:numel(allGroups)
            for r = 1:numel(responses)
                val = heatMatrix(i,j,r);
                if ~isnan(val)
                    yOffset = 0.25 * (r - 2);  % -0.25 for STD, 0 for CTR, +0.25 for DEV
                    switch responses{r}
                        case 'DEV', color = 'r';
                        case 'CTR', color = 'g';
                        case 'STD', color = 'b';
                    end
                    text(j, i + yOffset, sprintf('%.3f', val), ...
                         'HorizontalAlignment', 'center', 'FontSize', 7, 'Color', color);
                end
            end
        end
    end
end



%% ---------------- First revision for comms biol ------------------- %%

%% ===== STEP 1 (REVISED): Subject-level aggregation + ID hygiene + diagnostics =====
% Load and filter data
dataStruct = load('UpdatedFinalTableRevision.mat');
tbl = dataStruct.UpdatedFinalTable;

assert(all(ismember({'DIVISION','Sex','Age','Exposure','Level','nSFR','SubjectID'}, tbl.Properties.VariableNames)), ...
    'Missing required variables.');

% Adult + High filter (as before)
tbl = tbl(strcmp(tbl.Age, 'Adult') & strcmp(tbl.Level, 'High'), :);

% Harmonize types
tbl.Sex      = categorical(tbl.Sex);
tbl.Age      = categorical(tbl.Age);
tbl.Exposure = categorical(tbl.Exposure);
tbl.DIVISION = categorical(tbl.DIVISION);
tbl.Level    = categorical(tbl.Level);

% ---------- NEW: SubjectID hygiene (trim whitespace, unify case) ----------
tbl.SubjectID = categorical(strtrim(string(tbl.SubjectID)));

% Safety check
reqVars = {'SubjectID','DIVISION','Sex','Exposure','nSFR'};
assert(all(ismember(reqVars, tbl.Properties.VariableNames)), 'Missing required variables for STEP 1.');

% ---------- Build subject-level table: one row per Subject × Division × Sex × Exposure ----------
T_subject = groupsummary(tbl, {'SubjectID','DIVISION','Sex','Exposure'}, {'median','numel'}, 'nSFR');
T_subject.Properties.VariableNames(strcmp(T_subject.Properties.VariableNames,'median_nSFR')) = {'SubjectMedian_nSFR'};
T_subject.Properties.VariableNames(strcmp(T_subject.Properties.VariableNames,'GroupCount'))  = {'NNeurons_inStratum'};
T_subject = T_subject(T_subject.NNeurons_inStratum > 0, :);

% ---------- Per-division subject counts (animals per Sex × Exposure within each Division) ----------
fprintf('\n[STEP 1] Subjects per group by Division (Adult, High):\n');
for dv = categories(T_subject.DIVISION)'
    Tdv = T_subject(T_subject.DIVISION == dv{1}, :);
    ST  = groupcounts(Tdv, {'Sex','Exposure'}); % returns GroupCount
    ST.Properties.VariableNames(strcmp(ST.Properties.VariableNames,'GroupCount')) = {'Nsubjects'};
    fprintf('  Division: %s\n', dv{1});
    disp(ST(:, {'Sex','Exposure','Nsubjects'}));
end

% ---------- One-to-one mapping check: each SubjectID must map to exactly one Sex and one Exposure ----------
[G_subj, subj_ids] = findgroups(tbl.SubjectID);

% Count unique Sex and unique Exposure per SubjectID
nUniqueSex      = splitapply(@(x) numel(unique(x)),      tbl.Sex,      G_subj);
nUniqueExposure = splitapply(@(x) numel(unique(x)), tbl.Exposure, G_subj);

U = table(subj_ids, nUniqueSex, nUniqueExposure, ...
          'VariableNames', {'SubjectID','nUniqueSex','nUniqueExposure'});

flag = U(U.nUniqueSex > 1 | U.nUniqueExposure > 1, :);
if ~isempty(flag)
    fprintf('\n[WARNING] SubjectID(s) map to >1 Sex and/or >1 Exposure (check labels):\n');
    disp(flag);
end

% ---------- Overall unique-subject totals by Sex × Exposure (across both divisions) ----------
G_overall = unique(tbl(:,{'SubjectID','Sex','Exposure'}),'rows');
T_overall = groupcounts(G_overall, {'Sex','Exposure'});
T_overall.Properties.VariableNames(strcmp(T_overall.Properties.VariableNames,'GroupCount')) = {'Nsubjects_overall'};
fprintf('\n[STEP 1] Overall unique-subject totals by Sex × Exposure (Adult, High; across CNIC+RCIC):\n');
disp(T_overall);

% Store for downstream steps
T_subject_all = T_subject;   % subject-level medians for permutation tests


%% ===== STEP 2: Cluster-based permutation tests (subject-level shuffling, neuron-level data) =====
% Purpose: keeps neurons as units but shuffles subjects as clusters to test robustness to nesting.
% Output: same contrasts as before (Sex matched by Exposure, Exposure matched by Sex).
% Each p_cluster is based on subject-level permutation (10,000 by default).

fprintf('\n====================== STEP 2: CLUSTER-BASED PERMUTATION TESTS ======================\n');

% Parameters
nPerm = 10000;   % number of subject-level permutations
divisions = categories(tbl.DIVISION);
exposures = {'Control','GRIN2B','VPA'};
sexes     = {'Female','Male'};

Results_clusterPerm = struct();

for d = 1:numel(divisions)
    dv = divisions{d};
    tbl_div = tbl(tbl.DIVISION==dv, :);

    fprintf('\n[STEP 2] Division: %s | Performing subject-level cluster permutations (%d iters)...\n', dv, nPerm);

    %% --------- Family A: Sex effects matched by Exposure (Female_x vs Male_x) ---------
    famA_rows = {};
    for e = 1:numel(exposures)
        ex = exposures{e};
        g1mask = tbl_div.Sex=="Female" & tbl_div.Exposure==ex;
        g2mask = tbl_div.Sex=="Male"   & tbl_div.Exposure==ex;

        if sum(g1mask)>0 && sum(g2mask)>0
            res = cluster_perm_neuronwise(tbl_div, 'Sex', 'nSFR', 'SubjectID', "Female", "Male", nPerm, ex);
            famA_rows(end+1,:) = {dv, ['Female_' ex], ['Male_' ex], res.MedianDiff, res.p_cluster, res.S_bits}; %#ok<AGROW>
        end
    end
    TA = cell2table(famA_rows, 'VariableNames', ...
        {'Division','Group1','Group2','MedianDiff','p_cluster','S_bits'});
    TA.p_holm = holm_adjust(TA.p_cluster);

    %% --------- Family B: Exposure effects matched by Sex (Sex_Control vs Sex_Model) ---------
    famB_rows = {};
    for s = 1:numel(sexes)
        sx = sexes{s};
        for e = 2:numel(exposures)
            ex = exposures{e}; % GRIN2B, VPA
            res = cluster_perm_neuronwise(tbl_div, 'Exposure', 'nSFR', 'SubjectID', "Control", ex, nPerm, sx);
            famB_rows(end+1,:) = {dv, [sx '_Control'], [sx '_' ex], res.MedianDiff, res.p_cluster, res.S_bits}; %#ok<AGROW>
        end
    end
    TB = cell2table(famB_rows, 'VariableNames', ...
        {'Division','Group1','Group2','MedianDiff','p_cluster','S_bits'});
    TB.p_holm = holm_adjust(TB.p_cluster);

    % Store & print
    Results_clusterPerm.(dv).Sex_byExposure = TA;
    Results_clusterPerm.(dv).Exposure_bySex = TB;

    fprintf('\n[STEP 2] %s — Sex (matched by Exposure):\n', dv);
    disp(TA);

    fprintf('\n[STEP 2] %s — Exposure (matched by Sex):\n', dv);
    disp(TB);
end

fprintf('\n====================== END OF STEP 2 ======================\n');


%% The permutation analyses went wrong because of non-significant results due to underpower (all neurons behave similarly
% within subjects, so the clustering goes under-powered

%% ----------- Checkings to reinforce our neural-level analysis -------------------- %%

%% STEP 3A — LME-based ICC & Design Effect (Adult × High)
% Uses REML random-intercept LME, consistent with Step 3B.
% Computes ICC = Var(SubjectID) / (Var(SubjectID) + Var(Residual))
% and DE ≈ 1 + (m̄ - 1)*ICC with m̄ = mean neurons per animal.

clear; clc;
fprintf('=== STEP 3A — LME-based ICC & Design Effect (Adult×High) [NOR_STD] ===\n');

%% ===== Load & filter =====
S = load('UpdatedFinalTable.mat');
T = S.UpdatedFinalTable;

% Adult + High only
T = T(strcmp(T.Age,'Adult') & strcmp(T.Level,'High'), :);

% Types + hygiene
T.DIVISION  = categorical(T.DIVISION);
T.Sex       = categorical(T.Sex);
T.Exposure  = categorical(T.Exposure);
T.SubjectID = categorical(strtrim(string(T.SubjectID)));

% Keep finite response
T = T(isfinite(T.NOR_STD), :);

fprintf('Data after Adult×High filter: N rows = %d, N animals = %d\n', ...
    height(T), numel(categories(T.SubjectID)));

%% ===== Define group variable =====
T.Group = categorical(string(T.Sex) + "_" + string(T.Exposure));

%% ===== Loop over DIVISION × group pairs =====
divisions = categories(removecats(T.DIVISION));
minSubjectsPerGroup = 1;
minNeuronsPerGroup  = 5;

for d = 1:numel(divisions)
    div = divisions{d};
    Td  = T(T.DIVISION == div, :);

    ug = categories(removecats(Td.Group));
    if numel(ug) < 2
        warning('DIVISION %s has <2 groups after filtering; skipping.', div);
        continue;
    end

    grpPairs = nchoosek(1:numel(ug),2);

    for pp = 1:size(grpPairs,1)
        g1 = ug{grpPairs(pp,1)};
        g2 = ug{grpPairs(pp,2)};

        % Subset to exactly two groups within this DIVISION
        S2 = Td(Td.Group==g1 | Td.Group==g2, :);
        S2.Group = removecats(S2.Group);
        S2.SubjectID = removecats(S2.SubjectID);

        cats = categories(S2.Group);
        if numel(cats) ~= 2
            warning('Pairing yielded %d groups (expected 2). Skipping %s vs %s in %s.', ...
                numel(cats), g1, g2, div);
            continue;
        end

        % Minimal data checks
        nNeurons = zeros(1,2);
        nSubjects = zeros(1,2);
        for ii = 1:2
            idx = S2.Group == cats{ii};
            nNeurons(ii)  = sum(idx);
            nSubjects(ii) = numel(categories(S2.SubjectID(idx)));
        end
        if any(nSubjects < minSubjectsPerGroup) || any(nNeurons < minNeuronsPerGroup)
            warning('Low counts in %s: %s vs %s (subjects/group: %s; neurons/group: %s). Skipping.', ...
                div, cats{1}, cats{2}, mat2str(nSubjects), mat2str(nNeurons));
            continue;
        end

        fprintf('\nDIV=%s | %s vs %s | Nneur=%d | Nsub=%d\n', ...
            div, g1, g2, height(S2), numel(categories(S2.SubjectID)));

        % ---- Fit REML model with random intercept ----
        try
            m = fitlme(S2, 'NOR_STD ~ 1 + Group + (1|SubjectID)', ...
                'FitMethod','REML', 'DummyVarCoding','effects');

            % ----- Extract variance components robustly -----
            VarID  = NaN;
            VarRes = NaN;

            % Preferred: 2-output form (SubjectID variance in vc{1}, residual in details)
            try
                [vc, details] = covarianceParameters(m);
                if iscell(vc) && ~isempty(vc) && isnumeric(vc{1}) && isscalar(vc{1})
                    VarID = vc{1};
                end
                if isnumeric(details) && isscalar(details) && details > 0
                    VarRes = details;
                end
            catch
                % Fallback: 1-output form + MSE for residual variance
                vc = covarianceParameters(m);
                if iscell(vc) && ~isempty(vc) && isnumeric(vc{1}) && isscalar(vc{1})
                    VarID = vc{1};
                end
            end

            % Hard fallback for residual variance (R2022b has MSE)
            if ~(isfinite(VarRes) && VarRes > 0)
                try
                    VarRes = m.MSE;
                catch
                    % leave as NaN
                end
            end

            % Clamp tiny negatives due to numerics
            if isfinite(VarID)  && VarID  < 0, VarID  = 0; end
            if isfinite(VarRes) && VarRes < 0, VarRes = 0; end

            % ----- Compute ICC and DE -----
            % Mean neurons per animal (cluster size), no groupsummary (R2022b-safe)
            [G,~]  = findgroups(S2.SubjectID);
            counts = splitapply(@numel, S2.NOR_STD, G);
            mbar   = mean(counts);

            ICC = NaN; DE = NaN;
            denom = VarID + VarRes;
            if isfinite(VarID) && isfinite(denom) && denom > 0
                ICC = VarID / denom;
            end
            if isfinite(ICC) && isfinite(mbar)
                DE  = 1 + (mbar - 1) * ICC;
            end

            fprintf('  Var(SubjectID)=%.6f | Var(Residual)=%.6f | ICC_LME=%.6f\n', ...
                VarID, VarRes, ICC);
            fprintf('  Mean neurons/animal = %.2f | Design Effect (DE)=%.3f\n', ...
                mbar, DE);

        catch ME
            fprintf('  Model failed (%s)\n', ME.message);
            continue;
        end
    end
end

fprintf('\nDone (Step 3A, fixed LME version — NOR_STD).\n');


%% STEP 3B — Random-intercept model check ONLY (Adult × High)
% This script assesses whether including SubjectID as a random intercept
% improves model fit for NOR_STD across DIVISION × group pairs.

clear; clc;
fprintf('=== Random-intercept model check (Adult × High) ===\n');

%% ===== Load + Adult/High filter + ID hygiene =====
dataStruct = load('UpdatedFinalTable.mat');
tbl = dataStruct.UpdatedFinalTable;

% --- Required variables ---
reqVars = {'DIVISION','Sex','Age','Exposure','Level','NOR_STD','SubjectID'};
assert(all(ismember(reqVars, tbl.Properties.VariableNames)), ...
    'Missing required variables. Required: %s', strjoin(reqVars, ', '));

% --- Adult + High only ---
tbl = tbl(strcmp(tbl.Age,'Adult') & strcmp(tbl.Level,'High'), :);

% --- Harmonize types ---
tbl.Sex       = categorical(tbl.Sex);
tbl.Age       = categorical(tbl.Age);
tbl.Exposure  = categorical(tbl.Exposure);
tbl.DIVISION  = categorical(tbl.DIVISION);
tbl.Level     = categorical(tbl.Level);
tbl.SubjectID = categorical(strtrim(string(tbl.SubjectID)));
tbl.AnimalID  = categorical(tbl.SubjectID); % alias

% --- Clean NOR_STD values ---
tbl = tbl(isfinite(tbl.NOR_STD), :);

fprintf('Data after Adult×High filter: N rows = %d, N animals = %d\n', ...
    height(tbl), numel(categories(tbl.AnimalID)));

%% ===== Define grouping for pairwise contrasts =====
groupingMode = 'SexExposure'; % or 'Exposure', 'Sex'

switch groupingMode
    case 'SexExposure'
        tbl.Group = categorical(string(tbl.Sex) + "_" + string(tbl.Exposure));
    case 'Exposure'
        tbl.Group = categorical(tbl.Exposure);
    case 'Sex'
        tbl.Group = categorical(tbl.Sex);
    otherwise
        error('Unknown groupingMode: %s', groupingMode);
end

tbl = tbl(~isundefined(tbl.DIVISION) & ~isundefined(tbl.Group) & ~isundefined(tbl.AnimalID), :);

%% ===== Run Random-Intercept Model Checks =====
divisions = categories(removecats(tbl.DIVISION));
minSubjectsPerGroup = 1; 
minNeuronsPerGroup  = 5;

for d = 1:numel(divisions)
    div = divisions{d};
    Tdiv = tbl(tbl.DIVISION == div, :);

    ug = categories(removecats(Tdiv.Group));
    if numel(ug) < 2
        warning('DIVISION %s has <2 groups after filtering; skipping.', div);
        continue;
    end

    grpPairs = nchoosek(1:numel(ug), 2);

    for pp = 1:size(grpPairs,1)
        g1 = ug{grpPairs(pp,1)};
        g2 = ug{grpPairs(pp,2)};

        % Subset and re-level
        S = Tdiv(Tdiv.Group==g1 | Tdiv.Group==g2, :);
        S.Group = removecats(categorical(string(S.Group)));
        cats = categories(S.Group);
        if numel(cats) ~= 2
            warning('Pairing yielded %d groups (expected 2). Skipping %s vs %s in %s.', ...
                numel(cats), g1, g2, div);
            continue;
        end

        % --- Ensure sufficient data per group ---
        nNeurons  = zeros(1,2);
        nSubjects = zeros(1,2);
        for ii = 1:2
            idx = S.Group == cats{ii};
            nNeurons(ii)  = sum(idx);
            nSubjects(ii) = numel(unique(S.SubjectID(idx)));
        end
        if any(nSubjects < minSubjectsPerGroup) || any(nNeurons < minNeuronsPerGroup)
            warning('Low counts in %s: %s vs %s (subjects/group: %s; neurons/group: %s). Skipping.', ...
                div, cats{1}, cats{2}, mat2str(nSubjects), mat2str(nNeurons));
            continue;
        end

        fprintf('\nDIV=%s | %s vs %s | Nneur=%s | Nsub=%s\n', ...
            div, cats{1}, cats{2}, mat2str(nNeurons), mat2str(nSubjects));

        % --- Run the random-intercept model check ---
        mcOut = random_intercept_model_check_NOR_STD(S, 'IDVar', 'SubjectID');

        % --- Compute m̄ = mean neurons/animal (R2022b-safe; no groupsummary) ---
        [~,~,g] = unique(S.SubjectID);      % integer group index per row
        counts  = accumarray(g,1);          % neurons per SubjectID
        mbar    = mean(counts);

        % --- Print results (with LME-based ICC and DE) ---
        if isnan(mcOut.VarID), varTxt = 'NaN'; else, varTxt = sprintf('%.6f', mcOut.VarID); end
        if isnan(mcOut.VarResidual), resTxt = 'NaN'; else, resTxt = sprintf('%.6f', mcOut.VarResidual); end
        if isnan(mcOut.ICC_LME), iccTxt = 'NaN'; else, iccTxt = sprintf('%.6f', mcOut.ICC_LME); end

        fprintf('  Var(Intercept|SubjectID) = %s | Var(Residual) = %s | ICC_LME = %s\n', varTxt, resTxt, iccTxt);
        fprintf('  AIC(no RE) = %.2f | AIC(RE) = %.2f | ΔAIC = %.2f\n', ...
            mcOut.AIC_noRE, mcOut.AIC_RE, mcOut.AIC_RE - mcOut.AIC_noRE);
        fprintf('  BIC(no RE) = %.2f | BIC(RE) = %.2f | ΔBIC = %.2f\n', ...
            mcOut.BIC_noRE, mcOut.BIC_RE, mcOut.BIC_RE - mcOut.BIC_noRE);
        fprintf('  LRT p-value (ML) = %.4f\n', mcOut.LRT_p);

        if ~isnan(mcOut.ICC_LME)
            DE_est = 1 + (mbar - 1) * mcOut.ICC_LME;
            fprintf('  Estimated DE from LME ICC: DE ≈ %.3f (m̄=%.2f)\n', DE_est, mbar);
        end

        if (mcOut.LRT_p < 0.05) && ~isnan(mcOut.VarID) && (mcOut.VarID > 0)
            fprintf('  → Evidence for random intercept (nonzero subject-level variance).\n');
        elseif (mcOut.LRT_p < 0.05) && ~isnan(mcOut.VarID) && (mcOut.VarID == 0)
            fprintf('  → LRT significant but variance at boundary (σ²≈0): no meaningful clustering.\n');
        elseif (mcOut.LRT_p < 0.05) && isnan(mcOut.VarID)
            fprintf('  → LRT significant but variance unavailable; used ICC_LME to gauge clustering.\n');
        else
            fprintf('  → No significant gain from random intercept (neurons independent across animals).\n');
        end
    end
end

fprintf('\nDone (Step 3B).\n');


%%

% === Minimal one-off diagnostic ===
% Assume S is already the subset table for one pair (two groups in one DIVISION)

ID = 'SubjectID';

m1 = fitlme(S, sprintf('NOR_STD ~ 1 + Group + (1|%s)', ID), ...
    'FitMethod','REML','DummyVarCoding','effects');

disp('>>> disp(m1) (truncated to show covariance part):');
disp(m1);   % MATLAB prints a section with "Estimated Covariance Parameters"

disp('>>> covarianceParameters(m1): 1-output form');
vc1 = covarianceParameters(m1);
whos vc1
if iscell(vc1)
    fprintf('vc1 is a cell with %d elements\n', numel(vc1));
    for i=1:numel(vc1)
        fprintf('  vc1{%d} class = %s\n', i, class(vc1{i}));
        if istable(vc1{i})
            disp(vc1{i});
        else
            disp(vc1{i});
        end
    end
elseif istable(vc1)
    disp(vc1)
else
    disp(vc1)
end

disp('>>> covarianceParameters(m1): 2-output form (if supported)');
try
    [vc2, details] = covarianceParameters(m1);
    whos vc2 details
    disp(vc2)
    disp(details)
catch ME
    fprintf('2-output form threw: %s\n', ME.message);
end
