FatimaPillosu
commented
3 years ago Hi @onyb ,
on the basis of what Tim mentioned, I have created the Matlab code to show you what we might need to have in Python. The Matlab considers all the cases in the node that you are analysing. Then, the code runs the K-S for a pre-determined number of breakpoints determining the K-S's D statistic and the pvalue of the test (i.e. significance levels of the K-S test). Notice that the smaller the pvalue, the more significant the test result is. Finally, the code plots the Dstat and the (-1*log(pvalue)), in the same diagram for each of the considered breakpoints (see figure below). We have implemented a diagram with two different scales. The code provides you with the layout for the diagram that we might want to have, but if you have better ideas, please feel free to comment on them.
I have included in the code two different cases (i.e. case1: node at the top of the decision tree; case2: a node at the bottom of the decision tree). I did that to highlight the differences that you will encounter in these two situations. For case1, the pvalues are so so so so small that, even applying the log transformation, you cannot see anything, they are just set to 0. However, when you are at the bottom of the tree (where the node has a smaller number of cases), the pvalues start increasing and then you can see something in the diagram. Matlab applying a super high precision, so you can see that the pvalues, although very very very very small, they are not 0. I don't know how that might behold by Python. We need to see.
I'm putting this in priority 1 as, from what @ATimHewson told me, @EstiGascon and @AugustinVintzileos need this implementation as soon as possible. Thus, if there is something that it is not clear from Tim's explanation or my code, please, let's have a meeting as soon as you can. I think it would be quicker to discuss any problems directly rather than via email or via issues.
I add directly in the issue the Matlab code and the diagram for case2. I have also added in google drive the PDT.csv file that I used to create the plots, so we both start from the same data, and if results are different we can discuss them.
Cheers,
Fatima
Matlab Code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% TEST ON NEW WAY ON CREATING THE BREAKPOINTS FOR DECISION TREE %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The used K-S function is described here:
% https://uk.mathworks.com/help/stats/kstest2.html#btno0gd-ks2stat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUTS
InputData = 'PDT.mat';
NumBP = 20;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(InputData)
% Selection of the decision tree node to analyse and the predictor to consider
% % CASE 1
% predictor = "cpr";
% node = "Top of the tree";
% units = "-";
% VAR = cpr;
% CASE 2
predictor = "sr24h";
node = "(cpr>0.25 & cpr<0.75) & (tp>10 & tp<=inf) & (wspd700>20 & wspd700<=inf) & (cape>-inf & cape<=50)";
units = "J/kg";
VAR = sr24h(cpr>0.25 & cpr<0.75 & tp_acc>10 & tp_acc<=9999 & wspd700>20 & wspd700<=9999 & cape_wa>-9999 & cape_wa<=50);
[VAR, ind] = sort(VAR);
SizeVAR = length(VAR);
Len_SplitsVAR = round(SizeVAR/(NumBP+1));
BPs_ind = (1:NumBP)' * Len_SplitsVAR;
BPs = VAR(BPs_ind);
BPs = unique (BPs);
indBP_1 = min(VAR);
indBP_3 = max(VAR);
pvalue_vec = zeros(length(BPs),1);
Dstat_vec = zeros(length(BPs),1);
for i = 1 : length(BPs)
disp(i)
indBP_2 = BPs(i);
data2test_1 = FER(VAR>=indBP_1 & VAR<indBP_2);
data2test_2 = FER(VAR>=indBP_2 & VAR<=indBP_3);
if ~isempty(data2test_1) && ~isempty(data2test_2)
[h,pvalue,Dstat] = kstest2(data2test_1,data2test_2, 'Alpha', 0.05);
pvalue_vec(i) = pvalue;
Dstat_vec(i) = Dstat;
end
end
figure
yyaxis left
plot(BPs,Dstat_vec, "o-", "LineWidth", 1.5)
title({["Node"], node})
xlabel(strcat("Breakpoints for ", predictor, " [", units, "]"))
ylabel("K-S Dstat")
yyaxis right
plot(BPs, (-1 * log(pvalue_vec)), "o-", "LineWidth", 1.5)
ylabel("log(pvalue)")
AugustinVintzileos
onyb
We need a more complete way of using K-S tests to advise the user on potential breakpoints.
The key outputs we would like to send to the screen, to advise the user, are the standard K-S test outputs:
A. We would like to see 1 and 2, from the outset, across a range of possible breakpoints for all the data (this is different to what is currently done, where the code assumes successive removal of left-of-breakpoint data prior to each re-application of the K-S test). So the first set of K-S test runs like this might result in the user selecting one breakpoint.
B. After A the user will likely want to repeat the above for all data on one side of the selected breakpoint, to see if another breakpoint could be utilised there. And then they might look at the portion of data that lies on the other side of the first breakpoint in the same way. And then they might want to further divide up in the same way.
To start A we need a way of deciding how many independent K-S tests to run, and a way to output the results. It is difficult to advise on how many tests would be optimal - it is a compromise between run time, and dataset size. Perhaps 20 would be good value to start out with (this could perhaps be a user-defined value?).
Then how do we decide which 20 potential breakpoints to test? There are two main options. One would be to evenly divide up the governing variable range, between max and min, into 21 subsets and use those 20 breakpoints. A second would be to rank all governing variable values and divide into 21 equally sized subsets, and accordingly use the dividing points that correspond; these would of course be less evenly spread than in the first option. Probably the second case is preferable to give a general overview.
As regards conveying the output of the 20(?) tests, we could have i) tabular text, and/or ii) a graph. The graph is visually more appealing but would be more work. In each case we need to convey what the breakpoint is (in governing variable units), and what the values of 1 and 2 above are. If we went for the graphical option we would need two y-axes, a linear one for the D-statistics and non-linear for the P-value (some form of logarithmic?). Values for the P-value may be limited to a maximum of 99.99 (2 decimal places) if what I remember from my limited python experience is correct. This may prove a bit of a limitation, because as I understand it very high P values which could be useful to us strictly need a supercomputer to derive (!), but let's see.
With the above implemented the process of semi-subjective but informed selection of breakpoints by the user should be much improved (with clear graphical output as supporting evidence, if we can achieve that).
Priority for this issue should be to get something working with tabular output; the graphical output may take quite a bit longer and so is lower priority unless it can be done quickly.