diff --git a/+nla/+edge/+permutationMethods/FCWildBootstrap.m b/+nla/+edge/+permutationMethods/FCWildBootstrap.m
index 42999c83..e91ff511 100755
--- a/+nla/+edge/+permutationMethods/FCWildBootstrap.m
+++ b/+nla/+edge/+permutationMethods/FCWildBootstrap.m
@@ -4,7 +4,12 @@
             permuted_input_struct = orig_input_struct;
             permuted_input_struct.func_conn = orig_input_struct.func_conn.copy();
             fcData = permuted_input_struct.func_conn.v';
-            permuted_fcData = nla.helpers.wildBootstrap(fcData, [orig_input_struct.behavior, orig_input_struct.covariates], orig_input_struct.contrasts);
+            if isfield(orig_input_struct,'behavior') && (orig_input_struct.behavior ~= 0)
+                allCovariates = [orig_input_struct.behavior, orig_input_struct.covariates];
+            else
+                allCovariates = orig_input_struct.covariates;
+            end
+            permuted_fcData = nla.helpers.wildBootstrap(fcData, allCovariates, orig_input_struct.contrasts);
             permuted_input_struct.func_conn.v = permuted_fcData';
         end
     end
diff --git a/+nla/+edge/+result/MultiContrast.m b/+nla/+edge/+result/MultiContrast.m
new file mode 100644
index 00000000..8a2f9460
--- /dev/null
+++ b/+nla/+edge/+result/MultiContrast.m
@@ -0,0 +1,55 @@
+classdef MultiContrast < nla.edge.result.Base
+    
+    properties
+        contrastNames %cell array of names of each contrast
+        contrastValues % m x p matrix. Each row is the vector of values for one contrast
+        contrastResultsMap % containers.Map object that stores each contrast with the key of its name
+    end
+    
+    
+    methods
+        
+        function obj = MultiContrast()
+            obj.contrastNames = {};
+            obj.contrastValues = [];
+            obj.contrastResultsMap = containers.Map();
+        end
+        
+        function addNamedResult(obj, contrastName, resultObj)
+            obj.contrastNames{end+1} = contrastName;
+            if isempty(obj.contrastValues)
+                obj.contrastValues = resultObj.contrasts;
+            else
+                obj.contrastValues(end+1,:) = resultObj.contrasts;
+            end
+            obj.contrastResultsMap(contrastName) = resultObj;
+        end
+        
+        function resultObj = getNamedResult(obj, contrastName)
+            if obj.contrastResultsMap.isKey(contrastName)
+                resultObj = obj.contrastResultsMap(contrastName);
+            else
+                resultObj = [];
+            end
+        end
+        
+        function clearResults(obj)
+            obj.contrastNames = {};
+            obj.contrastValues = [];
+            obj.contrastResultsMap = containers.Map();
+        end
+        
+        function outStrArrays = contrastsAsStrings(obj)
+            outStrArrays = obj.contrastResultsMap.keys';
+        end
+        
+        function numberOfContrasts = numContrasts(obj)
+            numberOfContrasts = obj.contrastResultsMap.Count;
+        end
+        
+    end
+    
+    
+    
+    
+end
\ No newline at end of file
diff --git a/+nla/+edge/+test/OrdinaryLeastSquares.m b/+nla/+edge/+test/OrdinaryLeastSquares.m
new file mode 100644
index 00000000..6a920a1d
--- /dev/null
+++ b/+nla/+edge/+test/OrdinaryLeastSquares.m
@@ -0,0 +1,22 @@
+classdef OrdinaryLeastSquares < nla.edge.test.SandwichEstimator
+        
+    
+    methods
+        %Hacky way to change name and coeff_name in superclass is to change
+        %properties in constructor
+        function obj = OrdinaryLeastSquares()
+            
+            obj = obj@nla.edge.test.SandwichEstimator();
+            obj.name = "Ordinary Least Squares";
+            obj.coeff_name = "Contrast T -value";
+        end
+        
+    end
+    
+    methods (Static)        
+        function inputs = requiredInputs()
+            inputs = {nla.inputField.Number('prob_max', 'Edge-level P threshold <', 0, 0.05, 1),...
+                nla.inputField.NetworkAtlasFuncConn(), nla.inputField.OrdinaryLeastSquares()};
+        end
+    end
+end
\ No newline at end of file
diff --git a/+nla/+edge/+test/SandwichEstimator.m b/+nla/+edge/+test/SandwichEstimator.m
new file mode 100755
index 00000000..a16d90b5
--- /dev/null
+++ b/+nla/+edge/+test/SandwichEstimator.m
@@ -0,0 +1,254 @@
+classdef SandwichEstimator < nla.edge.BaseTest
+    %SANDWICHESTIMATOR Summary of this class goes here
+    %   Detailed explanation goes here
+    properties
+        name = "Sandwich Estimator"
+        coeff_name = 'SwE Contrast T-value'
+    end
+    
+    properties
+        % test specific properties go here (things that will persist
+        % over multiple runs) + aren't specific to a given data set)
+        
+        nonpermRegressCoeffs
+        nonpermResiduals
+    end
+    
+    methods
+        function obj = SandwichEstimator()
+            obj@nla.edge.BaseTest();            
+            
+        end
+        
+        function result = run(obj, input_struct)
+            
+            sweInput = obj.sweInputObjFromStruct(input_struct);
+            
+            if ~isfield(input_struct, 'fit_intercept') | input_struct.fit_intercept
+                sweInput = obj.addInterceptCovariateToInputIfNone(sweInput); %Forces model to include fitting an intercept term. Adds column to covariates and contrasts                
+            end    
+            
+            numContrasts = size(sweInput.contrasts,1);
+            if  numContrasts == 1
+                result = obj.fitModel(sweInput);                
+            else
+                result = cell(numContrasts,1);
+                for i = 1:numContrasts
+                    thisInput = sweInput;
+                    thisInput.contrasts = sweInput.contrasts(i,:);
+                    result{i} = obj.fitModel(thisInput);
+                end
+            end
+            
+            
+        end
+        
+        
+    end
+    
+    
+    methods (Access = private)
+        
+        
+        function sweInputStruct = sweInputObjFromStruct(obj, input_struct)
+            
+            sweInputStruct = struct();
+            
+            if isfield(input_struct, 'subjId')
+                sweInputStruct.subjId = input_struct.subjId;
+            else
+                numObs = size(input_struct.covariates,1);
+                sweInputStruct.subjId = (1:numObs)';
+            end
+            
+            if isfield(input_struct, 'groupId')
+                sweInputStruct.groupId = input_struct.groupId;
+            else
+                sweInputStruct.groupId = ones(length(input_struct.subjId),1);
+            end
+            
+            if isfield(input_struct, 'visitId')
+                sweInputStruct.visitId = input_struct.visitId;
+            else
+                sweInputStruct.visitId = ones(length(sweInputStruct.subjId),1);
+            end
+            
+            sweInputStruct.fcData = input_struct.func_conn.v';
+            if isfield(input_struct, 'behavior') & (input_struct.behavior ~= 0)
+                sweInputStruct.covariates = [input_struct.behavior, input_struct.covariates];
+            else
+                sweInputStruct.covariates = input_struct.covariates;
+            end
+            %sweInput.scanMetadata = scanMetadata;
+            sweInputStruct.prob_max = input_struct.prob_max;
+            sweInputStruct.contrasts = input_struct.contrasts;
+            
+            if isfield(input_struct, 'stdErrCalcObj')
+                sweInputStruct.stdErrCalcObj = input_struct.stdErrCalcObj; %How will this really be passed in?
+            else
+                sweInputStruct.stdErrCalcObj = nla.helpers.stdError.UnconstrainedBlocks();
+            end
+            
+            if isfield(input_struct, 'fit_intercept')
+                sweInputStruct.fit_intercept = input_struct.fit_intercept;
+            end
+            
+        end
+        
+        function sweInput = addInterceptCovariateToInputIfNone(obj, sweInput)
+            
+            
+            columnIsAllSameValue = ~any(diff(sweInput.covariates,1),1);
+            if ~any(columnIsAllSameValue)
+                %If no column in covariates is currently all same value,
+                %add a column so that an intercept term is fit by the
+                %linear model
+                numObs = size(sweInput.covariates,1);
+                sweInput.covariates = [sweInput.covariates, ones(numObs,1)];
+                sweInput.contrasts = [sweInput.contrasts, 0];
+                
+            end
+                
+            
+            
+        end
+        
+        
+        function sweRes = fitModel(obj, input)
+                    
+                        
+            numFcEdges = size(input.fcData,2);
+            
+            %the data for each scan in fcData (ie each row) represents the
+            %flattened lower triangle of fc data. Use the number of fc
+            %edges represented to compute the size of the original
+            %non-flattened matrix, which is needed to construct a result
+            %object
+            fcMatSize = (1 + sqrt(1 + 4*(numFcEdges*2))) / 2;
+            
+            sweRes = nla.edge.result.SandwichEstimator(fcMatSize, input.prob_max);
+            
+            designMtx = input.covariates;
+            %designMtx = obj.zeroMeanUnitVarianceByColumn(designMtx); %Make covariates zero mean unit variance
+                               
+            [regressCoeffs, residual] = obj.fitLinearModel(input.fcData, designMtx);
+                        
+            
+            %Build input and pass to one of several methods for calculating
+            %standard error
+            %TODO: Should this calculate beta covariance instead?
+            
+            stdErrCalcObj = input.stdErrCalcObj;
+            %stdErrCalcObj = nlaEckDev.sweStdError.Guillaume(); %If you want to hard code the std err calc object
+            
+            stdErrInput = struct();% nlaEckDev.sweStdError.SwEStdErrorInput();
+            scanMetadata = struct();
+            scanMetadata.subjId = input.subjId; % vec [numScans]
+            scanMetadata.groupId = input.groupId;%groupId; % vec [numScans]        
+            scanMetadata.visitId = input.visitId; % vec [numScans] 
+            
+            stdErrInput.scanMetadata = scanMetadata;
+            stdErrInput.residual = residual;
+            stdErrInput.pinvDesignMtx = pinv(designMtx);
+            
+            %sweRes.stdError = stdErrCalcObj.calculate(stdErrInput);
+            stdError = stdErrCalcObj.calculate(stdErrInput);
+            
+            
+            contrastCalc = input.contrasts * regressCoeffs;
+            contrastSE = sqrt((input.contrasts.^2) * (stdError.^2));
+            
+            dof = obj.calcDegreesOfFreedom(designMtx);
+            
+            %tVals = regressCoeffs ./ stdError;
+            contrastTVal = contrastCalc ./ contrastSE;
+            
+            %pVals = zeros(size(tVals));  
+            contrastPVal = 2*(1-cdf('T',abs(contrastTVal),dof));
+                
+            %sweRes.tVals = tVals;
+            %sweRes.pVals = pVals;
+            %sweRes.regressCoeffs = regressCoeffs;             
+            sweRes.coeff.v = contrastTVal';
+            sweRes.prob.v = contrastPVal';
+            sweRes.prob_sig.v = (sweRes.prob.v < input.prob_max);            
+            sweRes.avg_prob_sig = sum(sweRes.prob_sig.v) ./ numel(sweRes.prob_sig.v);
+            sweRes.contrasts = input.contrasts;
+            
+            %Change expected coefficient range to be more accurate to Sandwich
+            %Estimator ranges
+            sweRes.coeff_range = [-3 3];
+            
+        end
+        
+        
+        function outResidual = calcResidual(obj, X, pInvX, Y)
+            
+            hat = X * pInvX;
+            
+            residCorrectFactor = (1 - diag(hat)).^(-1); %Type 3 residual correction defined in Guillaume 2014
+            
+            %compute residuals
+            
+            % Goal is to multiply row for each subject by respective
+            % correction factor in vector residCorrectFactor.
+            % Below equation is equivalent to row wise multiplying Y-hat*Y
+            % with elements of residCorrectFactor.
+            % It is exactly equivalent to:
+            %    outResidual = diag(residCorrectFactor) * (Y - hat*Y);
+            % but over twice as fast by avoiding constructing sparse
+            % [subjxsubj] matrix of diag(residCorrectFactor)   
+            
+            %outResidual = diag(residCorrectFactor) * (Y - hat*Y);
+            outResidual = residCorrectFactor .* (Y - hat*Y);
+
+            
+        end
+       
+        function outMtx = zeroMeanUnitVarianceByColumn(obj, inMtx)
+            
+            outMtx = zeros(size(inMtx));
+            
+            for colIdx = 1:size(inMtx,2)
+                
+                colMean = mean(inMtx(:,colIdx));
+                colStd = std(inMtx(:,colIdx));
+                if colStd ~=0
+                    outMtx(:,colIdx) = (inMtx(:,colIdx) - colMean) / colStd;
+                else
+                    %If all values of column are same, do not change them
+                    %TODO: handle this here??? or somewhere else?
+                    %TODO: set to zeros? ones? orig values?
+                    outMtx(:,colIdx) = inMtx(:,colIdx);
+                end
+                
+            end
+            
+        end
+        
+        function degOfFree = calcDegreesOfFreedom(obj, designMtx)
+            
+            degOfFree = size(designMtx,1) - size(designMtx,2) - 1;
+            
+        end
+        
+        function [regressCoeffs, residual] = fitLinearModel(obj, fcData, designMtx)
+            
+            pinvDesignMtx = pinv(designMtx);
+            regressCoeffs = pinvDesignMtx * fcData;            
+            residual = obj.calcResidual(designMtx, pinvDesignMtx, fcData);
+            
+        end        
+        
+        
+        
+    end % end private methods
+    
+    methods (Static)
+        function inputs = requiredInputs()
+            inputs = {nla.inputField.Number('prob_max', 'Edge-level P threshold <', 0, 0.05, 1),...
+                nla.inputField.NetworkAtlasFuncConn(), nla.inputField.SandwichEstimator()};
+        end
+    end
+end
+
diff --git a/+nla/+edge/ContrastInput.m b/+nla/+edge/ContrastInput.m
new file mode 100644
index 00000000..f839522b
--- /dev/null
+++ b/+nla/+edge/ContrastInput.m
@@ -0,0 +1,32 @@
+classdef ContrastInput < handle & matlab.mixin.Copyable
+    
+    properties
+        dataTable = [];
+        name = '';
+        contrastVector = [];        
+    end
+    
+    methods
+        function isValidFlag = isValid(obj)
+            if any(~isspace(obj.name))
+                isValidFlag = true;
+            else
+                isValidFlag = false;
+            end
+        end
+        
+        function str = asDetailedString(obj)        
+            str = sprintf('%s (%s)',obj.name, obj.contrastVectorAsString());
+        end
+        
+        function str = contrastVectorAsString(obj)
+            str = num2str(obj.contrastVector);
+            str = regexprep(str, '\s+', ' ');
+        end
+    end
+    
+    
+    
+        
+    
+end
\ No newline at end of file
diff --git a/+nla/+edge/ContrastInputGUI.mlapp b/+nla/+edge/ContrastInputGUI.mlapp
new file mode 100644
index 00000000..30125ad5
Binary files /dev/null and b/+nla/+edge/ContrastInputGUI.mlapp differ
diff --git a/+nla/+edge/SandwichEstimatorInputGUI.mlapp b/+nla/+edge/SandwichEstimatorInputGUI.mlapp
new file mode 100644
index 00000000..27e8cdd6
Binary files /dev/null and b/+nla/+edge/SandwichEstimatorInputGUI.mlapp differ
diff --git a/+nla/+helpers/+stdError/Guillaume.m b/+nla/+helpers/+stdError/Guillaume.m
index f2f65cc3..616f4501 100755
--- a/+nla/+helpers/+stdError/Guillaume.m
+++ b/+nla/+helpers/+stdError/Guillaume.m
@@ -9,6 +9,9 @@
     properties
         totalSpectralCorrections;
     end
+    properties (SetAccess = protected)
+        REQUIRES_GROUP = true;
+    end
     
     methods
         
diff --git a/+nla/+helpers/+stdError/Heteroskedastic.m b/+nla/+helpers/+stdError/Heteroskedastic.m
index 0305632a..2dc7174b 100755
--- a/+nla/+helpers/+stdError/Heteroskedastic.m
+++ b/+nla/+helpers/+stdError/Heteroskedastic.m
@@ -1,10 +1,14 @@
 classdef Heteroskedastic < nla.helpers.stdError.AbstractSwEStdErrStrategy
+    
+    properties (SetAccess = protected)
+        REQUIRES_GROUP = true;
+    end
 
     methods
         
         function stdErr = calculate(obj, sweStdErrInput)
             
-            validateattributes(sweStdErrInput, 'nlaEckDev.sweStdError.SwEStdErrorInput', {});
+            
             %Calculation of standard error assuming heteroskedascticity
             %consistent errors
             [numCovariates, numObs] = size(sweStdErrInput.pinvDesignMtx);
diff --git a/+nla/+helpers/+stdError/Heteroskedastic_FAST.m b/+nla/+helpers/+stdError/Heteroskedastic_FAST.m
index 5b240894..ef672a86 100755
--- a/+nla/+helpers/+stdError/Heteroskedastic_FAST.m
+++ b/+nla/+helpers/+stdError/Heteroskedastic_FAST.m
@@ -1,4 +1,8 @@
 classdef Heteroskedastic_FAST < nla.helpers.stdError.AbstractSwEStdErrStrategy
+    
+    properties (SetAccess = protected)
+        REQUIRES_GROUP = true;
+    end
 
     methods
         
diff --git a/+nla/+helpers/+stdError/Homoskedastic.m b/+nla/+helpers/+stdError/Homoskedastic.m
index fce89f6b..0c7b0c1b 100755
--- a/+nla/+helpers/+stdError/Homoskedastic.m
+++ b/+nla/+helpers/+stdError/Homoskedastic.m
@@ -1,4 +1,8 @@
 classdef Homoskedastic < nla.helpers.stdError.AbstractSwEStdErrStrategy
+    
+    properties (SetAccess = protected)
+        REQUIRES_GROUP = false;
+    end
 
     methods
         
diff --git a/+nla/+helpers/+stdError/UnconstrainedBlocks.m b/+nla/+helpers/+stdError/UnconstrainedBlocks.m
index 18944e06..682848e9 100755
--- a/+nla/+helpers/+stdError/UnconstrainedBlocks.m
+++ b/+nla/+helpers/+stdError/UnconstrainedBlocks.m
@@ -5,6 +5,9 @@
         SPARSITY_THRESHOLD = 0.2;
         
     end
+    properties (SetAccess = protected)
+        REQUIRES_GROUP = true;
+    end
     
     methods
         
@@ -27,11 +30,11 @@
             FORCE_HALF_SW_ALGO = true;
             
             if FORCE_HALF_SW_ALGO
-                stdErrStrategy = nlaEckDev.sweStdError.UnconstrainedBlocks_BenKay();            
+                stdErrStrategy = nla.helpers.stdError.UnconstrainedBlocks_BenKay();            
             elseif vSparsity <= obj.SPARSITY_THRESHOLD
-                stdErrStrategy = nlaEckDev.sweStdError.UnconstrainedBlocks_Sparse();
+                stdErrStrategy = nla.helpers.stdError.UnconstrainedBlocks_Sparse();
             else
-                stdErrStrategy = nlaEckDev.sweStdError.UnconstrainedBlocks_Dense();
+                stdErrStrategy = nla.helpers.stdError.UnconstrainedBlocks_Dense();
             end
             
             stdErr = stdErrStrategy.calculate(sweStdErrInput);
diff --git a/+nla/+helpers/wildBootstrap.m b/+nla/+helpers/wildBootstrap.m
index d366a5e3..fb964a35 100755
--- a/+nla/+helpers/wildBootstrap.m
+++ b/+nla/+helpers/wildBootstrap.m
@@ -1,33 +1,51 @@
-function outY = wildBootstrap(inY, inX, inContrasts)
+function outY = wildBootstrap(inY, inX, inContrasts, origBeta, origResidual)
     
     %Inputs
     %inY - [n x m] matrix (n observations, m outputs per observation)
     %inX - [n x p] matrix (n observations, p covariates)
     %inContrasts - 1 x p matrix indicating which covariates we care about,
-    %and which are noise. Model and resulting residulas will be calculated
-    %by fitting a model only to noise
+    %and which are noise. Model and resulting residuals will be calculated
+    %by fitting a model only to noise    
+    %origBeta (optional) - p x 1 vector of beta weights if we've already fit the
+    %linear model
+    %origResidual (optional) - n x 1 vector of residuals if we've already fit the
+    %linear model
 
-    %Fit linear model to only noise variables
-    covariateIsNoise = inContrasts == 0;    
-    xNoise = inX(:,covariateIsNoise);    
-    pinvNoiseX = pinv(xNoise);    
-    beta = pinvNoiseX * inY;    
-    residual = inY - xNoise * beta;
     
-
-    %TODO: make this flexible to choose different random dists to
-    %hit residuals with
-    %Currently is hardcoded to use rademacher distribution (ie +1
-    %or -1 with equal likelihood
+    %If beta and residual not provided, calculate them now
+    
+    if nargin <=3
+        
+        %Fit linear model to only noise variables - incorrect! fit to all, then
+        %re-add betas for noise variables
+    %     covariateIsNoise = inContrasts == 0;    
+    %     xNoise = inX(:,covariateIsNoise);    
+    %     pinvNoiseX = pinv(xNoise);    
+    %     beta = pinvNoiseX * inY;    
+    %     residual = inY - xNoise * beta;
+
+        %Fit linear model to all variables
+        pinvX = pinv(inX);
+        beta = pinvX * inY;
+    else
+        beta = origBeta;
+    end
+    
+    %if residual not provided, calculate it from x, y, and beta
+    if nargin <=4
+        residual = inY - inX * beta;
+    else
+        residual = origResidual;
+    end
+    
+    %generate new values by using betas from noise only, and adding in
+    %random multiplier of residual
     numObs = size(inX,1);
     residualNoiseMultFactor = getRademacherVector(numObs);
-
-    %Applying this noise factor to all fc edges for one
-    %observation. Is this valid???
-
-    outY = xNoise * beta + (residualNoiseMultFactor .* residual); %MATLAB can interpret dot multiplication as row-wise multiplication!
-
-
+    
+    covariateIsNoise = inContrasts == 0;
+    outY = inX(:,covariateIsNoise) * beta(covariateIsNoise,:) + (residualNoiseMultFactor .* residual);
+    
 end
     
 function outVec = getRademacherVector(numPts)
diff --git a/+nla/+inputField/OrdinaryLeastSquares.m b/+nla/+inputField/OrdinaryLeastSquares.m
new file mode 100644
index 00000000..f664d8c2
--- /dev/null
+++ b/+nla/+inputField/OrdinaryLeastSquares.m
@@ -0,0 +1,57 @@
+classdef OrdinaryLeastSquares < nla.inputField.SandwichEstimator
+    %This class will act nearly identically to the SandwichEstimator input,
+    %with the only difference being that it will startup the sandwich
+    %estimator input GUI with an extra flag indicating that we will force
+    %the standard error calculation to be the homoskedastic one.
+    %
+    %It overrides the following functions from the base class:
+    %draw(...), giving the button a different label specific to OLS, and
+    %a different callback that starts the sandwich input GUI with a flag
+    %that forces OLS
+    
+    methods
+    
+        function obj = OrdinaryLeastSquares()
+            obj = obj@nla.inputField.SandwichEstimator();
+            obj.name = 'ordinary_least_squares';
+            obj.disp_name = 'Ordinary Least Squares';
+        end
+    end
+    
+    methods
+        function [w, h] = draw(obj, x, y, parent, fig)
+            import nla.inputField.LABEL_H nla.inputField.LABEL_GAP
+            
+            table_h = 300;
+            h = LABEL_H + LABEL_GAP + table_h + LABEL_GAP + LABEL_H + LABEL_GAP + LABEL_H;
+            h = 50;
+            
+            
+            %% 'Set Behavior' button
+            [obj.button_start_sandwich_gui, w3] = obj.createButton(obj.button_start_sandwich_gui, 'OLS Settings', parent, x,...
+                y - h + LABEL_H + LABEL_GAP + LABEL_H, @(h,e)obj.button_startSandwichGUI_withForceOLS_Callback());
+           
+            w = 0;
+            
+        end
+        
+    end
+    
+    methods (Access = protected)
+        
+        function button_startSandwichGUI_withForceOLS_Callback(obj, ~)
+
+            obj.settingsObj.FORCE_ORDINARY_LEAST_SQUARES = true;
+
+            settingsGUI = nla.edge.SandwichEstimatorInputGUI(obj.settingsObj);
+            uiwait(settingsGUI.SandwichEstimatorInputUIFigure);  
+            if obj.settingsObj.isValid()
+                obj.satisfied = true;
+            end
+
+        end
+    
+    end
+    
+    
+end
\ No newline at end of file
diff --git a/+nla/+inputField/SandwichEstimator.m b/+nla/+inputField/SandwichEstimator.m
new file mode 100644
index 00000000..1ffc313f
--- /dev/null
+++ b/+nla/+inputField/SandwichEstimator.m
@@ -0,0 +1,99 @@
+classdef SandwichEstimator < nla.inputField.InputField
+    properties
+        name = 'sandwich_estimator';
+        disp_name = 'Sandwich Estimator';
+    end
+    
+    properties
+        settingsObj
+    end
+    
+    properties (Access = protected)
+        button_start_sandwich_gui = false        
+    end
+    
+    
+    methods
+        function obj = SandwichEstimator()
+            obj.settingsObj = nla.edge.inputHelpers.SandwichEstimatorInput();
+        end
+        
+        function [w, h] = draw(obj, x, y, parent, fig)
+            import nla.inputField.LABEL_H nla.inputField.LABEL_GAP
+            
+            table_h = 300;
+            h = LABEL_H + LABEL_GAP + table_h + LABEL_GAP + LABEL_H + LABEL_GAP + LABEL_H;
+            h = 50;
+            
+            
+            %% 'Set Behavior' button
+            [obj.button_start_sandwich_gui, w3] = obj.createButton(obj.button_start_sandwich_gui, 'SWE Settings', parent, x,...
+                y - h + LABEL_H + LABEL_GAP + LABEL_H, @(h,e)obj.button_startSandwichGUICallback());
+           
+            w = 0;
+            
+        end
+        
+        function undraw(obj)
+            
+            if isgraphics(obj.button_start_sandwich_gui)
+                delete(obj.button_start_sandwich_gui)
+            end
+        end
+        
+        function read(obj, input_struct)
+            obj.settingsObj = nla.edge.inputHelpers.SandwichEstimatorInput();
+            
+            obj.settingsObj.loadFieldsFromStruct(input_struct);
+            
+            
+        end
+        
+        function [input_struct, error] = store(obj, input_struct)
+            
+            input_struct = obj.settingsObj.pushFieldsIntoStruct(input_struct);
+            error = false;
+        end
+    end
+    
+    methods (Access = protected)
+        function [button, w] = createButton(obj, button, label, parent, x, y, callback)
+            
+            %% Create button
+            if ~isgraphics(button)
+                button = uibutton(parent, 'push', 'ButtonPushedFcn', callback);
+            end
+            button_w = 100;
+            button.Position = [x, y - nla.inputField.LABEL_H, button_w, nla.inputField.LABEL_H];
+            
+            button.Text = label;
+            button.Position(3) = nla.inputField.widthOfString(button.Text, nla.inputField.LABEL_H) +...
+                nla.inputField.widthOfString('  ', nla.inputField.LABEL_H + nla.inputField.LABEL_GAP);
+            
+            w = button.Position(3);
+        end
+        
+        
+        function button_startSandwichGUICallback(obj, ~)
+                        
+            obj.settingsObj.FORCE_ORDINARY_LEAST_SQUARES = false;
+            settingsGUI = nla.edge.SandwichEstimatorInputGUI(obj.settingsObj);
+            uiwait(settingsGUI.SandwichEstimatorInputUIFigure);  
+            if obj.settingsObj.isValid()
+                obj.satisfied = true;
+            end
+            
+        end
+        
+        function update(obj)
+        end
+        
+        function settingsStruct = getSettingsObjAsStruct(obj)
+            warning off;
+            settingsStruct = struct(obj.settingsObj);
+            warning on;
+            
+        end
+    end
+end
+
diff --git a/+nla/+net/+mcc/Base.m b/+nla/+net/+mcc/Base.m
index ca4fceea..7cade860 100755
--- a/+nla/+net/+mcc/Base.m
+++ b/+nla/+net/+mcc/Base.m
@@ -1,6 +1,6 @@
 classdef Base
     methods (Abstract)
-        p_max = correct(obj, net_atlas, input_struct, prob)
+        [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
         correction_label = createLabel(obj, net_atlas, input_struct, prob)
     end
 end
\ No newline at end of file
diff --git a/+nla/+net/+mcc/BenjaminiHochberg.m b/+nla/+net/+mcc/BenjaminiHochberg.m
index a879373f..23063cad 100755
--- a/+nla/+net/+mcc/BenjaminiHochberg.m
+++ b/+nla/+net/+mcc/BenjaminiHochberg.m
@@ -4,21 +4,11 @@
     end
     
     methods
-        function p_max = correct(obj, net_atlas, input_struct, prob)
-            [~, p_max] = nla.lib.fdr_bh(prob.v, input_struct.prob_max, 'pdep');
+        function [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
+            [is_sig_vector, p_max] = nla.lib.fdr_bh(prob.v, input_struct.prob_max, 'pdep');
         end
         function correction_label = createLabel(obj, net_atlas, input_struct, prob)
-            p_max = obj.correct(net_atlas, input_struct, prob);
-            if p_max == 0
-                correction_label = sprintf('FDR_{BH} produced no significant nets');
-            else
-                format_specs = "%g/%d tests";
-                if isequal(input_struct.behavior_count, 1)
-                    format_specs = "%g/%d test";
-                end
-                correction_label = sprintf(strcat("FDR_{BH}(", format_specs, ")"), input_struct.prob_max * input_struct.behavior_count,...
-                    input_struct.behavior_count);
-            end
+            correction_label = sprintf("Benjamini-Hochberg (alpha = %g)",input_struct.prob_max);
         end
     end
 end
\ No newline at end of file
diff --git a/+nla/+net/+mcc/BenjaminiYekutieli.m b/+nla/+net/+mcc/BenjaminiYekutieli.m
index 5c8fb6f9..23e0de8d 100755
--- a/+nla/+net/+mcc/BenjaminiYekutieli.m
+++ b/+nla/+net/+mcc/BenjaminiYekutieli.m
@@ -4,21 +4,11 @@
     end
     
     methods
-        function p_max = correct(obj, net_atlas, input_struct, prob)
-            [~, p_max] = nla.lib.fdr_bh(prob.v, input_struct.prob_max, 'dep');
+        function [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
+            [is_sig_vector, p_max] = nla.lib.fdr_bh(prob.v, input_struct.prob_max, 'dep');
         end
         function correction_label = createLabel(obj, net_atlas, input_struct, prob)
-            p_max = obj.correct(net_atlas, input_struct, prob);
-            if p_max == 0
-                correction_label = sprintf('FDR_{BY} produced no significant nets');
-            else
-                format_specs = "%g/%d tests";
-                if isequal(input_struct.behavior_count, 1)
-                    format_specs = "%g/%d test";
-                end
-                correction_label = sprintf(strcat("FDR_{BY}(", format_specs, ")"), input_struct.prob_max * input_struct.behavior_count,...
-                    input_struct.behavior_count);
-            end
+            correction_label = sprintf("Benjamini-Yekutieli (alpha = %g)",input_struct.prob_max);
         end
     end
 end
\ No newline at end of file
diff --git a/+nla/+net/+mcc/Bonferroni.m b/+nla/+net/+mcc/Bonferroni.m
index de306236..8897b413 100755
--- a/+nla/+net/+mcc/Bonferroni.m
+++ b/+nla/+net/+mcc/Bonferroni.m
@@ -4,15 +4,18 @@
     end
     
     methods
-        function p_max = correct(obj, net_atlas, input_struct, prob)
+        function [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
+            
             p_max = input_struct.prob_max / net_atlas.numNetPairs();
+            is_sig_vector = prob.v < p_max;
         end
         function correction_label = createLabel(obj, net_atlas, input_struct, prob)
-            format_specs_tests = "%d tests";
+            format_specs_tests = "%d tests)";
             if isequal(input_struct.behavior_count, 1)
-                format_specs_tests = "%d test";
+                format_specs_tests = "%d test)";
             end
-            correction_label = sprintf(strcat("%g/%d net-pairs/", format_specs_tests), input_struct.prob_max * input_struct.behavior_count,...
+            p_max = input_struct.prob_max / net_atlas.numNetPairs();
+            correction_label = sprintf(strcat("P < %.2g (%.2g/%d net-pairs/", format_specs_tests), p_max, input_struct.prob_max * input_struct.behavior_count,...
                 net_atlas.numNetPairs(), input_struct.behavior_count);
         end
     end
diff --git a/+nla/+net/+mcc/HolmBonferroni.m b/+nla/+net/+mcc/HolmBonferroni.m
index 08efeda7..2f8c3995 100644
--- a/+nla/+net/+mcc/HolmBonferroni.m
+++ b/+nla/+net/+mcc/HolmBonferroni.m
@@ -4,22 +4,16 @@
     end
 
     methods
-        function p_max = correct(obj, net_atlas, input_struct, prob)
-            [is_significant, adjusted_pvals, ~] = nla.lib.bonferroni_holm(prob.v, input_struct.prob_max);
-            p_max = max(is_significant .* adjusted_pvals);
+        function [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
+            [is_sig_vector, adjusted_pvals, ~] = nla.lib.bonferroni_holm(prob.v, input_struct.prob_max);
+                        
+            p_max = max(is_sig_vector .* prob.v);
         end
         function correction_label = createLabel(obj, net_atlas, input_struct, prob)
-            p_max = obj.correct(net_atlas, input_struct, prob);
-            if p_max == 0
-                correction_label = sprintf('Holm-Bonferroni produced no significant networks');
-            else
-                format_specs = "%g/%d tests";
-                if isequal(input_struct.behavior_count, 1)
-                    format_specs = "%g/%d test";
-                end
-                correction_label = sprintf(strcat("Holm-Bonferroni(", format_specs, ")"), input_struct.prob_max * input_struct.behavior_count,...
-                    input_struct.behavior_count);
-            end
+            correction_label = sprintf("Holm-Bonferroni (alpha = %g)",input_struct.prob_max);
+            
+            %Since p threshold is variable with this test, exclude it and
+            %just use the initial 'alpha' term used in the algorithm.
         end
     end
 end
diff --git a/+nla/+net/+mcc/None.m b/+nla/+net/+mcc/None.m
index 8ef9dd66..f1bba3df 100755
--- a/+nla/+net/+mcc/None.m
+++ b/+nla/+net/+mcc/None.m
@@ -4,15 +4,16 @@
     end
     
     methods
-        function p_max = correct(obj, net_atlas, input_struct, prob)
+        function [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)            
             p_max = input_struct.prob_max;
+            is_sig_vector = prob.v < p_max;
         end
         function correction_label = createLabel(obj, net_atlas, input_struct, prob)
-            format_specs = "%g/%d tests";
+            format_specs = "P < %.2g (%g/%d tests)";
             if isequal(input_struct.behavior_count, 1)
-                format_specs = "%g/%d test";
+                format_specs = "P < %.2g (%g/%d test)";
             end
-            correction_label = sprintf(format_specs, input_struct.prob_max * input_struct.behavior_count,...
+            correction_label = sprintf(format_specs, input_struct.prob_max, input_struct.prob_max * input_struct.behavior_count,...
                 input_struct.behavior_count);
         end
     end
diff --git a/+nla/+net/+result/+plot/NetworkTestPlot.m b/+nla/+net/+result/+plot/NetworkTestPlot.m
index 4e00a3af..b3626707 100644
--- a/+nla/+net/+result/+plot/NetworkTestPlot.m
+++ b/+nla/+net/+result/+plot/NetworkTestPlot.m
@@ -89,8 +89,8 @@ function getPlotTitle(obj)
                 obj.title = sprintf("%s (D > %g)", obj.title, obj.network_test_options.d_max);
             end
             if ~isequal(obj.test_method, "no_permutations")
-                if isequal(obj.current_settings.ranking, nla.RankingMethod.WINKLER) 
-                    obj.title = strcat(obj.title, "\nRanking by Winkler Method");
+                if isequal(obj.current_settings.ranking, nla.RankingMethod.FREEDMAN_LANE) 
+                    obj.title = strcat(obj.title, "\nRanking by Freedman-Lane Method");
                 elseif isequal(obj.current_settings.ranking, nla.RankingMethod.WESTFALL_YOUNG)
                     obj.title = strcat(obj.title, "\nRanking by Westfall-Young Method");
                 else
@@ -228,8 +228,8 @@ function resizeFigure(obj, plot_width, plot_height)
             % All the options (buttons, pulldowns, checkboxes)
             scale_option = PullDown("plot_scale", "Plot Scale", ["Linear", "Log", "Negative Log10"],...
                 [ProbPlotMethod.DEFAULT, ProbPlotMethod.LOG, ProbPlotMethod.NEGATIVE_LOG_10]);
-            ranking_method = PullDown("ranking", "Ranking", ["Uncorrected", "Winkler", "Westfall-Young"],...
-                [RankingMethod.UNCORRECTED, RankingMethod.WINKLER, RankingMethod.WESTFALL_YOUNG]);
+            ranking_method = PullDown("ranking", "Ranking", ["Uncorrected", "Freedman-Lane", "Westfall-Young"],...
+                [RankingMethod.UNCORRECTED, RankingMethod.FREEDMAN_LANE, RankingMethod.WESTFALL_YOUNG]);
             cohens_d = CheckBox("cohens_d", "Cohen's D Threshold", true);
             centroids = CheckBox("centroids", "ROI Centroids in brain plots", false);
             multiple_comparison_correction = PullDown("mcc", "Multiple Comparison Correction",...
diff --git a/+nla/+net/+result/+plot/NetworkTestPlotApp.mlapp b/+nla/+net/+result/+plot/NetworkTestPlotApp.mlapp
index d2513df1..ea13c24c 100644
Binary files a/+nla/+net/+result/+plot/NetworkTestPlotApp.mlapp and b/+nla/+net/+result/+plot/NetworkTestPlotApp.mlapp differ
diff --git a/+nla/+net/+result/+plot/NetworkTestPlotApp_exported.m b/+nla/+net/+result/+plot/NetworkTestPlotApp_exported.m
index 88dfd8c6..c7fba833 100644
--- a/+nla/+net/+result/+plot/NetworkTestPlotApp_exported.m
+++ b/+nla/+net/+result/+plot/NetworkTestPlotApp_exported.m
@@ -88,8 +88,8 @@ function getPlotTitle(app)
                 app.title = sprintf("%s (D > %g)", app.title, app.CohensDThresholdEditField.Value);
             end
             if ~isequal(app.test_method, "no_permutations")
-                if isequal(app.RankingDropDown.Value, "nla.RankingMethod.WINKLER") % Look at me, I'm MATLAB. I have no idea why enums are beneficial or how to use them
-                    app.title = strcat(app.title, "\nRanking by Winkler Method");
+                if isequal(app.RankingDropDown.Value, "nla.RankingMethod.FREEDMAN_LANE") % Look at me, I'm MATLAB. I have no idea why enums are beneficial or how to use them
+                    app.title = strcat(app.title, "\nRanking by Freedman-Lane Method");
                 elseif isequal(app.RankingDropDown.Value, "nla.RankingMethod.WESTFALL_YOUNG")
                     app.title = strcat(app.title, "\nRanking by Westfall-Young Method");
                 else
@@ -439,8 +439,8 @@ function createComponents(app)
 
             % Create RankingDropDown
             app.RankingDropDown = uidropdown(app.Panel);
-            app.RankingDropDown.Items = {'Uncorrected', 'Winkler', 'Westfall-Young'};
-            app.RankingDropDown.ItemsData = {'nla.RankingMethod.UNCORRECTED', 'nla.RankingMethod.WINKLER', 'nla.RankingMethod.WESTFALL_YOUNG'};
+            app.RankingDropDown.Items = {'Uncorrected', 'Freedman-Lane', 'Westfall-Young'};
+            app.RankingDropDown.ItemsData = {'nla.RankingMethod.UNCORRECTED', 'nla.RankingMethod.FREEDMAN_LANE', 'nla.RankingMethod.WESTFALL_YOUNG'};
             app.RankingDropDown.ValueChangedFcn = createCallbackFcn(app, @PlotScaleValueChanged, true);
             app.RankingDropDown.Position = [291 297 100 22];
             app.RankingDropDown.Value = 'nla.RankingMethod.UNCORRECTED';
diff --git a/+nla/+net/+result/NetworkResultPlotParameter.m b/+nla/+net/+result/NetworkResultPlotParameter.m
index a72e8ad4..c2c81e85 100644
--- a/+nla/+net/+result/NetworkResultPlotParameter.m
+++ b/+nla/+net/+result/NetworkResultPlotParameter.m
@@ -27,7 +27,7 @@
         end
 
         function result = plotProbabilityParameters(obj, edge_test_options, edge_test_result, test_method, plot_statistic,...
-                plot_title, fdr_correction, significance_filter, ranking_method)
+                plot_title, fdr_correction, effect_size_filter, ranking_method)
             % plot_title - this will be a string
             % plot_statistic - this is the stat that will be plotted
             % significance filter - this will be a boolean or some sort of object (like Cohen's D > D-value)
@@ -36,9 +36,9 @@
 
             import nla.TriMatrix nla.TriMatrixDiag
             % We're going to use a default filter here
-            if isequal(significance_filter, false)
-                significance_filter = TriMatrix(obj.number_of_networks, "logical", TriMatrixDiag.KEEP_DIAGONAL);
-                significance_filter.v = true(numel(significance_filter.v), 1);
+            if isequal(effect_size_filter, false)
+                effect_size_filter = TriMatrix(obj.number_of_networks, "logical", TriMatrixDiag.KEEP_DIAGONAL);
+                effect_size_filter.v = true(numel(effect_size_filter.v), 1);
             end
 
             % Adding on to the plot title if it's a -log10 plot
@@ -53,20 +53,26 @@
             if isstring(fdr_correction) || ischar(fdr_correction)
                 fdr_correction = nla.net.mcc.(erase(fdr_correction, "-"))();
             end
-            p_value_max = fdr_correction.correct(obj.network_atlas, obj.updated_test_options, statistic_input);
+            
+            [is_sig_vector, p_value_max] = fdr_correction.correct(obj.network_atlas, obj.updated_test_options, statistic_input);
+            
             p_value_breakdown_label = fdr_correction.createLabel(obj.network_atlas, obj.updated_test_options,...
                 statistic_input);
-
-            name_label = sprintf("%s %s\nP < %.2g (%s)", obj.network_test_results.test_display_name, plot_title,...
-                p_value_max, p_value_breakdown_label);
-            if p_value_max == 0
-                name_label = sprintf("%s %s\nP = %.2g (%s)", obj.network_test_results.test_display_name, plot_title,...
-                    p_value_max, p_value_breakdown_label);
-            end
+            name_label = sprintf("%s %s\n%s", obj.network_test_results.test_display_name, plot_title,...
+                p_value_breakdown_label);
+%             else
+% 
+%                 name_label = sprintf("%s %s\nP < %.2g (%s)", obj.network_test_results.test_display_name, plot_title,...
+%                     p_value_max, p_value_breakdown_label);
+%                 if p_value_max == 0
+%                     name_label = sprintf("%s %s\nP = %.2g (%s)", obj.network_test_results.test_display_name, plot_title,...
+%                         p_value_max, p_value_breakdown_label);
+%                 end
+%             end
 
             % Filtering if there's a filter provided 
             significance_plot = TriMatrix(obj.number_of_networks, "logical", TriMatrixDiag.KEEP_DIAGONAL);
-            significance_plot.v = (statistic_input.v < p_value_max) & significance_filter.v;
+            significance_plot.v = is_sig_vector & effect_size_filter.v;
 
             % scale values very slightly for display so numbers just below
             % the threshold don't show up white but marked significant
@@ -195,8 +201,8 @@ function brainFigureButtonCallback(network1, network2)
             end
 
             switch ranking_method
-                case "nla.RankingMethod.WINKLER"
-                    ranking = "winkler_";
+                case "nla.RankingMethod.FREEDMAN_LANE"
+                    ranking = "freedman_lane_";
                 case "nla.RankingMethod.WESTFALL_YOUNG"
                     ranking = "westfall_young_";
                 otherwise
diff --git a/+nla/+net/+result/NetworkTestResult.m b/+nla/+net/+result/NetworkTestResult.m
index e8ff9792..a4068827 100644
--- a/+nla/+net/+result/NetworkTestResult.m
+++ b/+nla/+net/+result/NetworkTestResult.m
@@ -242,7 +242,7 @@ function createPValueTriMatrices(obj, number_of_networks, test_method)
 
             non_correlation_test = any(strcmp(obj.test_name, obj.noncorrelation_input_tests));
             uncorrected_names = ["uncorrected_", "legacy_"];
-            corrected_names = ["winkler_", "westfall_young_"];
+            corrected_names = ["freedman_lane_", "westfall_young_"];
 
             switch test_method
                 case "no_permutations"
@@ -285,15 +285,15 @@ function createPValueTriMatrices(obj, number_of_networks, test_method)
         % I don't really know what these do and haven't really thought about it. Hence the bad naming.
         function [sig, name] = singleSigMat(obj, network_atlas, edge_test_options, p_value, mcc_method, title_prefix)
             mcc_method = nla.net.mcc.(mcc_method)();
-            p_value_max = mcc_method.correct(network_atlas, edge_test_options, p_value);
+            [is_sig_vector, p_value_max] = mcc_method.correct(network_atlas, edge_test_options, p_value);
             p_breakdown_labels = mcc_method.createLabel(network_atlas, edge_test_options, p_value);
 
             sig = nla.TriMatrix(network_atlas.numNets(), 'double', nla.TriMatrixDiag.KEEP_DIAGONAL);
-            sig.v = (p_value.v < p_value_max);
-            name = sprintf("%s %s P < %.2g (%s)", title_prefix, obj.test_display_name, p_value_max, p_breakdown_labels);
-            if p_value_max == 0
-                name = sprintf("%s %s P = 0 (%s)", title_prefix, obj.test_display_name, p_breakdown_labels);
-            end
+            sig.v = is_sig_vector;
+            name = sprintf("%s %s %s", title_prefix, obj.test_display_name, p_value_max, p_breakdown_labels);
+%             if p_value_max == 0
+%                 name = sprintf("%s %s P = 0 (%s)", title_prefix, obj.test_display_name, p_breakdown_labels);
+%             end
         end
 
         function [number_of_tests, sig_count_mat, names] = appendSignificanceMatrix(...
@@ -337,11 +337,11 @@ function createPValueTriMatrices(obj, number_of_networks, test_method)
             % :return: The full name of the p-value. (example: "single_sammple_p_value")
 
             import nla.NetworkLevelMethod
-            noncorrelation_input_tests = ["chi_squared", "hypergeometric"];
-            non_correlation_test = any(strcmp(test_name, noncorrelation_input_tests));
+            reference_distribution_tests = ["chi_squared", "hypergeometric"];
+            reference_distribution_test = any(strcmp(test_name, reference_distribution_tests));
 
             probability = "two_sample_p_value";
-            if isequal(non_correlation_test, false)
+            if isequal(reference_distribution_test, false)
                 if isequal(test_method, "no_permutations") || isequal(test_method, "within_network_pair")
                     probability = "single_sample_p_value";
                 end
diff --git a/+nla/+net/ResultRank.m b/+nla/+net/ResultRank.m
index 4ab3594b..c9cd1788 100644
--- a/+nla/+net/ResultRank.m
+++ b/+nla/+net/ResultRank.m
@@ -42,8 +42,8 @@
                     ranking_result = obj.uncorrectedRank(test_method, permutation_results, no_permutation_results, ranking_statistic,...
                         probability, ranking_result);
 
-                    % Winkler Method ranking
-                    ranking_result = obj.winklerMethodRank(test_method, permutation_results, no_permutation_results, ranking_statistic,...
+                    % Freedman-Lane Method ranking
+                    ranking_result = obj.freedmanLaneMethodRank(test_method, permutation_results, no_permutation_results, ranking_statistic,...
                         probability, ranking_result);
 
                     % Westfall Young ranking
@@ -96,7 +96,7 @@
         end
 
 
-        function ranking = winklerMethodRank(obj, test_method, permutation_results, no_permutation_results, ranking_statistic,...
+        function ranking = freedmanLaneMethodRank(obj, test_method, permutation_results, no_permutation_results, ranking_statistic,...
             probability, ranking)
             % Ranks the observed result using method described by Winkler to correct for FWER
             %
@@ -108,7 +108,7 @@
             % :param ranking: The NetworkTestResult object to place the results
             % :return: The same NetworkTestResult object with ranking results
 
-            winkler_probability = strcat("winkler_", probability);
+            freedman_lane_probability = strcat("freedman_lane_", probability);
             % NET-278: The "max statistic" for hypergeometric is the most significant p-value. Which is the smallest p-value.
             max_statistic_array = max(abs(permutation_results.(strcat(ranking_statistic, "_permutations")).v));
             if isequal(ranking.test_name, "hypergeometric")
@@ -117,16 +117,16 @@
             
             for index = 1:numel(no_permutation_results.(strcat("uncorrected_", probability)).v)
                 if isequal(ranking.test_name, "hypergeometric")
-                    ranking.(test_method).(winkler_probability).v(index) = sum(...
+                    ranking.(test_method).(freedman_lane_probability).v(index) = sum(...
                         squeeze(max_statistic_array) <= abs(no_permutation_results.(ranking_statistic).v(index))...
                     );
                 else
-                    ranking.(test_method).(winkler_probability).v(index) = sum(...
+                    ranking.(test_method).(freedman_lane_probability).v(index) = sum(...
                         squeeze(max_statistic_array) >= abs(no_permutation_results.(ranking_statistic).v(index))...
                     );
                 end
             end
-            ranking.(test_method).(winkler_probability).v = ranking.(test_method).(winkler_probability).v ./ obj.permutations;
+            ranking.(test_method).(freedman_lane_probability).v = ranking.(test_method).(freedman_lane_probability).v ./ obj.permutations;
         end
 
         function ranking = westfallYoungMethodRank(obj, test_method, permutation_results, no_permutation_results, ranking_statistic,...
diff --git a/+nla/+net/unittests/ResultRankTestCase.m b/+nla/+net/unittests/ResultRankTestCase.m
index 0f65c4dd..8f75ffe4 100644
--- a/+nla/+net/unittests/ResultRankTestCase.m
+++ b/+nla/+net/unittests/ResultRankTestCase.m
@@ -114,11 +114,11 @@ function fullConnectomeUncorrectedRankingTest(testCase)
             testCase.verifyEqual(uncorrected_results, expected_uncorrected);
         end
 
-        function fullConnectomeWinklerRankingTest(testCase)
-            winkler_results = testCase.rank.full_connectome.winkler_two_sample_p_value.v;
-            expected_winkler = testCase.rank_results.full_connectome.winkler_two_sample_p_value.v;
+        function fullConnectomeFreedmanLaneRankingTest(testCase)
+            freedman_lane_results = testCase.rank.full_connectome.freedman_lane_two_sample_p_value.v;
+            expected_freedman_lane = testCase.rank_results.full_connectome.freedman_lane_two_sample_p_value.v;
 
-            testCase.verifyEqual(winkler_results, expected_winkler);
+            testCase.verifyEqual(freedman_lane_results, expected_freedman_lane);
         end
 
         function fullConnectomeWestfallYoungRankingTest(testCase)
@@ -140,11 +140,11 @@ function withinNetworkPairUncorrectedRankingTest(testCase)
             testCase.verifyEqual(uncorrected_results, expected_uncorrected);           
         end
 
-        function withinNetworkPairWinklerRankingTest(testCase)
-            winkler_results = testCase.rank.within_network_pair.winkler_single_sample_p_value.v;
-            expected_winkler = testCase.rank_results.within_network_pair.winkler_single_sample_p_value.v;
+        function withinNetworkPairFreedmanLaneRankingTest(testCase)
+            freedman_lane_results = testCase.rank.within_network_pair.freedman_lane_single_sample_p_value.v;
+            expected_freedman_lane = testCase.rank_results.within_network_pair.freedman_lane_single_sample_p_value.v;
 
-            testCase.verifyEqual(winkler_results, expected_winkler);
+            testCase.verifyEqual(freedman_lane_results, expected_freedman_lane);
         end
 
         function withinNetworkPairWestfallYoungRankingTest(testCase)
diff --git a/+nla/MultiContrastResult.m b/+nla/MultiContrastResult.m
new file mode 100644
index 00000000..de825be1
--- /dev/null
+++ b/+nla/MultiContrastResult.m
@@ -0,0 +1,67 @@
+classdef MultiContrastResult < matlab.mixin.Copyable
+    %Class for holding multiple result objects for multiple named contrasts
+    %
+    % Builds and accesses a containers.Map object of named objects
+    % objects, one per contrast
+    
+    properties
+        contrastNames %cell array of names of each contrast
+        contrastResultsMap % containers.Map object that stores each contrast with the key of its name
+    end
+    
+    
+    methods
+        
+        function obj = MultiContrastResult()
+            obj.contrastNames = {};
+            obj.contrastResultsMap = containers.Map();
+        end
+        
+        function addNamedResult(obj, contrastName, resultObj)
+            obj.contrastNames{end+1} = contrastName;
+            obj.contrastResultsMap(contrastName) = resultObj;
+        end
+        
+        function resultObj = getNamedResult(obj, contrastName)
+            if obj.contrastResultsMap.isKey(contrastName)
+                resultObj = obj.contrastResultsMap(contrastName);
+            else
+                resultObj = [];
+            end
+        end
+        
+        function clearResults(obj)
+            obj.contrastNames = {};
+            obj.contrastResultsMap = containers.Map();
+        end
+        
+        function outStrArrays = contrastsAsStrings(obj)
+            outStrArrays = obj.contrastResultsMap.keys';
+        end
+        
+        function numberOfContrasts = numContrasts(obj)
+            numberOfContrasts = obj.contrastResultsMap.Count;
+        end
+        
+        function merge(obj, next_multi_contrast_result)
+            %append results from a different multi contrast result to this
+            %one
+            contrast_names = next_multi_contrast_result.contrastsAsStrings();
+            for contrastIdx = 1:length(contrast_names)
+                this_contrast_name = contrast_names{contrastIdx};
+                this_contrast_results = obj.getNamedResult(this_contrast_name);
+                next_contrast_result = next_multi_contrast_result.getNamedResult(this_contrast_name);
+                for i = 1:length(this_contrast_results)
+                    this_contrast_results{i}.merge(next_contrast_result{i});
+                end
+                
+            end
+            
+        end
+        
+    end
+    
+    
+    
+    
+end
\ No newline at end of file
diff --git a/+nla/RankingMethod.m b/+nla/RankingMethod.m
index cb3f98b1..978f90ac 100644
--- a/+nla/RankingMethod.m
+++ b/+nla/RankingMethod.m
@@ -1,5 +1,5 @@
 classdef RankingMethod
     enumeration
-        UNCORRECTED, WINKLER, WESTFALL_YOUNG
+        UNCORRECTED, FREEDMAN_LANE, WESTFALL_YOUNG
     end
 end
\ No newline at end of file
diff --git a/+nla/TestPool.m b/+nla/TestPool.m
index fdca9536..862747bf 100755
--- a/+nla/TestPool.m
+++ b/+nla/TestPool.m
@@ -36,31 +36,88 @@
         
         function result = runPerm(obj, edge_input_struct, net_input_struct, network_atlas, nonpermuted_edge_test_results,...
             nonpermuted_network_test_results, num_perms, perm_seed, separate_network_and_edge_tests)
-            
+        
             if ~exist('perm_seed', 'var')
                 perm_seed = false;
             end
-            
+
             if ~exist('separate_network_and_edge_tests', 'var')
                 separate_network_and_edge_tests = false;
-            end
+            end                
+        
+%             if isa(nonpermuted_edge_test_results,'nla.edge.result.MultiContrast')
+%                 %If result has multiple contrasts, call the wrapper
+%                 %function that handles running them one at a time
+%                 result = obj.runPermMultiContrast(edge_input_struct, net_input_struct, network_atlas, nonpermuted_edge_test_results,...
+%                     nonpermuted_network_test_results, num_perms, perm_seed, separate_network_and_edge_tests);
+%             else
 
-            if isequal(separate_network_and_edge_tests, false)
-                [permuted_edge_test_results, permuted_network_test_results] = obj.runEdgeAndNetPerm(edge_input_struct,...
-                    net_input_struct, network_atlas, nonpermuted_edge_test_results, num_perms, perm_seed);
-            else
-                [permuted_edge_test_results, permuted_network_test_results] = obj.runPermSeparateEdgeAndNet(edge_input_struct,...
-                    net_input_struct, network_atlas, num_perms, perm_seed);
-            end
-            
-            ranked_permuted_network_test_results = obj.collateNetworkPermutationResults(nonpermuted_edge_test_results, network_atlas,...
-                nonpermuted_network_test_results, permuted_network_test_results, net_input_struct);
+                if isequal(separate_network_and_edge_tests, false)
+                    [permuted_edge_test_results, permuted_network_test_results] = obj.runEdgeAndNetPerm(edge_input_struct,...
+                        net_input_struct, network_atlas, nonpermuted_edge_test_results, num_perms, perm_seed);
+                else
+                    [permuted_edge_test_results, permuted_network_test_results] = obj.runPermSeparateEdgeAndNet(edge_input_struct,...
+                        net_input_struct, network_atlas, num_perms, perm_seed);
+                end
+
+                if ~isa(permuted_network_test_results, 'nla.MultiContrastResult')
+                    ranked_permuted_network_test_results = obj.collateNetworkPermutationResults(nonpermuted_edge_test_results, network_atlas,...
+                        nonpermuted_network_test_results, permuted_network_test_results, net_input_struct);
 
-            result = nla.ResultPool(edge_input_struct, net_input_struct, network_atlas, nonpermuted_edge_test_results,...
-                nonpermuted_network_test_results, permuted_edge_test_results, ranked_permuted_network_test_results);
+                    result = nla.ResultPool(edge_input_struct, net_input_struct, network_atlas, nonpermuted_edge_test_results,...
+                        nonpermuted_network_test_results, permuted_edge_test_results, ranked_permuted_network_test_results);
+                else
+                    
+                    result = nla.MultiContrastResult();
+                    contrastNames = nonpermuted_edge_test_results.contrastNames;
+
+                    for i = 1:length(contrastNames)
+                        this_contrast_name = contrastNames{i};
+
+                        this_nonperm_edge_result = nonpermuted_edge_test_results.getNamedResult(this_contrast_name);
+                        this_nonperm_network_result = nonpermuted_network_test_results.getNamedResult(this_contrast_name);
+                        this_perm_network_result = permuted_network_test_results.getNamedResult(this_contrast_name);
+
+                        this_ranked_permuted_network_test_results = obj.collateNetworkPermutationResults(this_nonperm_edge_result, network_atlas,...
+                                            this_nonperm_network_result, this_perm_network_result, net_input_struct);
+
+                        this_result_pool = nla.ResultPool(edge_input_struct, net_input_struct, network_atlas, this_nonperm_edge_result,...
+                                            this_nonperm_network_result, this_perm_network_result, this_ranked_permuted_network_test_results);
+
+                        result.addNamedResult(this_contrast_name, this_result_pool);
+                    end
+                    
+                    
+                end
+            
+            %end
             
         end
         
+        function multi_contrast_result = runPermMultiContrast(obj, edge_input_struct, net_input_struct, network_atlas, nonpermuted_edge_test_results,...
+            nonpermuted_network_test_results, num_perms, perm_seed, separate_network_and_edge_tests)
+            %If the input is multi contrast, this will generate result pools for the edge and
+            %net results one at a time via loop
+            
+            contrast_names = nonpermuted_edge_test_results.contrastsAsStrings();
+            num_contrasts = length(contrast_names);
+            
+            multi_contrast_result = nla.MultiContrastResult();
+            
+            for i = 1:num_contrasts
+                this_contrast_name = contrast_names{i};
+                this_nonperm_edge_result = nonpermuted_edge_test_results.getNamedResult(this_contrast_name);
+                this_nonperm_net_results = nonpermuted_network_test_results.getNamedResult(this_contrast_name);
+                
+                this_contrast_result_pool = obj.runPerm(edge_input_struct, net_input_struct, network_atlas, this_nonperm_edge_result,...
+                    this_nonperm_net_results, num_perms, perm_seed, separate_network_and_edge_tests);
+                
+                multi_contrast_result.addNamedResult(this_contrast_result_pool);
+                
+            end
+        
+        end
+        
         function ranked_results = collateNetworkPermutationResults(obj, nonpermuted_edge_test_results, network_atlas, nonpermuted_network_test_results,...
             permuted_network_test_results, network_test_options)
             
@@ -111,6 +168,7 @@
             % get current parallel pool or start a new one
             [number_of_processes, blocks] = obj.initializeParallelPool(num_perms);
             
+            
             parfor process = 1:number_of_processes
                 network_result_block = obj.runEdgeAndNetPermBlock(edge_input_struct, net_input_struct, net_atlas,...
                     blocks(process), blocks(process+1), perm_seed);
@@ -123,9 +181,13 @@
             permuted_network_results = network_result_blocks{1};
             for process = 2:number_of_processes
                 current_network_test_results = network_result_blocks{process};
-                for test_index = 1:numNetTests(obj)
-                    current_test_network_result = current_network_test_results(test_index);
-                    permuted_network_results{test_index}.merge(current_test_network_result);
+                if ~isa(current_network_test_results, 'nla.MultiContrastResult')
+                    for test_index = 1:numNetTests(obj)
+                        current_test_network_result = current_network_test_results(test_index);
+                        permuted_network_results{test_index}.merge(current_test_network_result);
+                    end
+                else
+                    permuted_network_results.merge(current_network_test_results);
                 end
             end
         end
@@ -143,13 +205,21 @@
                 
                 % Ugh, this is so horrible. Have to do this due to Matlab not being able to index 2D arrays separately among
                 % indexes
-                if iteration - block_start + 1 == 1
-                    for test = 1:numNetTests(obj)
-                        network_result_block{test} = copy(network_results{test});
+                if isa(network_results, 'nla.MultiContrastResult')
+                    if iteration - block_start + 1 == 1
+                        network_result_block = network_results;
+                    else
+                        network_result_block.merge(network_results);
                     end
-                else
-                    for test = 1:numNetTests(obj)
-                        network_result_block{test}.merge(network_results{test});
+                else                    
+                    if iteration - block_start + 1 == 1
+                        for test = 1:numNetTests(obj)
+                            network_result_block{test} = copy(network_results{test});
+                        end
+                    else
+                        for test = 1:numNetTests(obj)
+                            network_result_block{test}.merge(network_results{test});
+                        end
                     end
                 end
 
@@ -214,7 +284,31 @@
                 input_struct.iteration = 0;
             end
             
-            edge_result = obj.edge_test.run(input_struct);
+            is_multicontrast = isfield(input_struct, 'contrasts') && (size(input_struct.contrasts,1)>1);
+            
+            if ~is_multicontrast
+                edge_result = obj.edge_test.run(input_struct);                  
+            else
+                %Generate edge results for all contrasts
+                num_contrasts = size(input_struct.contrasts,1);
+                edge_result = nla.edge.result.MultiContrast();
+                has_contrast_names = isfield(input_struct,'contrastNames') & ~isempty(input_struct.contrastNames);
+                for contrast_idx = 1:num_contrasts
+                    
+                    this_contrast = input_struct.contrasts(contrast_idx,:);
+                    if has_contrast_names
+                        this_contrast_name = input_struct.contrastNames{contrast_idx};
+                    else
+                        this_contrast_name = ['contrast', num2str(contrast_idx)];
+                    end
+                    
+                    input_struct_this_contrast = input_struct;
+                    input_struct_this_contrast.contrasts = this_contrast;
+                    input_struct_this_contrast.contrastNames = this_contrast_name;
+                    edge_result_this_contrast = obj.edge_test.run(input_struct_this_contrast);
+                    edge_result.addNamedResult(this_contrast_name, edge_result_this_contrast);
+                end
+            end
         end
 
         function net_level_results = runNetTestsPerm(obj, net_input_struct, net_atlas, perm_edge_results)
@@ -267,9 +361,25 @@
         end
         
         function net_results = runNetTests(obj, input_struct, edge_result, net_atlas, permutations)
-            net_results = {};
-            for i = 1:numNetTests(obj)
-                net_results{i} = obj.net_tests{i}.run(input_struct, edge_result, net_atlas, permutations);
+            
+            if ~isa(edge_result, 'nla.edge.result.MultiContrast')                              
+                net_results = {};
+                for i = 1:numNetTests(obj)
+                    net_results{i} = obj.net_tests{i}.run(input_struct, edge_result, net_atlas, permutations);
+                end            
+            else
+                net_results = nla.MultiContrastResult();
+                contrast_names = edge_result.contrastsAsStrings();
+                num_contrasts = length(contrast_names);
+                for contrastIdx = 1:num_contrasts
+                    this_contrast_name = contrast_names{contrastIdx};
+                    this_contrast_edge_result = edge_result.getNamedResult(this_contrast_name);
+                    this_contrast_net_results = {};
+                    for i = 1:numNetTests(obj)
+                        this_contrast_net_result{i} = obj.net_tests{i}.run(input_struct, this_contrast_edge_result, net_atlas, permutations);
+                    end
+                	net_results.addNamedResult(this_contrast_name,this_contrast_net_result)
+                end
             end
         end
         
diff --git a/NLAResult.mlapp b/NLAResult.mlapp
index c21cb8cf..c5023f33 100644
Binary files a/NLAResult.mlapp and b/NLAResult.mlapp differ
diff --git a/NLAResult_exported.m b/NLAResult_exported.m
index 62e9bb1d..05be62dd 100644
--- a/NLAResult_exported.m
+++ b/NLAResult_exported.m
@@ -6,15 +6,17 @@
         FileMenu                   matlab.ui.container.Menu
         SaveButton                 matlab.ui.container.Menu
         SaveSummaryTableMenu       matlab.ui.container.Menu
-        OpenDiagnosticPlotsButton  matlab.ui.control.Button
-        OpenTriMatrixPlotButton    matlab.ui.control.Button
-        BranchLabel                matlab.ui.control.Label
-        RunButton                  matlab.ui.control.Button
-        NetLevelLabel              matlab.ui.control.Label
-        ViewEdgeLevelButton        matlab.ui.control.Button
-        EdgeLevelLabel             matlab.ui.control.Label
-        FlipNestingButton          matlab.ui.control.Button
         ResultTree                 matlab.ui.container.Tree
+        FlipNestingButton          matlab.ui.control.Button
+        EdgeLevelLabel             matlab.ui.control.Label
+        ViewEdgeLevelButton        matlab.ui.control.Button
+        NetLevelLabel              matlab.ui.control.Label
+        RunButton                  matlab.ui.control.Button
+        BranchLabel                matlab.ui.control.Label
+        OpenTriMatrixPlotButton    matlab.ui.control.Button
+        OpenDiagnosticPlotsButton  matlab.ui.control.Button
+        SelectContrastLabel        matlab.ui.control.Label
+        SelectContrastDropdown     matlab.ui.control.DropDown
     end
 
     
@@ -29,6 +31,10 @@
         net_adjustable_fields
         cur_iter = 0
         old_data = false % data from old result objects, not NetworkTestResult
+        
+        multi_contrast_edge_result = false
+        multi_contrast_net_result = false
+        selected_contrast_name = ''
     end
     
     methods (Access = private)
@@ -174,6 +180,17 @@ function initFromInputs(app, test_pool, input_struct, net_input_struct, ~)
             
             app.edge_result = test_pool.runEdgeTest(input_struct);
             
+            if app.edgeResultIsMultiContrast(app.edge_result)
+                app.multi_contrast_edge_result = app.edge_result;
+                app.initializeContrastSelectorFromResult(app.multi_contrast_edge_result);
+                
+                %set edge result to the one corresponding to the first
+                %contrast
+                all_contrast_names = app.multi_contrast_edge_result.contrastsAsStrings();
+                app.edge_result = app.multi_contrast_edge_result.getNamedResult(all_contrast_names{1});                
+                
+            end
+            
             app.input_struct = input_struct;
             app.net_input_struct = net_input_struct;
             app.test_pool = test_pool;
@@ -212,6 +229,17 @@ function initFromResult(app, result, file_name, ~, old_data)
                 app.net_input_struct.prob_max = app.net_input_struct.prob_max_original;
             end
             
+            if app.edgeResultIsMultiContrast(app.edge_result)
+                app.multi_contrast_edge_result = app.edge_result;
+                app.initializeContrastSelectorFromResult(app.multi_contrast_edge_result);
+                
+                %set edge result to the one corresponding to the first
+                %contrast
+                all_contrast_names = app.multi_contrast_edge_result.contrastsAsStrings();
+                app.edge_result = app.multi_contrast_edge_result.getNamedResult(all_contrast_names{1});                
+                
+            end
+            
             app.results = result;
             app.old_data = old_data;
        
@@ -233,6 +261,37 @@ function initFromResult(app, result, file_name, ~, old_data)
             
         end
         
+        function initializeContrastSelectorFromResult(app, edge_result)
+            if app.edgeResultIsMultiContrast(edge_result)
+                app.SelectContrastDropdown.Items = edge_result.contrastsAsStrings();
+                app.SelectContrastDropdown.Enable = true;
+                app.SelectContrastDropdown.Value = app.SelectContrastDropdown.Items{1};
+                app.SelectContrastLabel.Enable = true;
+            else
+                app.SelectContrastDropdown.Enable = false;
+                app.SelectContrastLabel.Enable = false;
+                app.SelectContrastDropdown.Items = {};
+                app.SelectContrastDropdown.Value = {};
+            end
+        
+        end
+        
+        function is_multi_contrast = edgeResultIsMultiContrast(app, edge_result)
+            if isa(edge_result, 'nla.edge.result.MultiContrast')
+                is_multi_contrast = true;
+            else
+                is_multi_contrast = false;
+            end
+        end
+        
+        function is_multi_contrast = resultPoolIsMultiContrast(app, result_pool)
+            if isa(result_pool, 'nla.MultiContrastResult')
+                is_multi_contrast = true;
+            else
+                is_multi_contrast = false;
+            end        
+        end
+        
         function enableNetButtons(app, val)
             % buttons that need net-level data to be used
             net_buttons = {app.ResultTree, app.FlipNestingButton, app.SaveSummaryTableMenu};
@@ -300,7 +359,19 @@ function RunButtonPushed(app, event)
             prog.Value = 0.02;
             drawnow;
             
-            net_results = app.test_pool.runNetTests(app.net_input_struct, app.edge_result, app.input_struct.net_atlas, false);
+            if app.multi_contrast_edge_result == false
+                edge_result = app.edge_result;
+            else
+                %If multiple contrast result exists, we'll
+                %run all of them, rather than just the currently selected
+                %one
+                edge_result = app.multi_contrast_edge_result;
+            end
+                
+                
+            
+            net_results = app.test_pool.runNetTests(app.net_input_struct, edge_result, app.input_struct.net_atlas, false);
+            
             
             if app.net_input_struct.full_connectome
                 prog.Message = 'Starting parallel pool...';
@@ -322,7 +393,16 @@ function RunButtonPushed(app, event)
 
                 % Run permuted statistics
                 app.test_pool.data_queue = data_queue;
-                app.results = app.test_pool.runPerm(app.input_struct, app.net_input_struct, app.input_struct.net_atlas, app.edge_result, net_results, app.net_input_struct.perm_count);
+                resultPool = app.test_pool.runPerm(app.input_struct, app.net_input_struct, app.input_struct.net_atlas, edge_result, net_results, app.net_input_struct.perm_count);
+                
+                if ~app.resultPoolIsMultiContrast(resultPool)
+                    app.results = resultPool;
+                else
+                    app.multi_contrast_net_result = resultPool;
+                    contrastName = app.SelectContrastDropdown.Value;
+                    app.results = app.multi_contrast_net_result.getNamedResult(contrastName);
+                end
+                
                 
                 % Delete the data queue
                 delete(data_queue);
@@ -330,7 +410,7 @@ function RunButtonPushed(app, event)
                 % Unset handle reference
                 app.prog_bar = false;
             else
-                app.results = ResultPool(app.input_struct, app.net_input_struct, app.input_struct.net_atlas, app.edge_result, net_results, false, false);
+                app.results = ResultPool(app.input_struct, app.net_input_struct, app.input_struct.net_atlas, edge_result, net_results, false, false);
             end
             
             prog.Value = 0.98;
@@ -561,6 +641,16 @@ function CohensDthresholdchordplotsCheckBoxValueChanged(app, event)
             d_thresh = app.CohensDthresholdchordplotsCheckBox.Value;
             app.net_input_struct.d_thresh_chord_plot = d_thresh;
         end
+
+        % Value changed function: SelectContrastDropdown
+        function SelectContrastDropdownValueChanged(app, event)
+            app.selected_contrast_name = app.SelectContrastDropdown.Value;
+            app.edge_result = app.multi_contrast_edge_result.getNamedResult(app.selected_contrast_name);
+            if app.resultPoolIsMultiContrast(app.multi_contrast_net_result)
+                app.results = app.multi_contrast_net_result.getNamedResult(app.selected_contrast_name);
+                app.setNesting();
+            end
+        end
     end
 
     % Component initialization
@@ -643,6 +733,20 @@ function createComponents(app)
             app.OpenDiagnosticPlotsButton.Position = [11 28 147 22];
             app.OpenDiagnosticPlotsButton.Text = 'Open Diagnostic Plots';
 
+            % Create SelectContrastLabel
+            app.SelectContrastLabel = uilabel(app.UIFigure);
+            app.SelectContrastLabel.Enable = 'off';
+            app.SelectContrastLabel.Position = [143 633 152 22];
+            app.SelectContrastLabel.Text = 'Select Contrast';
+
+            % Create SelectContrastDropdown
+            app.SelectContrastDropdown = uidropdown(app.UIFigure);
+            app.SelectContrastDropdown.Items = {};
+            app.SelectContrastDropdown.ValueChangedFcn = createCallbackFcn(app, @SelectContrastDropdownValueChanged, true);
+            app.SelectContrastDropdown.Enable = 'off';
+            app.SelectContrastDropdown.Position = [141 610 154 22];
+            app.SelectContrastDropdown.Value = {};
+
             % Show the figure after all components are created
             app.UIFigure.Visible = 'on';
         end
diff --git a/NLA_GUI.mlapp b/NLA_GUI.mlapp
index 36853e50..fd6e8f7e 100755
Binary files a/NLA_GUI.mlapp and b/NLA_GUI.mlapp differ
diff --git a/docs/source/methodology.rst b/docs/source/methodology.rst
index 2ea5220d..28026972 100644
--- a/docs/source/methodology.rst
+++ b/docs/source/methodology.rst
@@ -45,10 +45,6 @@ Benjamini-Hochberg, Westfall and Young :cite:p:`WestfallP`, and Freedman-Lane :c
     functional connectivity data are not shuffled in order to preserve the inherent covariant structure of the
     data across permutations
 
-.. note::
-    Freedman-Lane FWER correction is referred to as Winkler in NLA. Implementation was modeled after the algorithm described in the paper
-    by Winkler.
-
 Brain Network Map Selection
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
diff --git a/docs/source/network_level_results.rst b/docs/source/network_level_results.rst
index edc91855..17c43fc7 100644
--- a/docs/source/network_level_results.rst
+++ b/docs/source/network_level_results.rst
@@ -49,7 +49,7 @@ number of points above or below (depending on the test used) the non-permuted (o
 In NLA this is referred to as "ranking" since it is simply counting values in a sorted list. This basic
 ranking is referred to as the "uncorrected" *p*-value. There are two other options for ranking in NLA. These
 account for :abbr:`FWER (family-wise error rate)`. The first method is based off the "randomise" method :cite:p:`FreedmanD,WinklerA`.
-This is referred to as the "Winkler method". The second method is called "Westfall-Young" in NLA described by
+The second method is called "Westfall-Young" in NLA described by
 an alogrithm :cite:p:`WestfallP` by Westfall and Young.
 
 Result Rank
@@ -66,7 +66,7 @@ Result Rank
 
     .. mat:automethod:: uncorrectedRank
 
-    .. mat:automethod:: winklerMethodRank
+    .. mat:automethod:: freedmanLaneMethodRank
 
     .. mat:automethod:: westfallYoungMethodRank