% 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%
% Code to tidy up the output
   for k=1:K,
     eval([ 'if k==1, theta' num2str(k) ' = tb( 1 : nzcum(k),nt); else theta' num2str(k) ' = tb( (nzcum(k-1)+1) : nzcum(k),nt);  end; ']);
   end;

beta=tb((nzcum(K)+1):(nzcum(K)+1),nt);     
scale=tb((nzcum(K)+2):(nzcum(K)+2),nt); 
sc1 = tb(size(tb,1), nt);

cmatrix= (inv(-d2l));
save cmatrix cmatrix;
display('ML based cov matrix of Loc, Scale and Scale Covariate')
cn=size(cmatrix,1);
disp(cmatrix((cn-2):cn , (cn-2):cn))

   se=sqrt(diag(inv(-d2l)));
   st=length(se);
   sebeta=se((nzcum(K)+1): (st-2) );
   sescale=se((nzcum(K)+nx(1)+1):(st-1));
   sesc1 = se(st);

   disp('Number of iterations'); disp(nt);
   disp('  PROBIT LINK'); 
   disp('  Deviance'); disp(dev);

for k=1:K,
   eval([ 'disp(''  Study ' num2str(k) '''); ']);
   eval([ 'disp(''     Theta     SE''); ']);
   theta=[];
   eval([ 'theta=(theta' num2str(k) '); ']);
   lt=length(theta);
   disp([theta   se(1:lt)]);
end;

      disp('     Beta       SE');
      disp([beta sebeta ]);
        disp('     Scale       SE');
      disp([scale sescale ]);
      disp([sc1 sesc1 ]);
  

for k=1:K,
      theta=[]; x=[]; e=[];
      eval([ 'theta=(theta' num2str(k) '); ']);
      eval([ 'x=x' num2str(k) '; ']);
      eval([ 'y=y' num2str(k) '; ']);
 if ismember(k,[3 4 5 6 7 10 11 12 15 17 18 19 20]),
      e=( (x*(beta))*ones(1,nz(k)) )+( ones(my(k),ny(k))*theta' );
      ex=exp( (x*(scale+sc1)) * ones(1,nz(k)) );
 else,
      e=( (x*(beta))*ones(1,nz(k)) )+( ones(my(k),ny(k))*theta' );
      ex=exp( (x*(-sc1+scale)) * ones(1,nz(k)) );
 end;
      e=e./ex;   

   p=[];
   p=diff([zeros(my(k),ny(k))   cdf('norm',e,0,1) ones(my(k),ny(k))]')';
   disp('  ')
   disp('  ')
   eval([ 'disp(''  P(Y_i=j) in Study ' num2str(k) '''); ']);
   disp([['   disease   ' ; 'no disease   '], num2str([unique(p,'rows')])]);
   eval([ 'disp(''  Predicted Counts in Each Category in Study ' num2str(k) ':''); ']);
   disp(sum(p))
   eval([ 'disp(''  Observed Counts in Each Category in Study ' num2str(k) ':''); ']);
   eval([' disp( sum( crosstab(x' num2str(k) '+1 , y' num2str(k) ') )) ']);

end;

