[NET-193] New Rankings...and other junk?

+nla/+net/+result/+plot/PermutationTestPlotter.m

@@ -53,7 +53,11 @@ function plotProbabilityVsNetworkSize(obj, parameters, axes, plot_title)
             second_title = sprintf('Check if P-values correlate with net-pair size\n(corr: p = %.2f, r = %.2f)', p_values, rho);
             setTitle(axes, second_title, true);
             lims = ylim(axes);
-            ylim(axes, [0 lims(2)]);
+            if lims(2) < 0
+                ylim(axes, [lims(2) 0]);
+            else
+                ylim(axes, [0 lims(2)]);
+            end
         end
 
         function plotProbabilityHistogram(obj, axes, histogram_data, statistic_input, no_permutations_network_result, test_method,...

+nla/+net/+result/NetworkTestResult.m

@@ -29,7 +29,7 @@
         ranking_statistic = "" 
         within_network_pair = false % Results for within-network-pair tests
         full_connectome = false % Results for full connectome tests (formerly 'experiment wide')
-        no_permutations = false % Results for the network tests with no permutations (the 'observed' results)    
+        no_permutations = false % Results for the network tests with no permutations (the 'observed' results)
         permutation_results = struct() % Results for each permutation test used to calculate p-values for the test methods   
     end
 
@@ -203,12 +203,13 @@ function concatenateResult(obj, other_object)
         function createResultsStorage(obj, test_options, number_of_networks, test_specific_statistics)
             %CREATERESULTSSTORAGE Create the substructures for the methods chosen
 
-            % We're just doing them all! The ranking is so short compared to the data collection, just do them all
-            setup_test_methods = ["no_permutations", "full_connectome", "within_network_pair"];
-
             % create the results containers. This replaces the false boolean with a struct of TriMatrices
-            for test_method_index = 1:numel(setup_test_methods)
-                obj.createPValueTriMatrices(number_of_networks, setup_test_methods(test_method_index));
+            for test_method_index = 1:numel(obj.test_methods)
+                if isequal(obj.(obj.test_methods(test_method_index)), false) &&...
+                    isequal(test_options.(obj.test_methods(test_method_index)), true)
+                    obj.(obj.test_methods(test_method_index)) = struct();
+                    obj.createPValueTriMatrices(number_of_networks, obj.test_methods(test_method_index));
+                end
             end
             % This creates all the permutations and test specific stats (chi2, t, w, etc)
             obj.createTestSpecificResultsStorage(number_of_networks, test_specific_statistics);
@@ -220,8 +221,10 @@ function createTestSpecificResultsStorage(obj, number_of_networks, test_specific
             import nla.TriMatrix nla.TriMatrixDiag
 
             obj.permutation_results.p_value_permutations = TriMatrix(number_of_networks, TriMatrixDiag.KEEP_DIAGONAL);
-            obj.permutation_results.single_sample_p_value_permutations = TriMatrix(number_of_networks,...
-                TriMatrixDiag.KEEP_DIAGONAL);
+            if ~any(strcmp(obj.test_name, obj.noncorrelation_input_tests))
+                obj.permutation_results.single_sample_p_value_permutations = TriMatrix(number_of_networks,...
+                    TriMatrixDiag.KEEP_DIAGONAL);
+            end
 
             for statistic_index = 1:numel(test_specific_statistics)
                 test_statistic = test_specific_statistics(statistic_index);
@@ -236,13 +239,11 @@ function createPValueTriMatrices(obj, number_of_networks, test_method)
 
             import nla.TriMatrix nla.TriMatrixDiag
 
-            % I could've looped this, too. Just copy/paste from earlier, so it stays. Plus, this is in every test 
-            % regardless of test or method
-            obj.(test_method) = struct();
-            obj.(test_method).p_value = TriMatrix(number_of_networks, TriMatrixDiag.KEEP_DIAGONAL);
-            obj.(test_method).single_sample_p_value = TriMatrix(number_of_networks, TriMatrixDiag.KEEP_DIAGONAL);
+            obj.(test_method).p_value = TriMatrix(number_of_networks, TriMatrixDiag.KEEP_DIAGONAL); % p-value by statistic rank
             obj.(test_method).statistic_p_value = TriMatrix(number_of_networks, TriMatrixDiag.KEEP_DIAGONAL); % p-value by statistic rank
-            obj.(test_method).statistic_single_sample_p_value = TriMatrix(number_of_networks, TriMatrixDiag.KEEP_DIAGONAL); % p-value by statistic rank
+            if ~isequal(test_method, "full_connectome") && ~any(strcmp(obj.test_name, obj.noncorrelation_input_tests))
+                obj.(test_method).single_sample_p_value = TriMatrix(number_of_networks, TriMatrixDiag.KEEP_DIAGONAL);
+            end
             %Cohen's D results
             obj.(test_method).d = TriMatrix(number_of_networks, TriMatrixDiag.KEEP_DIAGONAL);
         end
@@ -458,7 +459,9 @@ function withinNetworkPairPlotting(obj, network_atlas, edge_test_options, edge_t
             name = sprintf("%s %s P < %.2g (%s)", title_prefix, obj.test_display_name, p_value_max, p_breakdown_labels);
         end
 
-        function [number_of_tests, sig_count_mat, names] = appendSignificanceMatrix(obj, number_of_tests, sig_count_mat, names, sig, name)
+        function [number_of_tests, sig_count_mat, names] = appendSignificanceMatrix(...
+            obj, number_of_tests, sig_count_mat, names, sig, name...
+        )
             number_of_tests = number_of_tests + 1;
             sig_count_mat.v = sig_count_mat.v + sig.v;
             names = [names name];

+nla/+net/CohenDTest.m

@@ -30,8 +30,12 @@
                         (std(edge_test_results.coeff.v).^2)));
 
                     result_object.no_permutations.d.set(row, column, single_sample_d);
-                    result_object.within_network_pair.d.set(row, column, single_sample_d);
-                    result_object.full_connectome.d.set(row, column, d);
+                    if isprop(result_object, "full_connectome") && isequal(result_object.full_connectome, true)
+                        result_object.full_connectome.d.set(row, column, d);
+                    end
+                    if isprop(result_object, "within_network_pair") && isequal(result_object.within_network_pair, true)
+                        result_object.within_network_pair.d.set(row, column, single_sample_d);
+                    end
                 end
             end
         end

+nla/+net/ResultRank.m

@@ -12,116 +12,152 @@
 
     properties (Dependent)
         permutations
+        number_of_networks
     end
 
     methods
-        function obj = ResultRank(nonpermuted_network_results, permuted_network_results, number_of_network_pairs)
+        function obj = ResultRank(permuted_network_results, number_of_network_pairs)
             if nargin > 0
-                obj.nonpermuted_network_results = nonpermuted_network_results;
+                obj.nonpermuted_network_results = permuted_network_results.no_permutations;
                 obj.permuted_network_results = permuted_network_results;
                 obj.number_of_network_pairs = number_of_network_pairs;
             end
         end
 
         function ranking_result = rank(obj)
+            import nla.TriMatrix nla.TriMatrixDiag
 
             ranking_result = obj.permuted_network_results.copy();
+
+            for test_type = obj.permuted_network_results.test_methods
+                if ~isequal(test_type, "no_permutations") && ~isequal(obj.permuted_network_results.test_display_name, "Cohen's D")
+
+                    [ranking_statistic, probability, denominator] = obj.getTestParameters(test_type);
+                    permutation_results = obj.permuted_network_results.permutation_results;
+                    no_permutation_results = obj.nonpermuted_network_results;
+
+                    % Eggebrecht ranking
+                    ranking_result = obj.eggebrechtRank(test_type, permutation_results, no_permutation_results, ranking_statistic,...
+                        probability, denominator, ranking_result);
 
-            % Ranking
-            if obj.permuted_network_results.test_display_name ~= "Cohen's D"
-               ranking_result = obj.basicRank(ranking_result);
+                    if ~isequal(obj.permuted_network_results.test_name, "hypergeometric")
+                        % Winkler Method ranking
+                        ranking_result.(test_type).winkler_p_value = TriMatrix(...
+                            obj.number_of_networks, TriMatrixDiag.KEEP_DIAGONAL...
+                        );
+                        ranking_result = obj.winklerMethodRank(test_type, permutation_results, no_permutation_results, ranking_statistic,...
+                            probability, denominator, ranking_result);
+
+                        % Westfall Young ranking
+                        ranking_result.(test_type).westfall_young_p_value = TriMatrix(...
+                            obj.number_of_networks, TriMatrixDiag.KEEP_DIAGONAL...
+                        );
+                        ranking_result = obj.westfallYoungMethodRank(test_type, permutation_results, no_permutation_results, ranking_statistic,...
+                            probability, denominator, ranking_result);
+                    end
+                end
             end
         end
 
-        function ranking = basicRank(obj, ranking)
-            ranking_statistic = false;
-            if obj.permuted_network_results.test_display_name ~= "Hypergeometric" % Hypergeomtric has no stat to rank
-                ranking_statistic = obj.permuted_network_results.ranking_statistic;
-            end
-            % Full Connectome ranking
+        function ranking = eggebrechtRank(obj, test_type, permutation_results, no_permutation_results, ranking_statistic,...
+                        probability, denominator, ranking)
 
-            ranking = obj.fullConnectomeRank(ranking, ranking_statistic);
-
-            % Network Pair ranking
-            ranking = obj.withinNetworkPairRank(ranking, ranking_statistic);
-        end
-
-        function ranking = fullConnectomeRank(obj, ranking, ranking_statistic)
-
-            probability = "p_value";
-            no_permutation_result = obj.nonpermuted_network_results.no_permutations;
-            permutation_results = obj.permuted_network_results.permutation_results;
-
-            for index = 1:numel(no_permutation_result.(probability).v)
-                % statistic ranking
-                if obj.permuted_network_results.test_display_name ~= "Hypergeometric"
-                    combined_statistics = [...
-                        permutation_results.(strcat(ranking_statistic, "_permutations")).v(:);...
-                        no_permutation_result.(ranking_statistic).v(index)...
+            for index = 1:numel(no_permutation_results.(probability).v)
+                if isequal(test_type, "full_connectome")
+                    combined_probabilities = [...
+                        permutation_results.(strcat((probability), "_permutations")).v(:);...
+                        no_permutation_results.(probability).v(index)...
+                    ];
+                else
+                    combined_probabilities = [...
+                        permutation_results.(strcat((probability), "_permutations")).v(index, :),...
+                        no_permutation_results.(probability).v(index)...
                     ];
-                    ranking.full_connectome.statistic_p_value.v(index) = sum(...
-                        abs(squeeze(combined_statistics)) >= abs(no_permutation_result.(ranking_statistic).v(index))...
-                        ) / (1 + obj.permutations * obj.number_of_network_pairs);
                 end
-                % p-value ranking
-                combined_probabilities = [...
-                    permutation_results.(strcat(probability, "_permutations")).v(:);...
-                    no_permutation_result.(probability).v(index)...
-                ];
                 [~, sorted_combined_probabilites] = sort(combined_probabilities);
-                ranking.full_connectome.p_value.v(index) = find(...
-                    squeeze(sorted_combined_probabilites) == 1 + obj.permutations * obj.number_of_network_pairs...
-                    ) / (1 + obj.permutations * obj.number_of_network_pairs);
+                ranking.(test_type).(probability).v(index) = find(...
+                    squeeze(sorted_combined_probabilites) == 1 + denominator...
+                    ) / (1 + denominator);
+
+                if ~isequal(obj.permuted_network_results.test_name, "hypergeometric")
+                    combined_statistics = [permutation_results.(strcat((ranking_statistic), "_permutations")).v(:); no_permutation_results.(ranking_statistic).v(index)];
+                    [~, sorted_combined_statistics] = sort(combined_statistics);
+                    ranking.(test_type).(strcat("statistic_", (probability))).v(index) = find(...
+                        squeeze(sorted_combined_statistics) == 1 + denominator...
+                        ) / (1 + denominator);
+                end
             end
-            ranking.full_connectome.d.v = obj.permuted_network_results.full_connectome.d.v;
         end
 
-        function ranking = withinNetworkPairRank(obj, ranking, ranking_statistic)
-
-            if ~any(strcmp(obj.permuted_network_results.test_name, obj.permuted_network_results.noncorrelation_input_tests))
-                single_sample_probability = "single_sample_p_value";
-                single_sample_statistic = strcat("single_sample_", ranking_statistic);
-
-                no_permutation_result = obj.nonpermuted_network_results.no_permutations;
-                permutation_results = obj.permuted_network_results.permutation_results;
-
-                if obj.permuted_network_results.test_name == "wilcoxon"
-                    single_sample_statistic = "single_sample_ranksum_statistic";
-                end
-
-                for index = 1:numel(no_permutation_result.(single_sample_probability).v)
-                    % statistic ranking
-                    combined_statistics = [...
-                        permutation_results.(strcat(single_sample_statistic, "_permutations")).v(index, :),...
-                        no_permutation_result.(single_sample_statistic).v(index)...
-                    ];
-                    ranking.within_network_pair.statistic_single_sample_p_value.v(index) = sum(...
-                        abs(squeeze(combined_statistics)) >= abs(no_permutation_result.(single_sample_statistic).v(index))...
-                        ) / (1 + obj.permutations);
-
-                        % p-value ranking
-                    combined_probabilities = [...
-                        permutation_results.(strcat(single_sample_probability, "_permutations")).v(index, :),...
-                        no_permutation_result.(single_sample_probability).v(index)...
-                    ];
-                    [~, sorted_combined_probabilites] = sort(combined_probabilities);
-                    ranking.within_network_pair.single_sample_p_value.v(index) = find(...
-                        squeeze(sorted_combined_probabilites) == 1 + obj.permutations...
-                        ) / (1 + obj.permutations);
+        function ranking = winklerMethodRank(obj, test_type, permutation_results, no_permutation_results, ranking_statistic,...
+                        probability, denominator, ranking)
+
+            max_statistic_array = max(abs(permutation_results.(strcat(ranking_statistic, "_permutations")).v));
+            for index = 1:numel(no_permutation_results.(probability).v)
+                ranking.(test_type).winkler_p_value.v(index) = sum(...
+                    squeeze(max_statistic_array) >= abs(no_permutation_results.(ranking_statistic).v(index))...
+                );
+            end
+            ranking.(test_type).winkler_p_value.v =...
+                ranking.(test_type).winkler_p_value.v ./ denominator;
+        end
+
+        function ranking = westfallYoungMethodRank(obj, test_type, permutation_results, no_permutation_results, ranking_statistic,...
+                        probability, denominator, ranking)
+
+            % sort statistics in ascending order
+            [sorted_no_permutation_results, sorted_statistic_indexes] = sort(...
+                abs(no_permutation_results.(ranking_statistic).v)...
+            );
+            permutations_sorted_by_non_permuted = abs(...
+                permutation_results.(strcat((ranking_statistic), "_permutations")).v(sorted_statistic_indexes, :)...
+            );
+
+            % Get max value of each permutation starting from max value of non-permuted statistics
+            % Remove each row of permutations associated with non-permuted statistic
+            % Get max of remaining. The last row of permutations should be with the smallest non-permuted statistic
+            % This value gets ranked against its own permutations
+            max_per_permutation_reducing_rows = zeros(...
+                size(permutations_sorted_by_non_permuted, 1), size(permutations_sorted_by_non_permuted, 2)...
+            );
+            for row_index = size(permutations_sorted_by_non_permuted, 1):-1:2
+                max_per_permutation_reducing_rows(row_index, :) = max(permutations_sorted_by_non_permuted(1:row_index, :));
+            end
+            max_per_permutation_reducing_rows(1, :) = permutations_sorted_by_non_permuted(1, :);
+
+            ranking.(test_type).westfall_young.p_value.v = mean(...
+                sorted_no_permutation_results < max_per_permutation_reducing_rows, 2);
+            ranking.(test_type).westfall_young.p_value.v(sorted_statistic_indexes) =...
+                ranking.(test_type).westfall_young.p_value.v;
+        end 
+
+        function [ranking_statistic, probability, denominator] = getTestParameters(obj, test_type)
+
+            ranking_statistic = obj.permuted_network_results.ranking_statistic;
+            probability = "p_value";
+            denominator = obj.permutations * obj.number_of_network_pairs;
+            % Only use these for within network pair and not Chi-Squared and Hypergeometric. 
+            if isequal(test_type, "within_network_pair")  && ~any(...
+                strcmp(obj.permuted_network_results.test_name, obj.permuted_network_results.noncorrelation_input_tests)...
+            )
+                ranking_statistic = strcat("single_sample_", ranking_statistic);
+                if isequal(obj.permuted_network_results.test_name, "wilcoxon")
+                    ranking_statistic = "single_sample_ranksum_statistic"
                 end
-
-            elseif isstruct(obj.permuted_network_results.within_network_pair) &&...
-                any(strcmp(obj.permuted_network_results.test_name, obj.permuted_network_results.noncorrelation_input_tests))
-                % This condition catches Chi-Squared and Hypergeometric tests. We do not do within network ranking for them, we just copy
-                % the full connectome ranking over. 
-                ranking.within_network_pair.single_sample_p_value = ranking.full_connectome.p_value;
+                probability = strcat("single_sample_", probability);
+                denominator = obj.permutations;
             end
-            ranking.within_network_pair.d.v = obj.permuted_network_results.within_network_pair.d.v;
         end
 
+        %% Getters for dependent properties
         % This takes the above statistic and gets the property to use its size to find the number of permutations
         function value = get.permutations(obj)
             value = size(obj.permuted_network_results.permutation_results.p_value_permutations.v, 2); 
         end
+
+        function value = get.number_of_networks(obj)
+            value = obj.permuted_network_results.no_permutations.p_value.size;
+        end
+        %%
     end
 end