
selecCriterion1=1; % 0=99.5%CI, 1=99%CI, 2=two most significant, -1=95%CI
selecCriterion2=1;
selecCriterion3=1;
selecCriterion4=1;
alpha=0.01;
intercept=1;
[n,d]=size(rawX);
lX=log(rawX);  
mLX=mean(lX,1);
meanLX=repmat(mLX,n,1);
dmLX=lX-meanLX;
X=dmLX;

dimHC5=ceil(nthroot(d,5));
nearestHypercube5=dimHC5^5;
remainderHC5=nearestHypercube5-d;
hypercube5=reshape([(1:d)';zeros(remainderHC5,1)],dimHC5,dimHC5,dimHC5,dimHC5,dimHC5);
hypercube5(:,:,:,:,find(all(hypercube5(:,1,1,1,:)==0,1)))=[]; 
hypercubeSelect=zeros(size(hypercube5));
for ind5=1:size(hypercube5,5)
    for ind4=1:size(hypercube5,4)
        for ind3=1:size(hypercube5,3)
            for indR=1:size(hypercube5,1)
                if size(nonzeros(hypercube5(indR,:,ind3,ind4,ind5)),1)>0
                    subsetX=X(:,nonzeros(hypercube5(indR,:,ind3,ind4,ind5)));
                    [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion1,intercept);
                    if count==1
                        hypercube5(indR,:,ind3,ind4,ind5)
                    end
                    if size(idxSelected,1)>0
                        hypercubeSelect(indR,idxSelected,ind3,ind4,ind5)=hypercubeSelect(indR,idxSelected,ind3,ind4,ind5)+1;
                    end %if
                end %if
            end %indR
            for indC=1:size(hypercube5,2)
                if size(nonzeros(hypercube5(:,indC,ind3,ind4,ind5)),1)>0
                    subsetX=X(:,nonzeros(hypercube5(:,indC,ind3,ind4,ind5)));
                    [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion1,intercept);
                    if count==1
                        hypercube5(:,indC,ind3,ind4,ind5)
                    end
                    if size(idxSelected,1)>0
                        hypercubeSelect(idxSelected,indC,ind3,ind4,ind5)=hypercubeSelect(idxSelected,indC,ind3,ind4,ind5)+1;
                    end %if
                end %if
            end %indC
        end %ind3
        for indR=1:size(hypercube5,1)
            for indC=1:size(hypercube5,2)
                if size(nonzeros(hypercube5(indR,indC,:,ind4,ind5)),1)>0
                    subsetX=X(:,nonzeros(hypercube5(indR,indC,:,ind4,ind5)));
                    [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion1,intercept);
                    if count==1
                        hypercube5(indR,indC,:,ind4,ind5)
                    end
                    if size(idxSelected,1)>0
                        hypercubeSelect(indR,indC,idxSelected,ind4,ind5)=hypercubeSelect(indR,indC,idxSelected,ind4,ind5)+1;
                    end %if
                end %if
            end %indC
        end %indR
    end %ind4
    for ind3=1:size(hypercube5,3)
         for indR=1:size(hypercube5,1)
             for indC=1:size(hypercube5,2)
                if size(nonzeros(hypercube5(indR,indC,ind3,:,ind5)),1)>0
                    subsetX=X(:,nonzeros(hypercube5(indR,indC,ind3,:,ind5)));
                    [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion1,intercept);
                    if count==1
                        hypercube5(indR,indC,ind3,:,ind5)
                    end
                    if size(idxSelected,1)>0
                        hypercubeSelect(indR,indC,ind3,idxSelected,ind5)=hypercubeSelect(indR,indC,ind3,idxSelected,ind5)+1;
                    end %if
                end %if
            end %indR
         end
    end
end %ind5
% now traverse in the 5th dimension
for ind4=1:size(hypercube5,4)
    for ind3=1:size(hypercube5,3)
        for indR=1:size(hypercube5,1)
            for indC=1:size(hypercube5,2)
                if size(nonzeros(hypercube5(indR,indC,ind3,ind4,:)),1)>0
                    subsetX=X(:,nonzeros(hypercube5(indR,indC,ind3,ind4,:)));
                    [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion1,intercept);
                    if count==1
                        hypercube5(indR,indC,ind3,ind4,:)
                    end
                    if size(idxSelected,1)>0
                        hypercubeSelect(indR,indC,ind3,ind4,idxSelected)=hypercubeSelect(indR,indC,ind3,ind4,idxSelected)+1;
                    end %if
                end
            end %indC
        end %indR
    end %ind3
end
setSelected4Times5=hypercube5(find(hypercubeSelect>3)); 
setSelected3Times5=hypercube5(find(hypercubeSelect>2)); 
setSelectedTwice5=hypercube5(find(hypercubeSelect>1));
numSelected4times5=length(setSelected4Times5)
numSelected3times5=length(setSelected3Times5)
numSelectedTwice5=length(setSelectedTwice5)
assert(sum(ismember(setSelected3Times5,setSelectedTwice5))==numSelected3times5)


dimHC4=ceil(nthroot(numSelected3times5,4));
nearestHypercube4=dimHC4^4;
remainderHC4=nearestHypercube4-numSelected3times5;
hypercube4=reshape([setSelected3Times5;zeros(remainderHC4,1)],dimHC4,dimHC4,dimHC4,dimHC4);
hypercube4(:,:,:,find(all(hypercube4(:,1,1,:)==0,1)))=[]; 
hypercubeSelect4=zeros(size(hypercube4));
for ind4=1:size(hypercube4,4)
    for ind3=1:size(hypercube4,3)
        for indR=1:size(hypercube4,1)
            if size(nonzeros(hypercube4(indR,:,ind3,ind4)),1)>0
                subsetX=X(:,nonzeros(hypercube4(indR,:,ind3,ind4)));
                [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion2,intercept);
                if count==1
                    hypercube4(indR,:,ind3,ind4)
                end
                if size(idxSelected,1)>0
                    hypercubeSelect4(indR,idxSelected,ind3,ind4)=hypercubeSelect4(indR,idxSelected,ind3,ind4)+1;
                end %if
            end %if
        end %indR
        for indC=1:size(hypercube4,2)
            if size(nonzeros(hypercube4(:,indC,ind3,ind4)),1)>0
                subsetX=X(:,nonzeros(hypercube4(:,indC,ind3,ind4)));
                [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion2,intercept);
                if count==1
                    hypercube4(:,indC,ind3,ind4)
                end
                if size(idxSelected,1)>0
                    hypercubeSelect4(idxSelected,indC,ind3,ind4)=hypercubeSelect4(idxSelected,indC,ind3,ind4)+1;
                end %if
            end %if
        end %indC
    end %ind3
    for indR=1:size(hypercube4,1)
        for indC=1:size(hypercube4,2)
            if size(nonzeros(hypercube4(indR,indC,:,ind4)),1)>0
                subsetX=X(:,nonzeros(hypercube4(indR,indC,:,ind4)));
                [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion2,intercept);   
                if count==1
                    hypercube4(indR,indC,:,ind4)
                end
                if size(idxSelected,1)>0
                    hypercubeSelect4(indR,indC,idxSelected,ind4)=hypercubeSelect4(indR,indC,idxSelected,ind4)+1;
                end %if
            end %if
        end %indC
    end %indR
end %ind4
% now traverse in the 4th dimension
for ind3=1:size(hypercube4,3)
    for indR=1:size(hypercube4,1)
        for indC=1:size(hypercube4,2)
            if size(nonzeros(hypercube4(indR,indC,ind3,:)),1)>0
                subsetX=X(:,nonzeros(hypercube4(indR,indC,ind3,:)));
                [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion2,intercept);
                if count==1
                    hypercube4(indR,indC,ind3,:)
                end
                if size(idxSelected,1)>0
                    hypercubeSelect4(indR,indC,ind3,idxSelected)=hypercubeSelect4(indR,indC,ind3,idxSelected)+1;
                end %if
            end
        end %indC
    end %indR
end %ind3
setSelected4Times=hypercube4(find(hypercubeSelect4>3)); 
setSelected3Times=hypercube4(find(hypercubeSelect4>2)); 
setSelectedTwice=hypercube4(find(hypercubeSelect4>1));
numSelected4times=length(setSelected4Times)
numSelected3times=length(setSelected3Times)
numSelectedTwice=length(setSelectedTwice)
assert(sum(ismember(setSelected3Times,setSelectedTwice))==numSelected3times)



dimCube=ceil(nthroot(numSelectedTwice,3)); 
nearestCube=dimCube^3;
remainderCube=nearestCube-numSelectedTwice;
cube=reshape([setSelectedTwice;zeros(remainderCube,1)],dimCube,dimCube,dimCube);
cube(:,:,all(cube(:,1,:)==0,1))=[];
cubeSelect=zeros(size(cube));
for indL=1:size(cube,3) 
    for indR=1:size(cube,1) 
        if size(nonzeros(cube(indR,:,indL)),1)>0
            subsetX=X(:,nonzeros(cube(indR,:,indL)));
            [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion3,intercept);
            if count==1
                squeeze(cube(indR,:,indL))
            end 
            if size(idxSelected,1)>0
                cubeSelect(indR,idxSelected,indL)=cubeSelect(indR,idxSelected,indL)+1;
            end
        end
    end %indR
    for indC=1:size(cube,2) 
        if size(nonzeros(cube(:,indC,indL)),1)>0
            subsetX=X(:,nonzeros(cube(:,indC,indL)));
            [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion3,intercept);
            if count==1
                squeeze(cube(:,indC,indL))
            end
            if size(idxSelected,1)>0
                cubeSelect(idxSelected,indC,indL)=cubeSelect(idxSelected,indC,indL)+1;
            end
        end
    end % indC
end %indL
for indR=1:size(cube,1)
    for indC=1:size(cube,2)
        if size(nonzeros(cube(indR,indC,:)),1)>0
            subsetX=X(:,nonzeros(cube(indR,indC,:)));
            [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion3,intercept);
            if count==1
                squeeze(cube(indR,indC,:))
            end
            if size(idxSelected,1)>0
               cubeSelect(indR,indC,idxSelected)=cubeSelect(indR,indC,idxSelected)+1;
            end
        end
    end % indC
end %indR
setSelected322=cube(find(cubeSelect>1)); 
numSelected322=length(setSelected322)

if numSelected322<25
    setSelectedB=setSelected322;
else
    dimSquare2=ceil(sqrt(numSelected322)); 
    nearestSquare2=dimSquare2^2;
    remainderSquare2=nearestSquare2-numSelected322;
    square2=reshape([setSelected322;zeros(remainderSquare2,1)],dimSquare2,dimSquare2);
    square2(:,all(square2==0,1))=[];
    squareSelect2=zeros(size(square2));
    for indR=1:size(square2,1)
        if size(nonzeros(square2(indR,:)),1)>0
            subsetX=X(:,nonzeros(square2(indR,:)));
            idxSelected=glmRoutine(subsetX,Y,n,selecCriterion4,intercept);
            if size(idxSelected,1)>0
                squareSelect2(indR,idxSelected)=squareSelect2(indR,idxSelected)+1;
            end
        end
    end
    for indC=1:size(square2,2)
        if size(nonzeros(square2(:,indC)),1)>0
            subsetX=X(:,nonzeros(square2(:,indC)));
            idxSelected=glmRoutine(subsetX,Y,n,selecCriterion4,intercept);
            if size(idxSelected,1)>0
                squareSelect2(idxSelected,indC)=squareSelect2(idxSelected,indC)+1;
            end
        end
    end
    if any(squareSelect2(:)>1)
        setSelected3222=square2(find(squareSelect2>1)); 
        numSelected3222=length(setSelected3222)
    else
        setSelected3222=[];
        numSelected3222=0
    end
    if any(squareSelect2(:)>0)
        setSelected3221=square2(find(squareSelect2>0)); 
        numSelected3221=length(setSelected3221)
    else
        setSelected3221=[];
        numSelected3221=0
    end
    
    if (numSelected3222>10 && numSelected3221>20)
        setSelectedB=setSelected3222;
    else
        setSelectedB=setSelected3221;
    end
end
numSelectedB=length(setSelectedB)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% Exploratory phase - judgement required %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


corrSelected=corrcoef(X(:,setSelectedB));
count=0;
for j=1:numSelectedB-1
    for k=j+1:numSelectedB
        if corrSelected(j,k)>0.9
            sprintf('%d,%d',j,k)
            count=count+1
        end
    end
end

tStats=zeros(numSelectedB,2);
pVals=zeros(numSelectedB,2);
for j=1:numSelectedB
    xTmp=zeros(n,2);
    xTmp(:,1)=X(:,setSelectedB(j,1));
    xTmp(:,2)=(X(:,setSelectedB(j,1))).^2;
    [~,~,stats]=glmfit(xTmp,[Y,ones(n,1)],'binomial');
    tSt=stats.t;
    tSt(1,:)=[]; 
    tStats(j,:)=tSt;
    pVs=stats.p;
    pVs(1,:)=[];
    pVals(j,:)=pVs;
end

unifQuantiles=1:1:numSelectedB;
unifQuantiles=unifQuantiles'-0.5;
unifQuantiles=unifQuantiles./numSelectedB;
normQuantiles=norminv(unifQuantiles,0,1);
z1=[-4:0.01:4]';
figure()
hold on
scatter(sort(tStats(:,1)),normQuantiles, 60, 'MarkerEdgeColor', [0 0.4470 0.7410],'LineWidth',2);
xlabel('Empirical quantiles of t statistics', 'FontSize', 12)
ylabel('Theoretical quantiles of the standard normal distribution', 'FontSize', 12)
ylim([-4,4])
xlim([-4,4])
scatter(sort(tStats(:,2)),normQuantiles, 60, '^', 'MarkerEdgeColor', [0.47,0.67,0.19],'LineWidth',2);
plot(z1,z1,'-','Color',[0.85 0.33 0])
hold off
box on


z3=[0:0.01:1]';
confU5=zeros(numSelectedB,1);
confL5=zeros(numSelectedB,1);
confU1=zeros(numSelectedB,1);
confL1=zeros(numSelectedB,1);
for i=1:numSelectedB
    confU5(i,1)=betainv(1-0.05/2,i,numSelectedB-i+1);
    confL5(i,1)=betainv(0.05/2,i,numSelectedB-i+1);
    confU1(i,1)=betainv(1-0.01/2,i,numSelectedB-i+1);
    confL1(i,1)=betainv(0.01/2,i,numSelectedB-i+1);
end
figure()
hold on
scatter(sort(pVals(:,1)),unifQuantiles, 60, 'MarkerEdgeColor', [0 0.4470 0.7410],'LineWidth',2);
xlabel('Empirical quantiles of p values', 'FontSize', 12)
ylabel('Theoretical quantiles of the uniform distribution', 'FontSize', 12)
ylim([0,1])
xlim([0,1])
scatter(sort(pVals(:,2)),unifQuantiles, 60, '^', 'MarkerEdgeColor', [0.47,0.67,0.19],'LineWidth',2);
plot(z3,z3,'-','Color',[0.85 0.33 0])
plot(unifQuantiles,confU5,'--','Color',[0.85 0.33 0])
plot(unifQuantiles,confL5,'--','Color',[0.85 0.33 0])
plot(unifQuantiles,confU1,'--','Color',[0.49 0.18 0.56])
plot(unifQuantiles,confL1,'--','Color',[0.49 0.18 0.56])
hold off
box on

lcU1=norminv(confU1,0,1);
lcL1=norminv(confL1,0,1);
lcU5=norminv(confU5,0,1);
lcL5=norminv(confL5,0,1);
plZ=norminv(unifQuantiles,0,1);
figure()
hold on
scatter(sort(tStats(:,1)),normQuantiles, 60, 'MarkerEdgeColor', [0 0.4470 0.7410],'LineWidth',2);
xlabel('Empirical quantiles of t statistics', 'FontSize', 12)
ylabel('Theoretical quantiles of the standard normal distribution', 'FontSize', 12)
ylim([-5,5])
xlim([-5,5])
scatter(sort(tStats(:,2)),normQuantiles, 60, '^', 'MarkerEdgeColor', [0.47,0.67,0.19],'LineWidth',2);
plot(z1,z1,'-','Color',[0.85 0.33 0])
plot(plZ,lcU5,'--','Color',[0.85 0.33 0])
plot(plZ,lcL5,'--','Color',[0.85 0.33 0])
plot(plZ,lcU1,'--','Color',[0.49 0.18 0.56])
plot(plZ,lcL1,'--','Color',[0.49 0.18 0.56])
hold off
box on


setOriginalScale=setSelectedB; % This doesn't do anything any more.
numSetOrigScale=length(setOriginalScale);

tStatsIntVec=zeros(numSetOrigScale*numSetOrigScale,3);
pValsIntVec=zeros(numSetOrigScale*numSetOrigScale,3);
for j=1:numSetOrigScale-1
    for k=j+1:numSetOrigScale
        xTmp=zeros(n,3);
        xTmp(:,1)=X(:,setOriginalScale(j,1));
        xTmp(:,2)=X(:,setOriginalScale(k,1));
        xTmp(:,3)=X(:,setOriginalScale(j,1)).*X(:,setOriginalScale(k,1));
        [~,~,stats]=glmfit(xTmp,[Y,ones(n,1)],'binomial');
        linearInd=sub2ind([numSetOrigScale,numSetOrigScale],j,k);
        tSt=stats.t;
        tSt(1,:)=[]; % don't care about intercept
        tStatsIntVec(linearInd,:)=tSt;
        pVs=stats.p;
        pVs(1,:)=[];
        pValsIntVec(linearInd,:)=pVs;
    end
end

tStatsIntVecNoZeros=tStatsIntVec;
pValsIntVecNoZeros=pValsIntVec;
tStatsIntVecNoZeros(all(tStatsIntVecNoZeros==0,2),:)=[];
pValsIntVecNoZeros(all(pValsIntVecNoZeros==0,2),:)=[];

numComparisons=size(tStatsIntVecNoZeros,1);
unifQuantiles2=1:1:numComparisons;
unifQuantiles2=unifQuantiles2'-0.5;
unifQuantiles2=unifQuantiles2./numComparisons;
normQuantiles2=norminv(unifQuantiles2,0,1);

z1=[-4:0.01:4]';
figure()
hold on
xlabel('Empirical quantiles of t statistics', 'FontSize', 12)
ylabel('Theoretical quantiles of the standard normal distribution', 'FontSize', 12)
scatter(sort(tStatsIntVecNoZeros(:,3)),normQuantiles2, 60, 'v', 'MarkerEdgeColor', [0 0.4470 0.7410],'LineWidth',2);
ylim([-4,4])
xlim([-4,4])
plot(z1,z1,'-','Color',[0.85 0.33 0])
hold off
box on



grey=[0.8,0.8,0.8];
confU2_5=zeros(numComparisons,1);
confL2_5=zeros(numComparisons,1);
confU2_1=zeros(numComparisons,1);
confL2_1=zeros(numComparisons,1);
for i=1:numComparisons
    confU2_5(i,1)=betainv(1-0.05/2,i,numComparisons-i+1);
    confL2_5(i,1)=betainv(0.05/2,i,numComparisons-i+1);
    confU2_1(i,1)=betainv(1-0.01/2,i,numComparisons-i+1);
    confL2_1(i,1)=betainv(0.01/2,i,numComparisons-i+1);
end
figure()
hold on
scatter(sort(pValsIntVecNoZeros(:,3)),unifQuantiles2, 60, 'v', 'MarkerEdgeColor', [0 0.4470 0.7410],'LineWidth',2);
xlabel('Empirical quantiles of p values', 'FontSize', 12)
ylabel('Theoretical quantiles of the uniform distribution', 'FontSize', 12)
ylim([0,1])
xlim([0,1])
plot(z3,z3,'-','Color',[0.85 0.33 0])
plot(unifQuantiles2,confU2_5,'--','Color',[0.85 0.33 0])
plot(unifQuantiles2,confL2_5,'--','Color',[0.85 0.33 0])
plot(unifQuantiles2,confU2_1,'--','Color',[0.49 0.18 0.56])
plot(unifQuantiles2,confL2_1,'--','Color',[0.49 0.18 0.56])
hold off
box on


lcU2_1=norminv(confU2_1,0,1);
lcL2_1=norminv(confL2_1,0,1);
lcU2_5=norminv(confU2_5,0,1);
lcL2_5=norminv(confL2_5,0,1);
plZ2=norminv(unifQuantiles2,0,1);
z1=[-4:0.01:4]';
figure()
hold on
xlabel('Empirical quantiles of t statistics', 'FontSize', 12)
ylabel('Theoretical quantiles of the standard normal distribution', 'FontSize', 12)
scatter(sort(tStatsIntVecNoZeros(:,3)),normQuantiles2, 60, 'v', 'MarkerEdgeColor', [0 0.4470 0.7410],'LineWidth',2);
ylim([-4,4])
xlim([-4,4])
plot(z1,z1,'-','Color',[0.85 0.33 0])
plot(plZ2,lcU2_5,'--','Color',[0.85 0.33 0])
plot(plZ2,lcL2_5,'--','Color',[0.85 0.33 0])
plot(plZ2,lcU2_1,'--','Color',[0.49 0.18 0.56])
plot(plZ2,lcL2_1,'--','Color',[0.49 0.18 0.56])
hold off
box on

sort(pValsIntVecNoZeros(:,3))
idxInclude3=find(tStatsIntVecNoZeros(:,3)<-2.5);
numInteractions=length(idxInclude3); 
[pairIndicesj,pairIndicesk]=ind2sub([numSetOrigScale,numSetOrigScale],idxInclude3);
interactingPairs=[setOriginalScale(pairIndicesj,1),setOriginalScale(pairIndicesk,1)]

% visually inspect all interactions that are suspect on the basis of their
% probability plots, e.g.

figure()
scatter(X(:,25470),X(:,11643),50,Y,'filled')
box on
figure()
scatter(X(:,27642),X(:,44276),50,Y,'filled')
box on
figure(1)
xlabel('variable 25470', 'FontSize', 14)
ylabel('variable 11643', 'FontSize', 14)
figure(2)
xlabel('variable 27642', 'FontSize', 14)
ylabel('variable 44276', 'FontSize', 14)


% discard the ones that do not look visually plausible, e.g.

interactingPairs(1,:)=[]
numInteractions=numInteractions-1


idxIncludeExploreMore=find(abs(tStatsIntVecNoZeros(:,3))>2);
numInteractions=length(idxIncludeExploreMore);
[pairIndicesjExploreMore,pairIndiceskExploreMore]=ind2sub([numSetOrigScale,numSetOrigScale],idxIncludeExploreMore);
interactingPairsExploreMore=[setOriginalScale(pairIndicesjExploreMore,1),setOriginalScale(pairIndiceskExploreMore,1)] % i


figure()
scatter(X(1:105,27642),X(1:105,7235),80,'yellow','filled')
box on
hold on
scatter(X(106:end,27642),X(106:end,7235),80,'blue','filled','^')
hold on
scatter(X([116,117,118,122],27642),X([116,117,118,122],7235),150,'red','filled','s')
xlabel('variable 27642', 'FontSize', 18)
ylabel('variable 7235', 'FontSize', 18)


figure()
scatter(X(:,27920),X(:,27642),50,Y,'filled')
box on
hold on
scatter(X([116,117,118,122],27920),X([116,117,118,122],27642),50,'r','filled')
xlabel('variable 27920', 'FontSize', 18)
ylabel('variable 27642', 'FontSize', 18)


figure()
scatter(X(:,27642),X(:,44276),50,Y,'filled')
box on
hold on
scatter(X([116,117,118,122],27642),X([116,117,118,122],44276),50,'r','filled')
xlabel('variable 27642', 'FontSize', 18)
ylabel('variable 44276', 'FontSize', 18)


figure()
scatter(X(:,25744),X(:,25125),50,Y,'filled')
box on
hold on
scatter(X([116,117,118,122],25744),X([116,117,118,122],25125),50,'r','filled')
xlabel('variable 25744', 'FontSize', 18)
ylabel('variable 25125', 'FontSize', 18)

% these individuals look suspicious: they break an otherwise strong
% multivariate relationship. Discard them.
outliers=[116,117,118,122]';


interactingPairsExploreMore2=interactingPairsExploreMore([2,8,10,11],:)
numInteractions=size(interactingPairsExploreMore2,1);

xInteractions=[];
for i=1:numInteractions
    xInteractions=[xInteractions, X(:,interactingPairsExploreMore2(i,1)).*X(:,interactingPairsExploreMore2(i,2))]; 
end

nOutliers=size(outliers,1);
XNoOutliers=X;
XNoOutliers(outliers,:)=[];
YNoOutliers=Y;
YNoOutliers(outliers,:)=[];
xInteractNoOutliers=xInteractions;
xInteractNoOutliers(outliers,:)=[];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Model selection phase %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


alpha=0.01;
xAugment2NoOutliers=[XNoOutliers(:,setOriginalScale),xInteractNoOutliers];
idxXAugment=1:size(xAugment2NoOutliers,2);
goodModels=cell(size(xAugment2NoOutliers,2),1);
[~,devComprehensive,stats]=glmfit(xAugment2NoOutliers,[YNoOutliers,ones(n-nOutliers,1)],'binomial');
%warning('off','last') % surpress the warning from each of the following model fits for speed
%goodFit=0;
%k=1;
%while goodFit==0
%    k
%    combinationMatrix=combntns(idxXAugment,k);
%    k
%    logicFitVector=zeros(size(combinationMatrix,1),1);
%    for j=1:size(combinationMatrix,1)
%        XSelect=xAugment2NoOutliers(:,combinationMatrix(j,:));
%        [~,dev,stats]=glmfit(XSelect,[YNoOutliers,ones(n-nOutliers,1)],'binomial');
%        nParams=size(XSelect,2);
%        fCrit=finv(1-alpha,size(xAugment2NoOutliers,2)-nParams,(n-nOutliers)-size(xAugment2NoOutliers,2));
%        if ((dev-devComprehensive)/(size(xAugment2NoOutliers,2)-nParams))/(devComprehensive/((n-nOutliers)-size(xAugment2NoOutliers,2)))<=fCrit
%            logicFitVector(j,1)=1;
%        end
%        idxFitVec=find(logicFitVector);
%        goodModels{k}=combinationMatrix(idxFitVec,:);
%    end
%    if sum(logicFitVector)>0
%       goodFit=1;
%    end
%    k=k+1;
%end     


% terminate the loop below once k is large
% enough.
goodModels=cell(size(xAugment2NoOutliers,2),1);
for k=3:10
    k
    combinationMatrix=combntns(idxXAugment,k);
    k
    logicFitVector=zeros(size(combinationMatrix,1),1);
    for j=1:size(combinationMatrix,1)
        XSelect=xAugment2NoOutliers(:,combinationMatrix(j,:));
        [~,dev,~]=glmfit(XSelect,[YNoOutliers,ones(n-nOutliers,1)],'binomial');
        nParams=size(XSelect,2);
        fCrit=finv(1-alpha,size(xAugment2NoOutliers,2)-nParams,(n-nOutliers)-size(xAugment2NoOutliers,2));
        if ((dev-devComprehensive)/(size(xAugment2NoOutliers,2)-nParams))/(devComprehensive/(n-size(xAugment2NoOutliers,2)))<=fCrit
            logicFitVector(j,1)=1;
        end
        if sum(logicFitVector)>0
           goodFit=1
        end
        idxFitVec=find(logicFitVector);
        goodModels{k}=combinationMatrix(idxFitVec,:);
    end
end     




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% processesing the output in a form suitable for presentation %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


models5=goodModels{5};
mod5VarIndices=zeros(size(models5));
for k=1:size(models5,1)
    for j=1:size(models5,2)
        if models5(k,j)<size(setOriginalScale,1)+1
            mod5VarIndices(k,j)=setOriginalScale(models5(k,j),1);
        elseif models5(k,j)>size(setOriginalScale,1)
            mod5VarIndices(k,j)=models5(k,j)-size(setOriginalScale,1);
        end
    end
end

models6=goodModels{6};
mod6VarIndices=zeros(size(models6));
for k=1:size(models6,1)
    for j=1:size(models6,2)
        if models6(k,j)<size(setOriginalScale,1)+1
            mod6VarIndices(k,j)=setOriginalScale(models6(k,j),1);
        elseif models6(k,j)>size(setOriginalScale,1)
            mod6VarIndices(k,j)=models6(k,j)-size(setOriginalScale,1);
        end
    end
end


models7=goodModels{7};
mod7VarIndices=zeros(size(models7));
for k=1:size(models7,1)
    for j=1:size(models7,2)
        if models7(k,j)<size(setOriginalScale,1)+1
            mod7VarIndices(k,j)=setOriginalScale(models7(k,j),1);
        elseif models7(k,j)>size(setOriginalScale,1)
            mod7VarIndices(k,j)=models7(k,j)-size(setOriginalScale,1);
        end
    end
end

listUnique5=unique(unique(mod7VarIndices))

tableModels5=zeros(size(mod5VarIndices,1),size(listUnique5,1));
for i=1:size(mod5VarIndices,1)
    for k=1:size(listUnique5,1)
       for j=1:size(mod5VarIndices,2)
           if mod7VarIndices(i,j)==7235
               tableModels5(i,1)=1;
           elseif mod5VarIndices(i,j)==11643
               tableModels5(i,2)=2;
           elseif mod5VarIndices(i,j)==25125
               tableModels5(i,3)=3;
           elseif mod5VarIndices(i,j)==25470
               tableModels5(i,4)=4;
           elseif mod5VarIndices(i,j)==25744
               tableModels5(i,5)=5;
           elseif mod5VarIndices(i,j)==27642
               tableModels5(i,6)=6;
           elseif mod5VarIndices(i,j)==27920
               tableModels5(i,7)=7;
           elseif mod5VarIndices(i,j)==29679
               tableModels5(i,8)=8;
           elseif mod5VarIndices(i,j)==33385
               tableModels5(i,9)=9;
           elseif mod5VarIndices(i,j)==36409
               tableModels5(i,10)=10;
           elseif mod5VarIndices(i,j)==37443
               tableModels5(i,11)=11;
           elseif mod5VarIndices(i,j)==44276
               tableModels5(i,12)=12;
           elseif mod5VarIndices(i,j)==45991
               tableModels5(i,13)=13;
           elseif mod5VarIndices(i,j)==46771
               tableModels5(i,14)=14;
           elseif mod5VarIndices(i,j)==48415
               tableModels5(i,15)=15;
           elseif mod5VarIndices(i,j)==48433
               tableModels5(i,16)=16;
           elseif mod5VarIndices(i,j)==48549
               tableModels5(i,17)=17;
           elseif mod5VarIndices(i,j)==1
               tableModels5(i,18)=18;
           elseif mod5VarIndices(i,j)==2
               tableModels5(i,19)=19;
           elseif mod5VarIndices(i,j)==3
               tableModels5(i,20)=20;
           elseif mod5VarIndices(i,j)==4
               tableModels5(i,21)=21;
           end
       end
    end
end
      

listUnique6=unique(unique(mod6VarIndices))

tableModels6=zeros(size(mod6VarIndices,1),size(listUnique6,1));
for i=1:size(mod6VarIndices,1)
    for k=1:size(listUnique6,1)
       for j=1:size(mod6VarIndices,2)
           if mod6VarIndices(i,j)==7235
               tableModels6(i,1)=1;
           elseif mod6VarIndices(i,j)==11643
               tableModels6(i,2)=2;
           elseif mod6VarIndices(i,j)==25125
               tableModels6(i,3)=3;
           elseif mod6VarIndices(i,j)==25470
               tableModels6(i,4)=4;
           elseif mod6VarIndices(i,j)==25744
               tableModels6(i,5)=5;
           elseif mod6VarIndices(i,j)==27642
               tableModels6(i,6)=6;
           elseif mod6VarIndices(i,j)==27920
               tableModels6(i,7)=7;
           elseif mod6VarIndices(i,j)==29679
               tableModels6(i,8)=8;
           elseif mod6VarIndices(i,j)==33385
               tableModels6(i,9)=9;
           elseif mod6VarIndices(i,j)==36409
               tableModels6(i,10)=10;
           elseif mod6VarIndices(i,j)==37443
               tableModels6(i,11)=11;
           elseif mod6VarIndices(i,j)==44276
               tableModels6(i,12)=12;
           elseif mod6VarIndices(i,j)==45991
               tableModels6(i,13)=13;
           elseif mod6VarIndices(i,j)==46771
               tableModels6(i,14)=14;
           elseif mod6VarIndices(i,j)==48415
               tableModels6(i,15)=15;
           elseif mod6VarIndices(i,j)==48433
               tableModels6(i,16)=16;
           elseif mod6VarIndices(i,j)==48549
               tableModels6(i,17)=17;
           elseif mod6VarIndices(i,j)==1
               tableModels6(i,18)=18;
           elseif mod6VarIndices(i,j)==2
               tableModels6(i,19)=19;
           elseif mod6VarIndices(i,j)==3
               tableModels6(i,20)=20;
           elseif mod6VarIndices(i,j)==4
               tableModels6(i,21)=21;
           end
       end
    end
end


listUnique7=unique(unique(mod7VarIndices))

tableModels7=zeros(size(mod7VarIndices,1),size(listUnique7,1));
for i=1:size(mod7VarIndices,1)
    for k=1:size(listUnique7,1)
       for j=1:size(mod7VarIndices,2)
           if mod7VarIndices(i,j)==7235
               tableModels7(i,1)=1;
           elseif mod7VarIndices(i,j)==11643
               tableModels7(i,2)=2;
           elseif mod7VarIndices(i,j)==25125
               tableModels7(i,3)=3;
           elseif mod7VarIndices(i,j)==25470
               tableModels7(i,4)=4;
           elseif mod7VarIndices(i,j)==25744
               tableModels7(i,5)=5;
           elseif mod7VarIndices(i,j)==27642
               tableModels7(i,6)=6;
           elseif mod7VarIndices(i,j)==27920
               tableModels7(i,7)=7;
           elseif mod7VarIndices(i,j)==29679
               tableModels7(i,8)=8;
           elseif mod7VarIndices(i,j)==33385
               tableModels7(i,9)=9;
           elseif mod7VarIndices(i,j)==36409
               tableModels7(i,10)=10;
           elseif mod7VarIndices(i,j)==37443
               tableModels7(i,11)=11;
           elseif mod7VarIndices(i,j)==44276
               tableModels7(i,12)=12;
           elseif mod7VarIndices(i,j)==45991
               tableModels7(i,13)=13;
           elseif mod7VarIndices(i,j)==46771
               tableModels7(i,14)=14;
           elseif mod7VarIndices(i,j)==48415
               tableModels7(i,15)=15;
           elseif mod7VarIndices(i,j)==48433
               tableModels7(i,16)=16;
           elseif mod7VarIndices(i,j)==48549
               tableModels7(i,17)=17;
           elseif mod7VarIndices(i,j)==1
               tableModels7(i,18)=18;
           elseif mod7VarIndices(i,j)==2
               tableModels7(i,19)=19;
           elseif mod7VarIndices(i,j)==3
               tableModels7(i,20)=20;
           elseif mod7VarIndices(i,j)==4
               tableModels7(i,21)=21;
           end
       end
    end
end

tableModelsAll=[tableModels5;tableModels6;tableModels7];
logicTableModelsAll=(tableModelsAll>0);
occuranceRate=sum(logicTableModelsAll,1)./size(logicTableModelsAll,1)
[sortedOccuranceRate,idxSorted]=sort(occuranceRate,'descend')
tableModelsReordered=tableModelsAll(:,idxSorted)
% should have added the main effects for all interaction terms
tableModelsInteractions=tableModelsAll;
for i=1:size(tableModelsInteractions,1)
    if tableModelsInteractions(i,18)==18
        tableModelsInteractions(i,6)=6;
        tableModelsInteractions(i,1)=1;
    end
    if tableModelsInteractions(i,19)==19
        tableModelsInteractions(i,5)=5;
        tableModelsInteractions(i,3)=3;
    end
    if tableModelsInteractions(i,20)==20
        tableModelsInteractions(i,7)=7;
        tableModelsInteractions(i,6)=6;
    end
    if tableModelsInteractions(i,21)==21
        tableModelsInteractions(i,12)=12;
        tableModelsInteractions(i,6)=6;
    end
end
logicTableModelsInteractions=(tableModelsInteractions>0);
occuranceRateInteract=sum(logicTableModelsInteractions,1)./size(logicTableModelsInteractions,1)
[sortedOccuranceRateInteract,idxSortedInteract]=sort(occuranceRateInteract,'descend')
tableModelsReorderedInteract=tableModelsInteractions(:,idxSortedInteract)
logicTableModelsReorderedInteract=(tableModelsReorderedInteract>0)
idxNot16=find(logicTableModelsReorderedInteract(:,2)==0)
logicTableNot16=logicTableModelsReorderedInteract(idxNot16,:)
sum(logicTableNot16(:,11))./size(logicTableNot16,1)
sum(logicTableNot16(:,10))./size(logicTableNot16,1)
idxNot1=find(logicTableModelsReorderedInteract(:,1)==0);
logicTableNot1=logicTableModelsReorderedInteract(idxNot1,:);
indNot1617=find(logicTableNot16(:,11)==0)
logicTableNot1617=logicTableNot16(indNot1617,:)

sum(logicTableNot1(:,9))./size(logicTableNot1,1)

indNot113=find(logicTableNot1(:,9)==0)
logicTableNot113=logicTableNot1(indNot113,:)