first commit

BootstrapInteractions.m → DGSA_computations/BootstrapInteractions.m

@@ -40,7 +40,7 @@
 
 for j = 1:NbClusters
         idx_c =  find(Clustering.T == j); % find points in the cluster
-
+        q_prior = quantile(ParametersValues(idx_c,:),(1:1:99)./100,1);% distribution of parameter p in the entire cluster: F(p|c)
          % Bin the conditional parameter and store the indices in a vector
          for l = 1:NbBins
             if l == 1
@@ -55,12 +55,19 @@
                 BootInteractions(:,j,l) = NaN(NbParams,1);
             else
                 for i =1:NbParams 
-                    q_prior = quantile(ParametersValues(idx_c,i),(1:1:99)./100);% distribution of parameter p in the entire cluster: F(p|c)
                     boot = zeros(NbDraw,1);
+                    ParametersValuesi = ParametersValues(:,i);
+                    x_redraw = zeros(length(idx_cl),NbDraw);
+                    for iter = 1:NbDraw
+                        x_redraw(:,iter) = ParametersValuesi(idx_c(randsample(Clustering.weights(j),length(idx_cl),0)));  % sample points
+                    end
+
+                    %x_redraw = reshape(ParametersValuesi(idx_c(randsample(Clustering.weights(j),length(idx_cl)*NbDraw,'true'))),[],NbDraw);  % sample points
+                    q = quantile(x_redraw,(1:1:99)./100,1);  % bootstrapped distribution
                     for iter = 1:NbDraw % resampling procedure
-                        x_redraw = idx_c(randsample(1:Clustering.weights(j),length(idx_cl))'); % sampling from F(p|c)
-                        q = quantile(ParametersValues(x_redraw,i),(1:1:99)./100); % distribution of parameter p conditioned to pj in bin l: F(p|i(pj,tl),c)
-                        boot(iter) = norm(q_prior-q,1); % L1-norm
+                        %x_redraw = idx_c(randsample(1:Clustering.weights(j),length(idx_cl))'); % sampling from F(p|c)
+                        %q = quantile(ParametersValues(x_redraw,i),(1:1:99)./100); % distribution of parameter p conditioned to pj in bin l: F(p|i(pj,tl),c)
+                        boot(iter) = norm(q_prior(:,i)-q(:,iter),1); % L1-norm
                     end
                     BootInteractions(i,j,l) = quantile(boot,alpha);
                 end
@@ -69,4 +76,4 @@
 end
 BootInteractions(CondIdx,:,:) = [];
 
-end
+end

BootstrapMainFactors.m → DGSA_computations/BootstrapMainFactors.m

@@ -5,6 +5,9 @@
 % Author: Celine Scheidt
 % Date: August 2012
 
+% Later added by Jihoon Park (jhpark3@stanford.edu)
+% Date: 14 May 2016
+
 function BootMainFactors = BootstrapMainFactors(ParametersValues,Clustering,NbDraw,alpha)
 
 %% Input Parameters
@@ -24,18 +27,26 @@
         alpha = 0.95;
     end
 
-    BootMainFactors = zeros(nbParametersValues,nbclusters);
+%% Modified by Jihoon to evalualte multiple vaules of alpha
+
+    BootMainFactors = zeros(nbParametersValues,nbclusters,length(alpha));
+
+%% Resampling
 
     for i = 1:nbParametersValues
         q_prior = quantile(ParametersValues(:,i),(1:1:99)./100);  % prior distribution
         for  j = 1:nbclusters
             boot = zeros(NbDraw,1);
             for iter = 1:NbDraw
-                x_redraw = randsample(size(ParametersValues,1),Clustering.weights(j))';  % sample points
+                x_redraw = randsample(size(ParametersValues,1),Clustering.weights(j),0)';  % sample points
                 q = quantile(ParametersValues(x_redraw,i),(1:1:99)./100);  % bootstrapped distribution
                 boot(iter) = norm(q_prior-q,1); % L1-norm
             end
-            BootMainFactors(i,j) = quantile(boot,alpha);
+            if length(alpha)==1
+                BootMainFactors(i,j) = quantile(boot,alpha);
+            else
+                BootMainFactors(i,j,:) = quantile(boot,alpha);
+            end
         end
     end
 

DGSA_computations/ComputeConditionalEffects.m

@@ -0,0 +1,100 @@
+
+%
+% Distance-based generalized sensitivity analysis (dGSA)
+% Evaluation of sensitivity of the conditioanl effects
+
+% Author: Celine Scheidt
+% Date: August 2013
+
+% Updated by Jihoon Park (jhpark3@stanford.edu)
+% Date: 14 May 2016
+%   - Use data structures
+%   - Separate computations and visualizations
+%   - Display the conditioning parameter computed
+
+function DGSA_ConditionalEffects = ...
+   ComputeConditionalEffects(DGSA,alpha)%(Clustering,ParametersValues,ParametersNames,NbBins,alpha)
+
+
+%% Input Parameters
+
+% DGSA: (Stucture) Previous results of DGSA.
+%
+%   .Clustering: (Structure) Results of clustering
+%   .ParametersValues: (matrix, # of models x # of parameters) contains values of each parameter from Monte Carlo sampling. 
+%   .ConditionalEffects.NbBins: (vector, length=# of parameters) specifies each number of bins of parameters to compute conditional effects
+%   .ParametersNames: (cell, 1 x # of parameters) contains the names of parameters
+%
+% alpha: (optional) alpha-percentile (by default, alpha = 0.95) for bootstrap percentile
+
+%% Output Parameters 
+%
+% DGSA_ConditionalEffects: (Structure) Results of DGSA after computing conditional effects.
+%   .ConditionalEffects.SensitivityByClusterandBins:
+%        (4D array, (# of parameters) x (# of parameters-1) x (# of clusters) x (maximum value of bins))
+%        Element (i,j,c,l)  corresponds to conditional effect of  i-th parameter conditioned to  l-th bin of  j-th parameter from class c.  
+%
+%   .ConditionalEffects.StandardizedSensitivity: (vector, length=(# of parameters) x (# of parameters -1))
+%       Conditional effects averaged over bins and classes
+%
+%   .ConditionalEffects.H0accConditional (logical vector, length=(# of parameters) x (#  of parameters -1)) 
+%       Results of bootstrap hypothesis testing. If the value is 1, corresponding conditional effect is sensitive.
+%
+%   .ConditionalEffects.Names_ConditionalEffects: (cell, length=(# of parameters) x (#  of parameters -1)) 
+%       contains names of conditional effects
+
+%% Assign varaibles.
+DGSA_ConditionalEffects=DGSA;
+
+Clustering=DGSA.Clustering;
+ParametersValues=DGSA.ParametersValues;
+NbBins=DGSA.ConditionalEffects.NbBins;
+ParametersNames=DGSA.ParametersNames;
+
+
+%% Computations of conditional effects.
+
+NbParams = size(ParametersValues,2);
+NbClusters = size(Clustering.medoids,2);
+
+if nargin < 2; alpha = .95; end
+
+% Evaluate the normalized conditionnal interaction for each parameter, each
+% class and each bin
+L1Interactions = NaN(NbParams,NbParams-1,NbClusters,max(NbBins));  % array containing all the Interactions 
+BootInteractions = NaN(NbParams,NbParams-1,NbClusters,max(NbBins));  
+
+for params = 1:NbParams
+    fprintf('Computing sensitivity conditional to %s...\n',ParametersNames{params});
+    L1InteractionsParams = L1normInteractions(ParametersValues,params,Clustering,NbBins(params));
+    L1Interactions(params,:,:,1:NbBins(params)) = L1InteractionsParams(:,:,1:NbBins(params));
+    BootInteractionsParams = BootstrapInteractions(ParametersValues,params,Clustering,NbBins(params),2000,alpha);
+    BootInteractions(params,:,:,1:NbBins(params)) = BootInteractionsParams(:,:,1:NbBins(params));
+    NormalizedInteractions = L1Interactions./BootInteractions;
+end
+
+
+% Measure of conditional interaction sensitivity per class
+SensitivityPerClass = nanmean(NormalizedInteractions,4);
+
+% Average measure of sensitivity over all classes 
+SensitivityOverClass = nanmean(SensitivityPerClass,3);
+
+
+% Hypothesis test: Ho = 1 if at least one value > 1 (per cluster and per bin)
+SensitivityPerInteraction = reshape(permute(NormalizedInteractions,[2,1,4,3]),[],max(NbBins)*NbClusters);
+H0accInter = any(SensitivityPerInteraction >=1,2);
+
+StandardizedSensitivityInteractions = reshape(SensitivityOverClass',1,[]);
+%% Save variables to output data structure
+
+DGSA_ConditionalEffects.ConditionalEffects.SensitivityByClusterandBins=NormalizedInteractions;
+DGSA_ConditionalEffects.ConditionalEffects.StandardizedSensitivity=StandardizedSensitivityInteractions;
+DGSA_ConditionalEffects.ConditionalEffects.H0accConditional=H0accInter;
+DGSA_ConditionalEffects.ConditionalEffects.Names_ConditionalEffects=CreateNames_ConditionalEffect(ParametersNames);
+
+
+end
+
+
+

DGSA_computations/ComputeMainEffects.m

@@ -0,0 +1,140 @@
+function DataStructMainEffects=ComputeMainEffects(DataStruct,alpha,ColorType)
+% This function computes/display main effects from DGSA.
+% Author: Jihoon Park (jhpark3@stanford.edu) and Celine
+% Scheidt(scheidtc@stanford.edu)
+
+% Input:
+
+% DataSturc : Data structure having inputs required for computing main effects.
+%   .ParametersValues: (matrix, # of models x # of parameters) contains values of each parameter from Monte Carlo sampling.
+%   .ParametersNames: (cell, 1 x # of parameters) contains the names of parameters 
+%   .D: (matrix, # of models x # of models or vector, length =# of models x (# of models -1))
+%       Distance matrix or vector of distances between the model.
+%   .Clustering: (structure) results of clustering
+%   
+% 
+% Inputs for Display::
+%
+%   .MainEffects.Display.StandardizedSensitivity: (string) specifies the method to visualize standardized main effect. 
+%       'Pareto': Pareto plot
+%       'CI': Add a confidence interval to the Pareto plot
+%       'None': (default) do not display
+%   .MainEffects.Display.ParetoPlotbyCluster: 
+%
+% Optional:: 
+%   
+%   alpha: (optional) % quantile of bootstrapped idstribution. If only Pareto plot is used to display standardized main effects, It should be scalar. The default value is 0.95 
+%          If a confidence interval is displayed, it should be a vector of length 3. The value of sensitivity (length of Pareto bar) is based on the first element. 
+%   ColorType: (optional, string) specifies how the Pareto plot is colored when confidence interval is dispalyed. If 'RB', the bar is colord with Red/Blue scheme (default). 
+%          If 'jet', the bar is colored with 'jet' obejct from Matlab.
+
+
+
+% Output:
+
+% DataStructMainEffects.MainEffects: Results of Main effects are added to DataStruc
+%   .SensitivitybyCluster : (matrix, # of parameters x # of clusters) contains main effects by parameters (rows) and clusters (columns)
+%   .SensitivityStandardized : (vector, length=# of parameters) Main effects averaged over clusters.
+%   .H0accMain: (logical, length=# of parameters) Result of bootstrap hypothesis test. If true, a parameter is sensitive.
+
+
+%% Create variables
+DataStructMainEffects=DataStruct;
+if nargin<2
+    alpha=.95;
+end
+
+
+%% 1. Check the visualization method.
+if ~DataStruct.MainEffects.Display.StandardizedSensitivity
+    DataStruct.MainEffects.Display.StandardizedSensitivity='None'; % Default is not to visualize sensitivities
+end
+
+switch DataStruct.MainEffects.Display.StandardizedSensitivity
+
+    case 'CI'
+        if length(alpha)~=3 
+            error('In order to disply confidence interval, alpha should be a vector of length 3')
+        end
+
+        [~,idx]=sort(alpha);
+
+        if ~all(idx==[2,1,3])
+            error('In order to disply confidence interval, alpha(3)>alpha(1)>alpha(2)')
+        end
+
+
+
+    case 'Pareto'
+        if ~isscalar(alpha)
+            error('Significane level should be scalar for Parteo plot only')
+        end
+    case 'None'
+        if ~isscalar(alpha)
+            error('Significane level should be scalar when sensitivities not displayed')
+        end
+    otherwise
+        error('Enter the right type of display method');
+
+end
+
+
+%% 2. Create variables 
+
+ParametersValues=DataStruct.ParametersValues; % Values of Parameters
+ParametersNames=DataStruct.ParametersNames;
+Clustering=DataStruct.Clustering;
+
+
+
+
+%% 3. Compute Main Effects
+
+L1MainFactors = L1normMainFactors(Clustering,ParametersValues); % evaluate L1-norm
+L1MainFactors=repmat(L1MainFactors,1,1,length(alpha));
+
+BootMainFactors = BootstrapMainFactors(ParametersValues,Clustering,2000,alpha);  % bootstrap procedure
+
+SensitivityMainFactors = L1MainFactors./BootMainFactors; % normalization 
+SensitivitybyCluster=SensitivityMainFactors(:,:,1);
+StandardizedSensitivity = mean(SensitivityMainFactors,2); 
+
+
+%% 4. Hypothetis test: H0=1(reject, so sensitive) if at least one value > 1 (per cluster)
+H0accMain = any(SensitivitybyCluster >=1,2); 
+
+
+
+%% 5. Save Main effects to Data Structure
+SensitivityStandaridzed=StandardizedSensitivity(:,:,1);
+SensitivitybyCluster=SensitivityMainFactors(:,:,1);
+DataStructMainEffects.MainEffects.SensitivitybyCluster=SensitivitybyCluster;
+DataStructMainEffects.MainEffects.SensitivityStandardized=SensitivityStandaridzed;
+
+DataStructMainEffects.MainEffects.H0accMain=H0accMain;
+
+%% 6. Display main effects.
+
+% 6.1 Sensitivity by cluster.
+if DataStruct.MainEffects.Display.ParetoPlotbyCluster % display main effects by cluster  
+        ParetoMainFactors(SensitivitybyCluster,ParametersNames);    
+end
+
+% 6.2 Display Standardized Sensitivity
+
+switch DataStruct.MainEffects.Display.StandardizedSensitivity
+    case 'CI'
+        if nargin<3
+        ColorType='RB';
+        end
+        Pareto_GlobalSensitivity_CI(StandardizedSensitivity, H0accMain, ParametersNames,ColorType);
+
+    case 'Pareto'
+
+        Pareto_GlobalSensitivity(SensitivityStandaridzed,ParametersNames,H0accMain);
+
+
+end
+
+
+end

DGSA_computations/ComputeNetConditionalEffects.m

@@ -0,0 +1,75 @@
+function [NetConditionalEffects,Bin_MaxNetConditional]=ComputeNetConditionalEffects...
+    (ParameterSpecified, DGSA,YTicktxt)
+% This functions compute net interactions and visualizae in bar graphs
+% Author: Jihoon Park (jhpark3@stanford.edu)
+% Date: 16 May 2016
+%
+%Input: 
+% ParameterSpecified: string, a name of parameter that wants to see net conditioanl effects
+% DGSA: Datastructure contains results of DGSA.
+% YTicktxt: Labels of y-axis of horizontal bar graph.
+
+%Output:
+
+% NetConditionalEffects:matrix (NbBins x (# of Parameters-1)): contains net conditional effects of the specified parameter.
+% Bin_MaxNetConditional: scalar, the index of bin that has the highest conditional effects.
+
+
+
+%% Assgin varaiables
+
+ParametersNames=DGSA.ParametersNames;
+IdxParametersSpecified=find(ismember(ParametersNames,ParameterSpecified));
+NamesConditionalEffects=DGSA.ConditionalEffects.Names_ConditionalEffects;
+ParametersNames=DGSA.ParametersNames;
+ConditionalEffects_by_bins_clusters=DGSA.ConditionalEffects.SensitivityByClusterandBins;
+if isempty(IdxParametersSpecified)
+    error('Enter Right name of the parameter')
+end
+
+
+%% 
+NbParams=length(ParametersNames);
+
+SensitivityInteractionOverClass=nanmean(ConditionalEffects_by_bins_clusters, 3); % Average over cluster
+NetConditionalEffects_temp=SensitivityInteractionOverClass(IdxParametersSpecified,:,:,:); % Take according to varable specified.
+
+[~,NumConditionalPerParams,~,NbBins_Parameters]=size(NetConditionalEffects_temp);
+
+% Transform the form of variable NetConditionalEffects_temp
+
+NetConditionalEffects=zeros(NbBins_Parameters,NumConditionalPerParams); % Initialize
+
+for k=1:NbBins_Parameters
+    NetConditionalEffects(k,:)=NetConditionalEffects_temp(:,:,:,k);
+end
+
+%% Visualization
+
+
+figure;
+barh(NetConditionalEffects,'stacked');
+if nargin<3
+    if NbBins_Parameters==3
+        set(gca,'YTickLabel', {'Low','Med','High'});
+    end
+    if NbBins_Parameters==2
+        set(gca,'YTickLabel', {'Low','High'});
+    end
+else
+    set(gca,'YTickLabel', YTicktxt);
+
+end
+
+
+
+
+xlabel(ParametersNames(IdxParametersSpecified));
+ylabel('Net conditional effects');
+
+legend(NamesConditionalEffects((NbParams-1)*(IdxParametersSpecified-1)+1:(NbParams-1)*IdxParametersSpecified),'Location','EastOutside');        
+
+[~,Bin_MaxNetConditional]=max(sum(NetConditionalEffects,2));
+
+
+end