% 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


% input 20 non-viable pregnancy studies (s):
%---------------------------------------
% 3 4 5 6 7 8 9 10 12 13 14 17 18 20 21 22 23 24 25 27
%-----------------------------------------------------
% OR input 16 ectopic studies (e):
%---------------------------------------
% 3 4 5 6 7 8 9 10 12 13 14 17 18 20 21 22 
%-----------------------------------------------------


dl=[];dlt=[]; d2l=[]; d2lt=[]; db=0; dbb=0; db2=[]; tb=[]; tbb=[]; 
dev=0; L=[]; YY=[]; XX=[]; Z=[]; Z1=[];

stepr= 0.0001; tol =1e-4; incr =10; decr =2; nt=1;

pr=[-1 -1  1  1  1  1  1 -1 -1  1  1  1  -1 -1  1 -1  1  1  1  1];
% ECT: pr=[-1 -1  1  1  1  1  1 -1 -1  1  1  1  -1 -1  1 -1 ];

global Link;  % (logit = 1, probit = 2)
Link=2;

K=20;

snn=[]; 

for k=1:K,
eval([ 'fid=fopen(''s' num2str(k) ''' ); ']);
s=fscanf(fid, '%g', [2,inf]);  
s=s';
[sm,sn]=size(s);
snn(k)=range(s(:,1));  

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

[my(k)  ny(k)] = size(y);
[mx(k)  nx(k)] = size(x);


aa=[]; bb=[];
crosstab(x +1, y)
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 ]; 


% initial calculations
ymin(k)=min(y);
ymax(k)=max(y);
yrange(k) = ymax(k) - ymin(k);

np(k)=yrange(k)+3;	% no of parms

z = (y*ones(1,yrange(k)) )==( ones(my(k),ny(k))*( ymin(k):(ymax(k)-1)));
z = z(:,any(z));
Z= ( y*ones(1,yrange(k)) )==( ones(my(k),ny(k))*((ymin(k)+1): ymax(k)));
Z = Z(:,any(Z));
[mz(k) nz(k)]=size(z);
[mZ(k) nZ(k)]=size(Z);
eval([ 'z' num2str(k) ' = z; ']);
eval([ 'Z' num2str(k) ' = Z; ']);  % <-- WAS BUG: z1_1 overwrote z_11 !!!



load theta_pf;  % Again, this is a set of initial values for all cutpoints I found
if k==1,
theta=rka(1:snn(1));
eval([' theta' num2str(k) '=theta; ']);
else
theta=rka(sum(snn(1:(k-1)))+1:sum(snn(1:k)));
eval([' theta' num2str(k) '=theta; ']);
end;


tbt=[];
tbt = [theta; beta; scale; sc1];   %  all parms for k-th study
eval([ 'tb' num2str(k) '(:,nt) = tbt; ']);
tbb=[tbb;theta];		   %  all cutpoints


% likelihood  at the starting values

[Lk(k), devk(k)] = lik(y,x,x,tbt,pr(k));


% derivatives  at the starting values

dk=[]; d2k=[];
[dk,d2k] = derivs(y,x,x,tbt,pr(k));
eval([ 'dlk' num2str(k)  '  =  dk; ' ]);
eval([ 'd2lk' num2str(k) '  =  d2k; ' ]);


% first derivatives

dlt = [dlt;dk(1:nz(k))];		% d wrt cutpoints

db = db + dk( (nz(k)+1) : (np(k)) );	% d wrt loc and scale 

% second derivatives

zert1=[]; zert2=[]; zert3=[];
zert1 = zeros( size(d2lt,1), yrange(k) ); 
zert2 = zeros( yrange(k), size(d2lt,2) );
zert3 = d2k(1:yrange(k),1:yrange(k));

d2lt = [ [d2lt, zert1]; [zert2, zert3]  ];  % d2 wrt cut/cut


zert4=[]; 
zert4 = d2k( (yrange(k)+1): np(k) , 1:yrange(k) );

db2 = [db2, zert4];	    % d2 wrt (beta, scale, sc1)/cut


zert5=[]; 
zert5 = d2k( (yrange(k)+1): np(k) ,  (nz(k)+1):np(k) );

dbb = dbb + zert5;	% d2 wrt (b, sc, sc1)/(b, sc, sc1)

end;


tb(:,nt)=[tbb;beta; scale; sc1];
L=sum(Lk); dev=sum(devk);

dl  = [ dlt;  db ];
d2l = [ [d2lt, db2'] ; [db2 , dbb] ];


% The commands below are going to save all derivatives and 
% the likelihood in the current directory
save dl dl
save d2l d2l
save L L
