% Written by Vanja Dukic, 2000.
%
% Please cite: 
% Dukic and Gatsonis, “Meta-analysis of Diagnostic Test Accuracy Assessment Studies with Varying Number of Thresholds”, Biometrics, 59(4):936-946
% Any questions or comments, please contact me: vanja.dukic@colorado.edu
%
% Prof. Vanja Dukic
% Department of Applied Mathematics
% University of Colorado-Boulder
% Box 526; Boulder, CO 80309-0526
% tel: 303-735-0919; fax: 303-492-4066


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%        plot the outcome ROC on top of empirical points
%        a plot of ( 1 - P(Y<=j | no disease ) ) on X-axis (FP)
%        and ( 1 - P(Y<=j | disease ) ) on Y-axis (TP)


beta      =  -2.4668 ;
sebeta =     0.0890;
scale   =      0.5485;   
sescale=     0.0484;
sc1 = 0.1018;
sesc1 = 0.0261;


K=20;

figure
axis('square'); hold on;

for k=1:K,

eval([ 'fid=fopen(''s' num2str(k) ''' ); ']);
s=fscanf(fid, '%g', [2,inf]);
s=s';

      y=s(:,1);
      x=s(:,2);
    eval([ 'y' num2str(k) ' = y; ']);
    eval([ 'x' num2str(k) ' = x; ']);

    aa=[]; bb=[];
    aa = cumsum(crosstab(x+1, y),2);
    bb = (aa(:,size(aa,2))*ones(1,size(aa,2))); 
    bb = ((bb-aa)./bb); 
    bb =  [ [1;1], bb ]; 

if ismember(k,[3 4 5 6 7 10 11 12 15 17 18 19 20]), %prosp
    plot(bb(1,:),bb(2,:),'O'); 
else, 
    plot(bb(1,:),bb(2,:),'O'); 
    plot(bb(1,:),bb(2,:),'*'); 
end;
    plot(bb(1,:),bb(2,:),':');


end;


th=-10:0.05:10;
    P1=1-cdf('norm',(th+beta)/exp(scale+sc1),0,1);
    P2=1-cdf('norm',th,0,1);
    plot(P2,P1,'-');
 title('Pregnancy Failure Studies: \newline sROC curves and 95% confidence bands');

%...the area under ROC 
    A=(P2(2:length(P2))-P2(1:(length(P2)-1))).*(P1(2:length(P1))+P1(1:(length(P1)-1)));
    A=-A/2;
    area=sum(A);

   sA=[];
    sAb=normpdf(beta/(sqrt(1+exp(2*scale+2*sc1)))) * (1/(sqrt(1+exp(2*scale+2*sc1))));
    sAa=(-sAb*beta*exp(2*scale+2*sc1))/(1+exp(2*scale+2*sc1));
    sAp=[sAb, sAa, sAa];
    cv=cmatrix((cn-2):cn, (cn-2):cn);

    sdAp=sqrt((sAp*cv*sAp'));

      eval(['text(.4,.35,''location ='  num2str(-beta, '%2.3f') ' (' num2str(sebeta, '%2.3f') ') \newline scale = '  num2str(scale, '%2.3f') ' (' num2str(sqrt(sescale.^2), '%2.3f') ') \newline design effect= '  num2str(sc1, '%2.3f') ' (' num2str(sqrt(sesc1.^2), '%2.3f') ')  '' ) ']);


    xlabel('False Positive Ratio    P( T^{+}  |  D^{-} )');
    ylabel('True Positive Ratio    P( T^{+}  |  D^{+} )');

% delta method bounds

sdelta=[] ;      
sdelta=normpdf((th+beta)/exp(scale+sc1))*exp(-scale-sc1);
delta=[sdelta', sdelta'.*(th'+beta), sdelta'.*(th'+beta)];
    load cmatrix;

cv=cmatrix((cn-2):cn, (cn-2):cn);

sdROC=sqrt(diag(delta*cv*delta'))

ub=1.96*sdROC';
lb=-1.96*sdROC' ;
hold on;
P1u=ub+P1;
P1l=lb+P1;
%plot(P2, P1u, '-')
%plot(P2, P1l, '-')
    

%...the area under ROC 
    A=(P2(2:length(P2))-P2(1:(length(P2)-1))).*(P1u(2:length(P1))+P1u(1:(length(P1)-1)));
    A=-A/2;
    Uarea=sum(A);
 
%...the area under ROC 
    A=(P2(2:length(P2))-P2(1:(length(P2)-1))).*(P1l(2:length(P1))+P1l(1:(length(P1)-1)));
    A=-A/2;
    Larea=sum(A);
 
eval(['text(.05,.15,'' Prospective sROC area = '  num2str(area, '%2.3f') '  '' ) ']);


    P1=1-cdf('norm',(th+beta)/exp(scale-sc1),0,1);
    P2=1-cdf('norm',th,0,1);
    plot(P2,P1, '.-');
 
%...the area under ROC 
    A=(P2(2:length(P2))-P2(1:(length(P2)-1))).*(P1(2:length(P1))+P1(1:(length(P1)-1)));
    A=-A/2;
    area=sum(A);

   sA=[];
    sAb=normpdf(beta/(sqrt(1+exp(2*scale-2*sc1)))) * (1/(sqrt(1+exp(2*scale-2*sc1))));
    sAa=(-sAb*beta*exp(2*scale-2*sc1))/(1+exp(2*scale-2*sc1));
    sAr=[sAb, sAa, sAa];
    cv=cmatrix((cn-2):cn, (cn-2):cn);

    sdAr=sqrt((sAr*cv*sAr'));
sdAd=sqrt((sAr-sAp)*cv*(sAr'-sAp'))

eval(['text(.05,.07,'' Retrospective sROC area = '  num2str(area, '%2.3f') '  '' ) ']);


sdelta=[] ;      
sdelta=normpdf((th+beta)/exp(scale-sc1))*exp(-scale+sc1);
delta=[sdelta', sdelta'.*(th'+beta), sdelta'.*(th'+beta)];

cv=cmatrix((cn-2):cn, (cn-2):cn);

sdROC=sqrt(diag(delta*cv*delta'));

ub=1.96*sdROC';
lb=-1.96*sdROC';
hold on;
P1u=ub+P1;
P1l=lb+P1;
%plot(P2, P1u, '-')
%plot(P2, P1l, '-')
    

%...the area under ROC 
    A=(P2(2:length(P2))-P2(1:(length(P2)-1))).*(P1u(2:length(P1))+P1u(1:(length(P1)-1)));
    A=-A/2;
    Uarea=sum(A);
 
%...the area under ROC 
    A=(P2(2:length(P2))-P2(1:(length(P2)-1))).*(P1l(2:length(P1))+P1l(1:(length(P1)-1)));
    A=-A/2;
    Larea=sum(A);
 

 ylim([0,1]);
