%% analysis of cytotoxicity data (Section 5.2)
%  May 22, 2014
%  IMPORTANT: the analysis requires the imaging toolbox

%% set preferences

iptsetpref('ImshowAxesVisible', 'on');


%% DATA AND OTHER DEFINITIONS

% output from logistic regression

% fitted beta values
betahat     =   [-1.330 -0.084 0.159 0.003883 -0.001308];

% Sigma estimate (in the notation of the paper this is hat(Sigma_n)/n)
Sigma       =   zeros(5,5);
Sigma(1,:)  =   [0.01389 -0.00128 -0.001079 0.00003679 0.000009878];
Sigma(2,:)  =   [-0.00128 0.0007083 -0.00002017 -0.00002673 0.0000001613];
Sigma(3,:)  =   [-0.001079 -0.00002017 0.0002672 0.0000009442 -0.000002565];
Sigma(4,:)  =   [0.00003679 -0.00002673 0.0000009442 0.000001063 -0.000000007581];
Sigma(5,:)  =   [0.000009878 0.0000001613 -0.000002565 -0.000000007581 0.00000002485];


% the data itself
MMS         = [0 0 0 0 10 100 250 10 10 10 100 100 100 250 250 250]/10;
PMA         = [0 1 10 100 0 0 0 1 10 100 1 10 100 1 10 100];
dead        = [19 24 54 68 16 17 19 19 41 79 10 36 62 36 56 74];
total       = [98 87 91 87 91 90 92 94 90 93 83 85 83 92 93 87];

% other definitions 
nlat        =  401; % number of pixels is nlat^2
x           =  0:26/(nlat-1):26;
y           =  0:130/(nlat-1):130;
[XX,YY]     =  meshgrid(x, fliplr(y));


%% calculation and plot of ED50 estimate 

Fhat        =  betahat(1)+betahat(2).*XX+betahat(3).*YY+betahat(4).*XX.^2+betahat(5).*YY.^2;  % fitted response curve
eta0        =  0;
Ihat        =  Fhat >= eta0;  % ED50_plus estimate


imshow(Ihat+1-Ihat);
hold on
[cc,hh1]    = imcontour(Ihat,[0 0],'k'); % plot of ED50 estimate
hold off
set(hh1, 'LineWidth', 1.25) 
set(gcf,'Color', 'w');
set(get(gca,'XLabel'),'String','MMS/10');
set(get(gca,'YLabel'),'String','PMA');
set(gca,'XTick',401*([5:5:25])/26)
set(gca,'XTickLabel',[5:5:25])
set(gca,'YTick',401*([20:20:120])/130)
set(gca,'YTickLabel',[120:-20:20])



%% SCHEFFE'S BOUNDS (result is plotted)

V           =  Sigma;
Vhat        =  V(1,1)+2*V(1,2)*XX   +2*V(1,3)*YY    +2*V(1,4).*XX.^2   +2*V(1,5).*YY.^2 ...
                     +  V(2,2)*XX.^2+2*V(2,3)*XX.*YY+2*V(2,4)*XX.^3    +2*V(2,5).*XX.*YY.^2 ...
                                    +  V(3,3).*YY.^2+2*V(3,4)*XX.^2.*YY+2*V(3,5).*YY.^3 ...
                                                    +  V(4,4)*XX.^4    +2*V(4,5).*XX.^2.*YY.^2 ...
                                                                       +  V(5,5).*YY.^4;
% calculate 95% confidence region
chialpha    =  11.07050;
CR_SH_up    =  Fhat + sqrt(chialpha*Vhat) >= eta0; 
CR_SH_down  =  Fhat - sqrt(chialpha*Vhat) <= eta0;
CR_SH       =  CR_SH_up.*CR_SH_down;

% plot region along with ED50 estimate
[ii] = imshow(1-0.3*CR_SH);
hold on
[cc,hh1]    = imcontour(Ihat,[0 0],'k'); 
hold off
set(hh1, 'LineWidth', 1.25)
set(gcf,'Color', 'w');
set(get(gca,'XLabel'),'String','MMS/10');
set(get(gca,'YLabel'),'String','PMA');
set(gca,'XTick',401*([5:5:25])/26)
set(gca,'XTickLabel',[5:5:25])
set(gca,'YTick',401*([20:20:120])/130)
set(gca,'YTickLabel',[120:-20:20])




%%  NEW METHOD (takes about one minute, result is plotted)

B           =  1000;
A           =  chol(Sigma);
useseed     =  43814;
q           =  0.975;  %1-alpha/2


count1      =  zeros([nlat nlat]);
count2      =  zeros([nlat nlat]);
EMPcount1   =  zeros(1,B);
EMPcount2   =  zeros(1,B);

% find the intersection of nonempty sets
RandStream.setGlobalStream(RandStream('mt19937ar','seed', useseed));

for i = 1:B
    
    betastar    =  betahat + (A'*randn(5, 1))';
    Fstar       =  betastar(1)+betastar(2).*XX+betastar(3).*YY+betastar(4).*XX.^2+betastar(5).*YY.^2;
    Istar       =  Fstar >= eta0;
    EMPcount1(i)=  1-max(max(Istar));
    EMPcount2(i)=  1-max(max(~Istar));
    count1      =  count1+Istar;
    count2      =  count2+(~Istar);
        
end;

K1          =  count1 >= (B-sum(EMPcount1));
K2          =  count2 >= (B-sum(EMPcount2));
df1         =  bwdist(K1);
df2         =  bwdist(K2);
rank1       =  zeros(1,B);
rank2       =  zeros(1,B);


% calculate distances from all empty sets
RandStream.setGlobalStream(RandStream('mt19937ar','seed', useseed));

for i = 1:B
    
    betastar    =  betahat + (A'*randn(5, 1))';
    Fstar       =  betastar(1)+betastar(2).*XX+betastar(3).*YY+betastar(4).*XX.^2+betastar(5).*YY.^2;
    Istar       =  Fstar >= eta0;
    
    if ~isempty(find(Istar>0))
        rank1(i)    =  max(df1(Istar));
    end;
    
    if ~isempty(find(~Istar>0))
        rank2(i)    =  max(df2(~Istar));    
    end;   
        
end;

% find the cutoff value
ind1        =  find(EMPcount1==0);
ind2        =  find(EMPcount2==0);
rank11      =  rank1(ind1);
rank22      =  rank2(ind2);
sortrank1   =  sort(rank11);
sortrank2   =  sort(rank22);


cutoffInd1  =  floor(q*length(sortrank1));
cutoff1_m2  =  sortrank1(cutoffInd1);
cutoffInd2  =  floor(q*length(sortrank2));
cutoff2_m2  =  sortrank2(cutoffInd2);

count1_m2   =  zeros([nlat nlat]);
count2_m2   =  zeros([nlat nlat]);


% find union of all sets within cutoff
RandStream.setGlobalStream(RandStream('mt19937ar','seed',43814));

for i = 1:B
    
    betastar    =  betahat + (A'*randn(5, 1))';
    Fstar       =  betastar(1)+betastar(2).*XX+betastar(3).*YY+betastar(4).*XX.^2+betastar(5).*YY.^2;
    Istar       =  Fstar >= eta0;
    
    if(rank1(i)<=cutoff1_m2)
        count1_m2      =  count1_m2+Istar;
    end;
   
    if(rank(i)<=cutoff2_m2)
        count2_m2      =  count2_m2+(~Istar);
    end;
    
end;



% calculate 95% confidence region
CR_up       =  count1_m2>0;
CR_down     =  count2_m2>0;
CR          =  CR_up.*CR_down;


% plot region along with ED50 estimate
imshow(1-0.3*CR);
hold on
[cc,hh1]    = imcontour(Ihat,[0 0],'k'); 
hold off
set(hh1, 'LineWidth', 1.25)
set(gcf,'Color', 'w');
set(get(gca,'XLabel'),'String','MMS/10');
set(get(gca,'YLabel'),'String','PMA');
set(gca,'XTick',401*([5:5:25])/26)
set(gca,'XTickLabel',[5:5:25])
set(gca,'YTick',401*([20:20:120])/130)
set(gca,'YTickLabel',[120:-20:20])