/****************************************************************************************

  Macro SIZE computes sample size and power for clinical trails comparing a treatment
  and control group, incorporating drop-ins, drop-outs, losses to followup and other
  issues.  There are many parameters and options that allow the user to do many things.
  These are described in "Sample Size Estimation for Complex Clinical Trials" by
  Joanna Shih and outlined below. 
  An example is given using the most used parameters.

Required Arguments:
 np       = number of periods
 pc       = event rates in control group
 pcratio  = ratio of hazards (relative to first period) in control group over the np periods
 k        = hazard ratio (treatment/control)

1 of the following:
 n        = total sample size  OR
 power    = power


Optional Arguments: (Inclusion of some of these may make other parameters required):
 alpha	  = type I error for two-tailed alternatives. The default value is 0.05.
 din      = dropin rates (proportion in control group that switch to treatment group. 
	    default is 0)
 diratio  = ratio of hazard rates (relative to first period) of dropins over the np periods)
 distdin  = a discrete prior probability distribution for din.
 distdout = a discrete prior probability distribution for dout.
 distk    = a discrete prior probability distribution for the hazard ratio, k. Each k value 
	    is considered a mass point and is assigned a probability mass from distk. The 
	    number of masses in disk must be equal to the number of entries in k. Values of 
 	    distk are separated by spaces; the sum of values is 1.
 distloss = a discrete prior distribution for loss.
 distpc	  = a discrete prior probability distribution for pc. Each vector of pc is considered 
	    a mass point and is assigned a probability mass. The number of masses in distpc 
	    must be equal to the number of vectors in pc. Values of distpc are separated by 
	    spaces; the sum of values is 1.
 doratio  = ratio of hazards (relative to first period) of dropouts over the np periods).
 dout     = dropout rates (proportion in treatment group that switch to control group)
 	    The default is 0.
 lag 	  = lag time, the number of periods to achieve the full treatment effect. It adjusts 
	    the event rate for crossovers under the nonproportional hazards model. If lag is 
	    not specified, then participants return to the efficacy level comparable to the 
	    level in the opposite arm immediately after crossover.
 lagdout  = lagged effect of dropout. It equals 1 if crossovers from the treatment group 
	    return to the control group level in the same fashion as they reach the treatment 
	    effect level and equals 0 if the treatment effect vanishes immediately after 
	    crossoever. The default value is 1.
 logr 	  = indicator of the test statistic. If logr = 0, then the binomial formula is used. 
	    If logr = 1, then the log-rank test is used. If logr = a vector of np weights, 
	    then the weighted log-rank test is used. The default value is 1. 
 loratio  = ratio of hazard rates(relative to first period) of lost to followups over np 
	    periods
 loss     = proportion of participants lost to followup because of censoring
 nmin,nmax,ntol
 	  = minimum and maximum of the sample size to be used in the iterative bisection 
	    method to search for the sample size such that the calculated power is equal
	    to or greater than the specified power within the relative tolerance, ntol.
	    Default is 0, 20000, 0.001.
 npmin,npmax,nptol
	  = minimum and mzximum of the duration of time to be used in the iterative 
	    bisection method to search for the duration time such that the calculated power 
	    is equal to or greater than the specified power within the relative tolerance, 
	    nptol. Default is 0, 200, 0,001.
 output   = name of the unix file which stores results. Can then be used for futher
            analyses or plotting.
 prop	  = 0 or 1. If prop = 1 then the hazards are proportional and k remains the same 
	    for the duration of the study. If proportion = 0, then the hazards are 
	    nonproportional and k must be a vector of np entries.
 ratio	  = two numbers separated by spaces, referring to the ratio of the sample size in 
	    the control group to that in the treatment group. For example, if the allocation 
	    ratio is 2:3, then ratio = 2 3.
 rectime  = list of recruitment (accrual) periods (if simult = 0, then must be specified.)
 recratio = list of relative (to first period) recruitment rates over recruitment periods.
	    must be the same length as rectime.
 recrate  = the recruitment rate in each period. If both sample size and power are specified 
	    and the goal is to solve for the duration time, then recrate must be specified.
 sbdv	  = number of intervals to be divided in each period. Within a given period, the 
	    probability of an event is assumed the same across intervals of equal length
	    and equals 1 - (1-x)^(1/sbbdv), where x is the probability of an event in a 
	    period. If lag is not specified, then the default of sbdv is 20. If lag is 
	    specified, the default is 4. The product of lag and sbdv determines the number 
	    of intermediate states for the event rate in the treatment groups. The 
	    intermediate states account for the lag in treatment effect.
 simult   = set to 0 if participants are recruited over periods; set to 1 if participants
            enter simultaneously, if 0, rectime and recratio must be specified.
 printpar = if set to something then macro parameters are printed to listing file


EXAMPLE OF TYPICAL CALL

%size (
np = 6,                               [6 periods (years/months etc)]
pc = 0.15 5,                          [event rate of 15% in control group after 5 periods]
pcratio = 1 1.06 1.12 1.19 1.26 1.36, [event rate increases by 6% each period]
k = 0.84,                             [hazard rate treatment/control]
power = 0.80,                         [power is set so sample size will be calculated]
loss = 0.10 5,                        [lost to follow of 10% after 5 periods]
loratio = 1 1 1 1 1 1,                [lost to followup rate equal each period (2%)]
din = 0.20 5,                         [dropin rate of 20% after 5 years]
diratio = 1 1 1 2 2 2,                [dropin rate doubles last 3 periods]
dout = 0.20 5,                        [dropout rate of 20% after 5 periods]
doration = 2 1 1 1 1 1,               [dropout rate double in first period]
rectime = 1 2,                        [two periods to complete recruitment]
recratio = 1 1,                       [equal number recruited in 2 periods]
simult = 0 );                         [must be set to 0 when their is staggered entry]

******************************************************************************************/

options mprint;
%macro size(np=0,
            pc=0,
            pcratio =0,
            k= ,
            n= ,
            power= ,
            din=0,
            dout=0,
            loss=0,
            sbdv=20,
            alpha = .05,
            simult=1,
            rectime=0,
            recratio=0,
            recrate=0,
            ratio=1,
            prop = 1,
            lag = 0,
            distpc = 0,
            distk = 0,
            distdin = 0,
            distdout = 0,
            distloss = 0,
            diratio=0,
            doratio=0,
            loratio=0,
            nmin = 0,
            nmax = 20000,
            ntol = .001,
            npmin = 0,
            npmax = 200,
            nptol = .01,
            itmax = 9999,
            output = ,
            printpar =,
            logr = 1);
%lakprog;
%write;
%mend size;

%macro lakprog;
proc iml;
*reset noprint;
reset fw=7;

periods = &np;
sbdv = &sbdv;
alpha = α
pc = {&pc};
noncmpl = {&dout};
dropin = {&din};
loss = {&loss};
simult = &simult;
rectime = {&rectime};
recratio = {&recratio};
recrate = {&recrate};
ratio = {&ratio};
pcratio = {&pcratio};
diratio = {&diratio};
doratio = {&doratio};
loratio = {&loratio};
itmax = &itmax;
lag = &lag;
prop = ∝

%if &printpar ne %then %do;
print " np=&np, pc=&pc, pcratio=&pcratio, k=&k, n=&n, power=&power";
print " din=&din, diratio=&diratio, dout=&dout, doratio=&doratio" ;
print " loss=&loss, loratio=&loratio, rectime=&rectime recratio=&recratio, simult=&simult";
print " ratio=&ratio";
%end;

if periods = 0 then do;
 print " No. of periods were not specified.";
 stop;
end;

if (pc = 0) then do;
 print " Baseline event rates were not specified.";
 stop;
end;

if ratio = 1 then do; r1=1; r2=1; end;
else do; r1=ratio[,1]; r2=ratio[,2]; end;

if (pcratio = 0) then do;
 if(mod(ncol(pc),periods) ^= 0) then do;
  print " The length of event rates does not match no. of periods.";
  stop;
 end;
end;
if (pcratio ^= 0) then do;
 if (pc = 0) then do;
  print " event rates are not specified.";
  stop;
 end;
 if(mod(ncol(pcratio),periods) ^= 0 | mod(ncol(pc),2) ^= 0) then do;
  print "The length of event rates does not match no. of periods or";
  print "pc is inadequately specified.";
  stop;
 end;
end;

if (diratio = 0) then do;
 if(dropin ^= 0 & mod(ncol(dropin),periods) ^= 0) then do;
  print "The length of dropins does not match no. of periods.";
  stop;
 end;
end;
if (diratio ^= 0) then do;
 if (dropin = 0) then do;
  print "Dropin rates are not specified.";
  stop;
 end;
 if(mod(ncol(diratio),periods) ^= 0 | mod(ncol(dropin),2) ^= 0)then do;
  print "The length of dropins does not match no. of periods or";
  print "din is inadequately specified.";
  stop;
 end;
end;

if (doratio = 0) then do;
 if(noncmpl ^= 0 & mod(ncol(noncmpl),periods) ^= 0) then do;
  print "The length of dropouts does not match no. of periods.";
  stop;
 end;
end;
if (doratio ^= 0) then do;
 if (noncmpl = 0) then do;
  print "dropout rates are not specified.";
  stop;
 end;
 if(mod(ncol(doratio),periods) ^= 0 | mod(ncol(noncmpl),2) ^= 0)
  then do;
  print "The length of dropout rates does not match no. of periods or";
  print "dout is inadequately specified.";
  stop;
 end;
end;

if (loratio = 0) then do;
 if(loss ^= 0 & mod(ncol(loss),periods) ^= 0) then do;
  print "The length of loss of follow-up does not match no. of periods.";
  stop;
 end;
end;
if (loratio ^= 0) then do;
 if (loss = 0) then do;
  print "loss of follow-up rates are not specified.";
  stop;
 end;
 if(mod(ncol(loratio),periods) ^= 0 | mod(ncol(loss),2) ^= 0) then do;
  print "The length of loss does not match no. of periods or";
  print "loss is inadequately specified.";
  stop;
 end;
end;

%if %length(&power) ^= 0 %then %do;
 power = 999;
%end;
%else %do;
 power = .;
%end;

%if %length(&k) ^= 0 %then %do;
 k = 999;
%end;
%else %do;
 k = .;
 print "k was not specified.";
 stop;
%end;

%if %length(&n) ^= 0 %then %do;
 ntot = 999;
%end;
%else %do;
 ntot = .;
%end;

if power = . & ntot = . then do;
 print "Both power & sample size are unspecified.";
 stop;
end;

if (sbdv = 20 & (prop ^= 1 & lag ^= 0)) then sbdv = 4;
/* check k, pc, din, dout, loss are deterministic or from a distribution.
   d1 = distn for k; d2 = distn for pc; d3 = distn for dropins
   d4 = distn for dropouts; d5 = distn for loss to follow-up */
%let d1 = &distk; %let d2= &distpc; %let d3= &distdin; %let d4=&distdout;
%let d5 = &distloss;
d1 = {&d1}; d2={&d2}; d3={&d3}; d4={&d4}; d5={&d5};

/* max no. of masses */
m1max = 1000; m2max = 100; m3max = 100; m4max = 100; m5max = 100;
/* masses */
m1 = j(m1max,1,0); m2 = j(m2max,1,0); m3 = j(m3max,1,0);
m4=j(m4max,1,0); m5 = j(m5max,1,0);
/* no. of masses */
nom1 = 0; nom2 = 0; nom3 = 0; nom4 = 0; nom5 = 0;

%do ii = 1 %to 5;
%if &&d&ii ^= 0 %then %do;
 %let to = %index(%upcase(&&d&ii),TO);
 %if &to ^= 0 %then %do;
  %let len = %eval(&to-1);
  minm = %substr(&&d&ii,1,&len);
  %let by = %index(%upcase(&&d&ii),BY);
  %let len2 = %eval(&by - &to - 2);
  maxm = %substr(&&d&ii,%eval(&to+2),&len2);
  %let len3 = %length(&&d&ii);
  %let len3 = %eval(&len3 - &by-1);
  inc = %substr(&&d&ii,%eval(&by+2),&len3);
  nom&ii = int((maxm-minm)/inc)+1;
  m&ii = j(nom&ii,1,0); m&ii[1,] = minm;
  m&ii[2:nom&ii,]=j(nom&ii-1,1,inc);
 %end;
 %else %do;
 %let jj = 1;
 %mloop:
 %let value = %scan(&&d&ii,&jj,' ' );
 %if %length(&value) ^= 0  %then %do;
  m&ii[&jj,] = &value;
  %let jj = %eval(&jj+1);
  %goto mloop;
 %end;
 %else %do;
  nom&ii = %eval(&jj-1);
 %end;
 %end;
%end;
%end; /*ii*/

/*Scan k, n, power, pc, dropins, dropouts & loss to
  follow-up rates */
%let v1 = &k; %let v2 = &pc; %let v3 = &din; %let v4 = &dout;
%let v5 = &loss; %let v6 = &n; %let v7 = &power;
%let u2 = &pcratio; %let u3 = &diratio; %let u4 = &doratio;
%let u5 = &loratio;

/* max no. of mass points */
mp1max = 1000; mp3max = 1000; mp4max=1000; mp5max=1000;
mp6max = 1000; mp7max = 1000; mp2max = 1000;
/* mass points */
mp1 = j(mp1max,1,0); mp2=j(mp2max,1,0); mp3=j(mp3max,1,0);
mp4 = j(mp4max,1,0); mp5=j(mp5max,1,0); mp6=j(mp6max,1,0);
mp7 = j(mp7max,1,0);
/* initialize no. of mass points */
nomp1 = 0; nomp2 = 0; nomp3=0; nomp4=0;
nomp5=0; nomp6=0; nomp7=0;

%do ii = 1 %to 3; /* 1 = k, 2= n, 3 = power */
 %let kk = ⅈ
 %if &ii > 1 %then %let kk = %eval(&ii+4);
 %if &kk = 6 & &v6 =  %then %goto mplp;  /*n is not specified.*/
 %if &kk = 7 & &v7 =  %then %goto mplp; /*power is not specified. */
 %let to = %index(%upcase(&&v&kk),TO);
 %if &to ^= 0 %then %do;
  %let len = %eval(&to-1);
  minm = %substr(&&v&kk,1,&len);
  %let by = %index(%upcase(&&v&kk),BY);
  %let len2 = %eval(&by - &to - 2);
  maxm = %substr(&&v&kk,%eval(&to+2),&len2);
  %let len3 = %length(&&v&kk);
  %let len3 = %eval(&len3 - &by-1);
  inc = %substr(&&v&kk,%eval(&by+2),&len3);
  nomp&kk = int((maxm-minm)/inc)+1;
    mp&kk = j(nomp&kk,1,0); mp&kk[1,] = minm;
  do jj = 2 to nomp&kk;
   mp&kk[jj,] = mp&kk[jj-1,]+inc;
  end ;
 %end;
 %else %do;
  %let jj = 1;
  %mploop:
  %let value = %scan(&&v&kk,&jj,' ' );
  %if %length(&value) ^= 0 %then %do;
   mp&kk[&jj,] = &value;
   %let jj = %eval(&jj+1);
   %goto mploop;
  %end;
  %else %do;
   %put nomp %eval(&jj-1);
   nomp&kk = %eval(&jj-1);
  %end;
 %end;
%if &ii = 1 & &prop = 0 %then %do;
  if nomp1 ^= periods then do;
   print "number of hazard ratios & number of periods are not equal.";
   print nomp1;
   print periods;
   stop;
  end;
  else do;
   nomp1 = 1;
  end;
%end;

%mplp:
%end; /*ii*/

 if d1 ^= 0 & nomp1 ^= nom1 then do;
  print " number of masses & mass points for K are not equal.";
  print nom1;
  print nomp1;
  stop;
 end;

%do ii = 2 %to 5;
tempv = {&&v&ii};
tempu = {&&u&ii};
%if &&v&ii = 0 %then %do;
 nomp&ii = 1; mp&ii = j(1,periods,0);
%end;
%else %do;
 %if &&u&ii = 0 %then %do;
  nomp&ii = ncol(tempv)/periods;
  mp&ii = j(nomp&ii,periods,0);
  do jj = 1 to nompⅈ
  mp&ii[jj,] = tempv[,(jj-1)*periods+1:jj*periods];
  end;
  if nomp&ii ^= nom&ii & d&ii ^= 0 then do;
   print " number of masses & mass points for rates are not equal.";
   print nomⅈ print nompⅈ
   stop;
  end;
 %end;
 %else %do;
  no1= ncol(tempv)/2;
  no2= ncol(tempu)/periods;
  nomp&ii= no1#no2;
  mp&ii = j(nomp&ii,periods,0);
  nn = 1;
  do kk = 1 to no1;
  do ll = 1 to no2;
   sum= 0;
   do mm = 1 to tempv[,2#kk];
    sum= sum+tempu[,(ll-1)#periods+mm];
   end;
   rate = -log(1-tempv[,2#(kk-1)+1])/sum;
   do mm = 1 to periods;
    mp&ii[nn,mm] = 1-exp(-rate#tempu[,(ll-1)#periods+mm]);
   end;
   nn = nn + 1;
  end;
  end;
 %end;
%end;
%end;


/*Check the order of the DO loops to be executed below. Events with
a p.d.f. are of lower orders & are executed sooner in the
DO loops. */
%let high = 0; %let low = 6;
%do ii = 1 %to 5; /* 1=k 2=pc 3=din 4=dout 5=loss */
 %if &&d&ii = 0 %then %do; /*no p.d.f. */
  %let od&ii = %eval(&high+1);
  %let high = &&odⅈ
  %let odi&high = ⅈ
  %let random&high = 0;
 %end;
 %else %do; /* has a p.d.f. */
  %let od&ii = %eval(&low-1);
  %let low = &&odⅈ
  %let odi&low = ⅈ
  %let random&low = 1;
 %end;
%end;

%if &power = | &n = %then %do;
/* Adjust for staggered entry. */
n_intrvl=sbdv; /* # of subinterval in each period */
*THE NEXT 10 LINES SET UP A VECTOR FOR STAGGERED ENTRY;
nn=periods#n_intrvl; ad_cens=j(1,nn,0); rcrt_sum=0;
if simult=0 then do;
wk = 0||int(rectime#n_intrvl+0.5);
rcrt=j(1,nn,0);
do ii=1 to ncol(rectime);
 rcrt[1,(wk[1,ii]+1):wk[1,ii+1]]=j(1,wk[1,ii+1]-wk[1,ii],recratio[1,ii]);
   end;
do ii=1 to nn;
   rcrt_sum=rcrt_sum+rcrt[1,ii];
   ad_cens[1,nn+1-ii]=rcrt[1,ii]/rcrt_sum;
   end;
end; /* end staggered entry. */
%end;

markov:
free par2; /*used for storing results */
%if &power = %then %do;
do iv6 = 1 to nomp6;
 n=mp6[iv6,];
 n_lr=n;
  meanpow5 = 0; meandth5 = 0; meanpc5 = 0; meanpe5 = 0;
do iv1 = 1 to nomp&odi1;
 meanpow4 = 0; meandth4 = 0; meanpc4 = 0; meanpe4 = 0;
do iv2 = 1 to nomp&odi2;
 meanpow3 = 0; meandth3 = 0; meanpc3 = 0; meanpe3 = 0;
do iv3 = 1 to nomp&odi3;
 meanpow2 = 0; meandth2 = 0; meanpc2 = 0; meanpe2 = 0;
do iv4 = 1 to nomp&odi4;
 meanpow1 = 0; meandth1 = 0; meanpc1 = 0; meanpe1 = 0;
do iv5 = 1 to nomp&odi5;
 pc = mp2[iv&od2,];
 if &prop = 1 then
  k = mp1[iv&od1,];
 else
  k = mp1;
 dropin = mp3[iv&od3,];
 noncmpl = mp4[iv&od4,];
 loss = mp5[iv&od5,];
  temp1 = mp1[iv&od1,]; temp2 = iv&od2; temp3=iv&od3; temp4=iv&od4;
  temp5=iv&od5;
%markov; /* invoke markov chain model */
%if &logr = 1
  %then %do ;
    d_lr = n#(r1#pc2+r2#pe2)/(r1+r2);
    z_beta = (d_lr##0.5)#ed-z_alpha;
  %end ;
  %else %do ;
    p = (r1#pc2+r2#pe2)/(r1+r2);
    sd1 = ((r1+r2)#p#(1-p)/(r1#r2))##0.5;
    sd2 = (pc2#(1-pc2)/r1+pe2#(1-pe2)/r2)##0.5;
    z_beta = (((n##0.5)#abs((pc2-pe2))/((r1+r2)##0.5)-z_alpha#sd1)/sd2) ;
    d_lr = n#p ;
  %end ;
  d_lr = int(d_lr+0.5);
  power = probnorm(z_beta);

if &prop = 0 then temp1 = .;

if &random5 = 0 then do;
 n_lr = int(n+0.5); d_lr = int(d_lr+0.5);
 par = alpha||power||temp1||iv&od2||iv&od3||iv&od4||iv&od5||pc2||pe2||
      n_lr||d_lr||periods;
 par2 = par2//par;
end;
meanpow1 = meanpow1+power*m&odi5.[iv5,];
meandth1 = meandth1+d_lr*m&odi5.[iv5,];
meanpc1 = meanpc1+pc2*m&odi5.[iv5,];
meanpe1 = meanpe1+pe2*m&odi5.[iv5,];
end; /*v5 */

if &random4 = 0 & &random5 = 1 then do;
 temp&odi5 = .;
 n_lr = int(n+0.5); meandth1 = int(meandth1+0.5);
 par = alpha||meanpow1||temp1||temp2||temp3||temp4||temp5||pc2||pe2||
      n_lr||meandth1||periods;
 par2 = par2//par;
end;
meanpow2 = meanpow2+meanpow1*m&odi4.[iv4,];
meandth2 = meandth2+meandth1*m&odi4.[iv4,];
meanpc2 = meanpc2+meanpc1*m&odi4.[iv4,];
meanpe2 = meanpe2+meanpe1*m&odi4.[iv4,];
end; /*v4*/

if &random3 = 0 & &random4 = 1 then do;
 temp&odi5 = .; temp&odi4 = .;
 n_lr = int(n+0.5); meandth2 = int(meandth2+0.5);
 par = alpha||meanpow2||temp1||temp2||temp3||temp4||temp5||pc2||pe2||
      n_lr||meandth2||periods;
 par2 = par2//par;
end;
meanpow3 = meanpow3+meanpow2*m&odi3.[iv3,];
meandth3 = meandth3+meandth2*m&odi3.[iv3,];
meanpc3 = meanpc3+meanpc2*m&odi3.[iv3,];
meanpe3 = meanpe3+meanpe2*m&odi3.[iv3,];
end; /*v3*/

if &random2 = 0 & &random3 = 1 then do;
 temp&odi5 = .; temp&odi4 = .; temp&odi3 = .;
 n_lr = int(n+0.5); meandth3 = int(meandth3+0.5);
 par = alpha||meanpow3||temp1||temp2||temp3||temp4||temp5||pc2||pe2||
      n_lr||meandth3||periods;
 par2 = par2//par;
end;
meanpow4 = meanpow4 +meanpow3*m&odi2.[iv2,];
meandth4 = meandth4 +meandth3*m&odi2.[iv2,];
meanpc4 = meanpc4 +meanpc3*m&odi2.[iv2,];
meanpe4 = meanpe4 +meanpe3*m&odi2.[iv2,];
end; /*v2*/

if &random1 = 0 & &random2 = 1 then do;
 temp&odi5 = .; temp&odi4 = .; temp&odi3 = .; temp&odi2 = .;
 n = int(n+0.5); meandth4 = int(meandth4+0.5);
 par = alpha||meanpow4||temp1||temp2||temp3||temp4||temp5||pc2||pe2||
      n_lr||meandth4||periods;
 par2 = par2//par;
end;
meanpow5 = meanpow5 +meanpow4*m&odi1.[iv1,];
meandth5 = meandth5 +meandth4*m&odi1.[iv1,];
meanpc5= meanpc5 +meanpc4*m&odi1.[iv1,];
meanpe5 = meanpe5 +meanpe4*m&odi1.[iv1,];
end; /*v1*/
if &random1 = 1 then do;
 temp&odi5 = .; temp&odi4 = .; temp&odi3 = .;
 temp&odi2 = .; temp&odi1 = .;
 n_lr = int(n+0.5); meandth5 = int(meandth5+0.5);
 par = alpha||meanpow5||temp1||temp2||temp3||temp4||temp5||meanpc5||
       meanpe5||n_lr||meandth5||periods;
 par2 = par2//par;
end;
end; /*v6*/
%end;

%if &n = %then %do;
sum = &random1+&random2+&random3+&random4+&random5;
if sum = 0 then do;
/* all the parameters are fixed, not random variables */
do iv7 = 1 to nomp7; /*power*/
 power = mp7[iv7,];
do iv1 = 1 to nomp&odi1;
do iv2 = 1 to nomp&odi2;
do iv3 = 1 to nomp&odi3;
do iv4 = 1 to nomp&odi4;
do iv5 = 1 to nomp&odi5;
 if &prop = 1 then  k = mp1[iv&od1,];
 else k = mp1;
 pc = mp2[iv&od2,];
 dropin = mp3[iv&od3,];
 noncmpl = mp4[iv&od4,];
 loss = mp5[iv&od5,];
 temp1 = mp1[iv&od1,]; temp2 = iv&od2; temp3 = iv&od3; temp4 = iv&od4;
 temp5 = iv&od5;
 %markov;
%if &logr = 1
  %then %do ;
    z_beta=probit(power);
    d_lr=((z_alpha+z_beta)/ed)##2;
    n_lr=(r1+r2)#d_lr/(r1#pc2+r2#pe2);
  %end ;
  %else %do ;
    z_beta=probit(power);
    p = (r1#pc2+r2#pe2)/(r1+r2);
    sd1 = ((r1+r2)#p#(1-p)/(r1#r2))##0.5;
    sd2 = (pc2#(1-pc2)/r1+pe2#(1-pe2)/r2)##0.5;
    n_lr = ((z_alpha#sd1+z_beta#sd2)/(pc2-pe2))##2 ;
    n_lr = (r1+r2)#n_lr ;
    d_lr = n_lr#p ;
  %end ;
  n_lr = int(n_lr+0.5);
  d_lr = int(d_lr+0.5);
 /* d_lr=((z_alpha#sig+z_beta#sig)/((rho#gamma)[,+]))##2; */
 if &prop = 0 then temp1 = .;
 par = alpha||power||temp1||temp2||temp3||temp4||temp5||pc2||pe2||
      n_lr||d_lr||periods;
 par2 = par2//par;
end; /*iv5*/
end; /*iv4*/
end; /*iv3*/
end; /*iv2*/
end; /*iv1*/
end;/*iv7*/
end;

if sum ^= 0 then do;
/* use interval halving method to find n s.t.
       mean power >= specified one. */
do iv7 = 1 to nomp7;
power = mp7[iv7,];
nmax = &nmax; nmin = &nmin; tol = &ntol;
it = 0; n = (nmin+nmax)/2;
iter1:
meanpow5 = 0; meandth5 = 0; meanpc5 = 0; meanpe5 = 0;
do iv1 = 1 to nomp&odi1;
iter2:
 meanpow4 = 0; meandth4 = 0; meanpc4 = 0; meanpe4 = 0;
do iv2 = 1 to nomp&odi2;
iter3:
 meanpow3 = 0; meandth3 = 0; meanpc3 = 0; meanpe3 = 0;
do iv3 = 1 to nomp&odi3;
iter4:
 meanpow2 = 0; meandth2 = 0; meanpc2 = 0; meanpe2 = 0;
do iv4 = 1 to nomp&odi4;
iter5:
 meanpow1 = 0; meandth1 = 0; meanpc1 = 0; meanpe1 = 0;
do iv5 = 1 to nomp&odi5;
 if &prop = 1 then  k = mp1[iv&od1,];
 else k = mp1;
 pc = mp2[iv&od2,];
 dropin = mp3[iv&od3,];
 noncmpl = mp4[iv&od4,];
 loss = mp5[iv&od5,];
 temp1 = mp1[iv&od1,]; temp2=iv&od2; temp3=iv&od3; temp4=iv&od4;
 temp5 = iv&od5;
%markov;
%if &logr = 1
  %then %do ;
    d_lr = (n#r1#pc2)/(r1+r2)+(n#r2#pe2)/(r1+r2);
    z_beta = (d_lr##0.5)#ed-z_alpha;
  %end ;
  %else %do ;
    p = (r1#pc2+r2#pe2)/(r1+r2);
    sd1 = ((r1+r2)#p#(1-p)/(r1#r2))##0.5;
    sd2 = (pc2#(1-pc2)/r1+pe2#(1-pe2)/r2)##0.5;
    z_beta = (((n##0.5)#abs((pc2-pe2))/((r1+r2)##0.5)-z_alpha#sd1)/sd2) ;
    d_lr = n#p ;
  %end ;
/*  d_lr = int(d_lr+0.5);*/
newpow = probnorm(z_beta);
meanpow1 = meanpow1+newpow*m&odi5.[iv5,];
meandth1 = meandth1+d_lr*m&odi5.[iv5,];
meanpc1 = meanpc1+pc2*m&odi5.[iv5,];
meanpe1 = meanpe1+pe2*m&odi5.[iv5,];
end; /*v5*/
if &random4 = 0 then do;
 if abs(meanpow1-power)/power > tol | meanpow1 < power then do;
  it = it+1;
  if (it > itmax) then do;
   print " Number of iterations has exceeded the maximum number.";
   stop;
  end;
  else do;
  if meanpow1 > power then nmax = n;
  else nmin = n;
  n = (nmin+nmax)/2;
  goto iter5;
  end;
 end;
 n_lr = n; nmin=&nmin; nmax = &nmax;
 temp&odi5 = .; n_lr = int(n_lr+0.5); meandth1 = int(meandth1+0.5);
 par = alpha||meanpow1||temp1||temp2||temp3||temp4||temp5||pc2||pe2||
     n_lr||meandth1||periods;
 par2 = par2//par;
 it = 0;
end;
meanpow2 = meanpow2+meanpow1*m&odi4.[iv4,];
meandth2 = meandth2+meandth1*m&odi4.[iv4,];
meanpc2 = meanpc2+meanpc1*m&odi4.[iv4,];
meanpe2 = meanpe2+meanpe1*m&odi4.[iv4,];
end; /*v4*/
if &random3 = 0 & &random4 = 1 then do;
 if abs(meanpow2-power)/power > tol | meanpow2 < power then do;
  it = it+1;
  if (it > itmax) then do;
   print " Number of iterations has exceeded the maximum number.";
   stop;
  end;
  if meanpow2 > power then nmax=n;
  else nmin = n;
  n = (nmin+nmax)/2;
  goto iter4;
 end;
 n_lr = n; nmin=&nmin;nmax=&nmax; temp&odi5 = .; temp&odi4 = .;
 n_lr = int(n_lr+0.5); meandth2 = int(meandth2+0.5);
 par = alpha||meanpow2||temp1||temp2||temp3||temp4||temp5||pc2||pe2||
      n_lr||meandth2||periods;
 par2 = par2//par;
 it = 0;
end;
meanpow3 = meanpow3+meanpow2*m&odi3.[iv3,];
meandth3 = meandth3+meandth2*m&odi3.[iv3,];
meanpc3 = meanpc3+meanpc2*m&odi3.[iv3,];
meanpe3 = meanpe3+meanpe2*m&odi3.[iv3,];
end; /*v3*/
if &random2 = 0 & &random3 = 1 then do;
 if abs(meanpow3-power)/power > tol | meanpow3 < power then do;
  it = it+1;
  if (it > itmax) then do;
   print " Number of iterations has exceeded the maximum number.";
   stop;
  end;
  if meanpow3 > power then nmax=n;
  else nmin = n;
  n = (nmin+nmax)/2;
  goto iter3;
 end;
 n_lr = n; nmin=&nmin; nmax=&nmax;temp&odi5 = .; temp&odi4=.; temp&odi3=.;
 n_lr = int(n_lr+0.5); meandth3 = int(meandth3+0.5);
 par = alpha||meanpow3||temp1||temp2||temp3||temp4||temp5||pc2||pe2||
      n_lr||meandth3||periods;
 par2 = par2//par;
 it = 0;
end;
meanpow4 = meanpow4+meanpow3*m&odi2.[iv2,];
meandth4 = meandth4+meandth3*m&odi2.[iv2,];
meanpc4 = meanpc4+meanpc3*m&odi2.[iv2,];
meanpe4 = meanpe4+meanpe3*m&odi2.[iv2,];
end; /*v2*/
if &random1 = 0 & &random2 = 1 then do;
 if abs(meanpow4-power)/power > tol | meanpow4 < power then do;
  it = it+1;
  if (it > itmax) then do;
   print " Number of iterations has exceeded the maximum number.";
   stop;
  end;
  if meanpow4 > power then nmax=n;
  else nmin = n;
  n = (nmin+nmax)/2;
  goto iter2;
 end;
 n_lr = n; nmin=&nmin; nmax=&nmax;
 temp&odi5 = .; temp&odi4=.; temp&odi3=.; temp&odi2=.;
 n_lr = int(n_lr+0.5); meandth4 = int(meandth4+0.5);
 par = alpha||meanpow4||temp1||temp2||temp3||temp4||temp5||pc2||pe2||
      n_lr||meandth4||periods;
 par2 = par2//par;
 it = 0;
end;
meanpow5 = meanpow5+meanpow4*m&odi1.[iv1,];
meandth5 = meandth5+meandth4*m&odi1.[iv1,];
meanpc5 = meanpc5+meanpc4*m&odi1.[iv1,];
meanpe5 = meanpe5+meanpe4*m&odi1.[iv1,];
end; /*v1*/
if &random1 = 1 then do;
 if abs(meanpow5-power)/power > tol | meanpow5 < power then do;
  it = it+1;
  if (it > itmax) then do;
   print " Number of iterations has exceeded the maximum number.";
   stop;
  end;
  if meanpow5 > power then nmax = n;
  else nmin = n;
  n = (nmin+nmax)/2;
  goto iter1;
 end;
 n_lr = n; nmin=&nmin; nmax=&nmax; temp&odi5 = .; temp&odi4=.; temp&odi3=.;
 temp&odi2=.; temp&odi1=.; n_lr = int(n_lr+0.5); meandth5 = int(meandth5+0.5);
 par = alpha||meanpow5||temp1||temp2||temp3||temp4||temp5||pc2||pe2||
      n_lr||meandth5||periods;
 par2 = par2//par;
 it = 0;
end;
end; /*power*/
end;
%end;

%if &power ^= & &n ^= %then %do;
if recrate = 0 then do;
 print " Recruitment rate was not specified.";
 stop;
end;


/* Use interval halving method to find the duration time s.t.
mean power >= specified one. */
do iv7 = 1 to nomp7; /*power*/
do iv6 = 1 to nomp6; /* n */
power = mp7[iv7,];
n = mp6[iv6,];
npmin = &npmin; npmax = &npmax; tol = &nptol;
it = 0;
iter12:
meanpow5 = 0; meandth5 = 0; meanpc5 = 0; meanpe5 = 0;
do iv1 = 1 to nomp&odi1;
iter22:
 meanpow4 = 0; meandth4 = 0; meanpc4 = 0; meanpe4 = 0;
do iv2 = 1 to nomp&odi2;
iter32:
 meanpow3 = 0; meandth3 = 0; meanpc3 = 0; meanpe3 = 0;
do iv3 = 1 to nomp&odi3;
iter42:
 meanpow2 = 0; meandth2 = 0; meanpc2 = 0; meanpe2 = 0;
do iv4 = 1 to nomp&odi4;
iter52:
 meanpow1 = 0; meandth1 = 0; meanpc1 = 0; meanpe1 = 0;
do iv5 = 1 to nomp&odi5;
iter62:
 free con;
 nptmp = (npmin+npmax)/2; nptmp2 = ceil(nptmp);
 if nptmp2 > periods then do;
 con= j(1,(nptmp2-periods),mp2[1,periods]); /* pc */
 if nomp2 > 1 then do;
  do k = 2 to nomp2;
  con2= j(1,(nptmp2-periods),mp2[k,periods]); /* pc */
  con = con//con2;
  end;
 end;
 mp2 =mp2[1:nomp2,1:periods]||con; /*exp& the vect| of pc rates */
 end;

 if nptmp2 > periods then do;
 con= j(1,(nptmp2-periods),mp3[1,periods]); /* dropin */
 if nomp3 > 1 then do;
  do k = 2 to nomp3;
  con2= j(1,(nptmp2-periods),mp3[k,periods]); /* dropin */
  con = con//con2;
  end;
 end;
 mp3 =mp3[1:nomp3,1:periods]||con;
end;

 if nptmp2 > periods then do;
 con= j(1,(nptmp2-periods),mp4[1,periods]); /* dropout */
 if nomp4 > 1 then do;
  do k = 2 to nomp4;
  con2= j(1,(nptmp2-periods),mp4[k,periods]); /* dropout */
  con = con//con2;
  end;
 end;
  mp4 =mp4[1:nomp4,1:periods]||con;
 end;
 if nptmp2 > periods then do;
 con= j(1,(nptmp2-periods),mp5[1,periods]); /* loss of follow-up */
 if nomp5 > 1 then do;
  do k = 2 to nomp5;
  con2= j(1,(nptmp2-periods),mp5[k,periods]); /* loss of follow-up */
  con = con//con2;
  end;
 end;
  mp5 =mp5[1:nomp5,1:periods]||con;
 end;

 /* Adjust for staggered entry. */
n_intrvl=sbdv;
 nn = ceil(nptmp#n_intrvl); /* no. of subperiods */
 ad_cens=j(1,nn,0); rcrt_sum=0;
 nprec = n/recrate;
 minrec = nprec; if (nptmp < nprec) then minrec = nptmp;
 ntmp = minrec#recrate;
 nprec = minrec#n_intrvl;
 rcrt = j(1,nn,0);
 rcrt[,1:nprec] = j(1,nprec,recrate);
 do ii = 1 to nn;
  rcrt_sum = rcrt_sum+rcrt[,ii];
  ad_cens[,nn+1-ii] = rcrt[,ii]/rcrt_sum;
 end;
 if &prop = 1 then do;
 k = mp1[iv&od1,];
 end;
 else do;
 k = mp1;
 end;
 pc = mp2[iv&od2,];
 dropin = mp3[iv&od3,];
 noncmpl = mp4[iv&od4,];
 loss = mp5[iv&od5,];
 temp1 = mp1[iv&od1,]; temp2=iv&od2; temp3=iv&od3; temp4=iv&od4;
  temp5 = iv&od5;
%markov2;
%if &logr = 1
  %then %do ;
    d_lr = (ntmp#r1#pc2)/(r1+r2)+(ntmp#r2#pe2)/(r1+r2);
    z_beta = (d_lr##0.5)#ed-z_alpha;
  %end ;
  %else %do ;
    p = (r1#pc2+r2#pe2)/(r1+r2);
    sd1 = ((r1+r2)#p#(1-p)/(r1#r2))##0.5;
    sd2 = (pc2#(1-pc2)/r1+pe2#(1-pe2)/r2)##0.5;
    z_beta = (((ntmp##0.5)#abs((pc2-pe2))/((r1+r2)##0.5)-z_alpha#sd1)/sd2) ;
    d_lr = n#p ;
  %end ;
newpow = probnorm(z_beta);
if &prop = 0 then temp1= .;
if &random5 = 0 then do; /* not random */
  if abs(newpow-power)/power > tol | newpow < power then do;
   it = it+1;
   if (it > itmax) then do;
    print " no. of interations has exceeded the maximum number.";
    stop;
   end;
   if newpow > power then
    npmax = nptmp;
   else
    npmin = nptmp;
   goto iter62;
  end;
  n_lr = ceil(ntmp); d_lr = int(d_lr+0.5);
  par = alpha||newpow||temp1||iv&od2||iv&od3||iv&od4||iv&od5||pc2||pe2||
        n_lr||d_lr||nptmp;
  par2 = par2//par;
  it = 0;
  npmin = &npmin; npmax = &npmax;
 end;
 meanpow1 = meanpow1+newpow*m&odi5.[iv5,];
 meandth1 = meandth1+d_lr*m&odi5.[iv5,];
 meanpc1 = meanpc1+pc2*m&odi5.[iv5,];
 meanpe1 = meanpe1+pe2*m&odi5.[iv5,];
end; /*v5*/
 if &random4 = 0 & &random5 = 1 then do;
  if abs(meanpow1-power)/power | meanpow1 < power > tol then do;
   it = it+1;
   if (it > itmax) then do;
    print " no. of iterations has exceeded the maximum number.";
    stop;
   end;
   if meanpow1 > power then
    npmax = nptmp;
   else
    npmin = nptmp;
   goto iter52;
  end;
  n_lr = ceil(ntmp); npmin=&npmin; npmax = &npmax;
  temp&odi5 = .; meandth1 = int(meandth1+0.5);
  par = alpha||meanpow1||temp1||temp2||temp3||temp4||temp5||meanpc1||
        meanpe1||n_lr||meandth1||nptmp;
  par2 = par2//par;
  it = 0;
 end;
 meanpow2 = meanpow2+meanpow1*m&odi4.[iv4,];
 meandth2 = meandth2+meandth1*m&odi4.[iv4,];
 meanpc2 = meanpc2+meanpc1*m&odi4.[iv4,];
 meanpe2 = meanpe2+meanpe1*m&odi4.[iv4,];
end; /*v4*/
 if &random3 = 0 & &random4 = 1 then do;
  if abs(meanpow2-power)/power > tol | meanpow2 < power then do;
   it = it+1;
   if (it > itmax) then do;
    print " no. of iterations has exceeded the maximum number.";
    stop;
   end;
   if meanpow2 > power then
    npmax = nptmp;
   else
    npmin = nptmp;
   goto iter42;
  end;
  n_lr = ntmp; npmin=&npmin;npmax=&npmax;
  temp&odi5 = .; temp&odi4 = .;
  meandth2 = int(meandth2+0.5);
  if &prop = 0 then temp1 = .;
  par = alpha||meanpow2||temp1||temp2||temp3||temp4||temp5||meanpc2||
        meanpe2||n_lr||meandth2||nptmp;
  par2 = par2//par;
  it = 0;
 end;
 meanpow3 = meanpow3+meanpow2*m&odi3.[iv3,];
 meandth3 = meandth3+meandth2*m&odi3.[iv3,];
 meanpc3 = meanpc3+meanpc2*m&odi3.[iv3,];
 meanpe3 = meanpe3+meanpe2*m&odi3.[iv3,];
end; /*v3*/
 if &random2 = 0 & &random3 = 1 then do;
  if abs(meanpow3-power)/power > tol | meanpow3 < power then do;
   it = it+1;
   if (it > itmax) then do;
    print " no. of iterations has exceeded the maximum number.";
    stop;
   end;
   if meanpow3 > power then
    npmax = nptmp;
   else
    npmin = nptmp;
   goto iter32;
  end;
  n_lr = ntmp; npmin=&npmin; npmax=&npmax;
  temp&odi5 = .; temp&odi4=.; temp&odi3=.;
  meandth3 = int(meandth3+0.5);
  if &prop = 0 then temp1 = .;
  par = alpha||meanpow3||temp1||temp2||temp3||temp4||temp5||meanpc3||
        meanpe3||n_lr||meandth3||nptmp;
  par2 = par2//par;
  it = 0;
 end;
 meanpow4 = meanpow4+meanpow3*m&odi2.[iv2,];
 meandth4 = meandth4+meandth3*m&odi2.[iv2,];
 meanpc4 = meanpc4+meanpc3*m&odi2.[iv2,];
 meanpe4 = meanpe4+meanpe3*m&odi2.[iv2,];
end; /*v2*/
 if &random1 = 0 & &random2 = 1 then do;
  if abs(meanpow4-power)/power > tol | meanpow4 < power then do;
   it = it+1;
   if (it > itmax) then do;
    print " no. of interations has exceeded the maximum number.";
    stop;
   end;
   if meanpow4 > power then
    npmax = nptmp;
   else
    npmin = nptmp;
   goto iter62;
  end;
  n_lr = ntmp; npmin=&npmin; npmax=&npmax;
  temp&odi5 = .; temp&odi4=.; temp&odi3=.;
  temp&odi2=.;
  meandth4 = int(meandth4+0.5);
  if &prop = 0 then temp1 = .;
  par = alpha||meanpow4||temp1||temp2||temp3||temp4||temp5||meanpc4||
        meanpe4||n_lr||meandth4||nptmp;
  par2 = par2//par;
  it = 0;
 end;
 meanpow5 = meanpow5+meanpow4*m&odi1.[iv1,];
 meandth5 = meandth5+meandth4*m&odi1.[iv1,];
 meanpc5 = meanpc5+meanpc4*m&odi1.[iv1,];
 meanpe5 = meanpe5+meanpe4*m&odi1.[iv1,];
 end; /*v1*/
 if &random1 = 1 then do;
 if abs(meanpow5-power)/power > tol | meanpow5 < power then do;
  it = it+1;
   if (it > itmax) then do;
    print " no. of interations has exceeded the maximum number.";
    stop;
   end;
   if manpow5 > power then
    npmax = nptmp;
   else
    npmin = nptmp;
   goto iter62;
  end;
  n_lr = ntmp; npmin=&npmin; npmax=&npmax;
  temp&odi5 = .; temp&odi4=.; temp&odi3=.;
  temp&odi2=.;
  temp&odi1=.;
  meandth5 = int(meandth5+0.5);
  if &prop = 0 then temp1 = .;
  par = alpha||meanpow5||temp1||temp2||temp3||temp4||temp5||meanpc5||
        meanpe5||n_lr||meandth5||nptmp;
  par2 = par2//par;
  it = 0;
 end;
 end; /* n */
 end; /* power */
%end;
%mend lakprog;

%macro write;
name = {"alpha" "power" "k" "pcno" "din" "dout" "loss" "pc" "pe"
        "n" "noevent" "duration"};
nr = nrow(par2);
nc = ncol(par2);
rname = j(1,nr,"");
/*par2[,2] = round(par[,2],0.001);
par2[,8] = round(par[,8],0.001);
par2[,9] = round(par[,9],0.001);
*/

print par2[rowname=rname colname=name format=8.3];
%if &output ne %then %do; filename out "&output"; %end;
%else %do; filename out 'out.list'; %end;
/*filename out 'out.list';*/
file out;
 do i = 1 to nr;
  do j = 1 to nc;
   put (par2[i,j]) 8.3 +1 @;
  end;
  put /;
 end;
closefile out;
%mend write;

%macro markov;
if (&prop = 1 | (&prop = 0 & &lag = 0))then do;
distr_e={0,0,1,0}; distr_c={0,0,0,1};
free dstr_e dstr_c;
trans = i(4);
do year=1 to periods;
  ls=1-(1-loss[,year])##(1/n_intrvl);
  dro=1-(1-noncmpl[,year])##(1/n_intrvl);
  dri=1-(1-dropin[,year])##(1/n_intrvl);
  pc1=1-(1-pc[,year])##(1/n_intrvl);
  if &prop = 1 then
  pe1 = 1-exp(k#log(1-pc1));
  else
  pe1 = 1-exp(k[year,]#log(1-pc1));
  do ii=1 to n_intrvl;
    trans[,3]=ls//pe1//(1-(ls+pe1+dro))//dro;
    trans[,4]=ls//pc1//dri//(1-(ls+pc1+dri));
    distr_e=trans*distr_e; distr_c=trans*distr_c;
*the next 6 lines adjust for staggered entry;
    temp_e=distr_e[{3 4},1]#(1-ad_cens[,ii+(year-1)#n_intrvl]);
    distr_e[1,]=distr_e[1,]+(distr_e[{3 4},]-temp_e)[+,];
    distr_e[{3 4},]=temp_e;
    temp_c=distr_c[{3 4},1]#(1-ad_cens[,ii+(year-1)#n_intrvl]);
    distr_c[1,]=distr_c[1,]+(distr_c[{3 4},]-temp_c)[+,];
    distr_c[{3 4},]=temp_c;
  *end of transition matrix loop;   end;
    dstr_e=dstr_e||distr_e;  dstr_c=dstr_c||distr_c;
  *end of years loop;   end;
*end if clause for prop hazards or nonlag models; end;
else do; /* nonproportional hazards and lag is specified */
nactv = lag#n_intrvl;
nstates=2#nactv+2;
distr_e = 0//0//1//j(nactv-1,1,0)//j(nactv,1,0);
distr_c = 0//0//j(nactv,1,0)//1//j(nactv-1,1,0);
free dstr_e dstr_c;

do year=1 to periods;
  ls=1-(1-loss[,year])##(1/n_intrvl)#j(1,nactv,1);
  dro=1-(1-noncmpl[,year])##(1/n_intrvl)#j(1,nactv,1);
  dri=1-(1-dropin[,year])##(1/n_intrvl)#j(1,nactv,1);
  pc1=1-(1-pc[,year])##(1/n_intrvl);
  pe1 = j(1,nactv,0);
  do l = 1 to lag;
    pe1[,(l-1)#n_intrvl+1:l#n_intrvl] = 1-(exp(k[l,]#log(1-pc[,year])#
     j(1,n_intrvl,1)))##(1/n_intrvl);
  end;
  actv = pc1||pe1;
  a = j(nactv,nactv,0); b=a; d=a;
  *Start transition creationn and multiplication loop;
  c = diag(dro); b= diag(dri);
  a[2:nactv,1:(nactv-1)] =
   diag(1-((ls+dro+actv[,2:(nactv+1)])[,1:(nactv-1)]));
  d[1:(nactv-1),2:nactv] =
   diag(1-(((ls+dri)[,1:(nactv-1)]+actv[,2:nactv])));
  a[nactv,nactv] = 1-(ls[,1]+dro[,1]+actv[,nactv+1]);
  d[1,1] = 1-((ls+dri)[,1]+actv[,1]);
  do ii=1 to n_intrvl;
    trans=(i(2)||((ls||ls)//(actv[,2:nactv+1]||actv[,1:(nactv)])))
          //(j(2#nactv,2,0)||((a||b)//(c||d)));
    distr_e=trans*distr_e; distr_c=trans*distr_c;
*the next 6 lines adjust for staggered entry;
    temp_e=distr_e[3:nstates,1]#(1-ad_cens[,ii+(year-1)#n_intrvl]);
    distr_e[1,]=distr_e[1,]+(distr_e[3:nstates,]-temp_e)[+,];
    distr_e[3:nstates,]=temp_e;
    temp_c=distr_c[3:nstates,1]#(1-ad_cens[,ii+(year-1)#n_intrvl]);
    distr_c[1,]=distr_c[1,]+(distr_c[3:nstates,]-temp_c)[+,];
    distr_c[3:nstates,]=temp_c;
  *end of transition matrix loop;   end;
    dstr_e=dstr_e||distr_e;  dstr_c=dstr_c||distr_c;
  *end of years loop;   end;
  collaps1=3:(nactv+2); collaps2=(nactv+3):(2#nactv+2);
  dstr_e=dstr_e[{1 2},]//dstr_e[collaps1,][+,]//dstr_e[collaps2,][+,];
  dstr_c=dstr_c[{1 2},]//dstr_c[collaps1,][+,]//dstr_c[collaps2,][+,];

end; /*end nonproportional lag model */

if periods = 1 then do; event_c=dstr_c[2,];

event_e=dstr_e[2,];
loss_c=dstr_c[1,];
loss_e=dstr_e[1,];
end;
else do;
event_c=dstr_c[2,]-({0}||dstr_c[2,1:(ncol(dstr_c)-1)]);
event_e=dstr_e[2,]-({0}||dstr_e[2,1:(ncol(dstr_e)-1)]);
loss_c=dstr_c[1,]-({0}||dstr_c[1,1:(ncol(dstr_c)-1)]);
loss_e=dstr_e[1,]-({0}||dstr_e[1,1:(ncol(dstr_e)-1)]);
end;
atrisk_c=(dstr_c[{3 4},][+,]+loss_c+event_c);
atrisk_e=(dstr_e[{3 4},][+,]+loss_e+event_e);
phi=atrisk_c#r1/(atrisk_e#r2);
theta=log(1-event_c/atrisk_c)/log(1-event_e/atrisk_e);
/*theta=(event_c/atrisk_c)/(event_e/atrisk_e); */
rho=(event_c#r1+event_e#r2)/((event_c#r1+event_e#r2)[,+]);
gamma=abs(phi#theta/(1+phi#theta)-phi/(1+phi));
eta=phi/((1+phi)##2);
sig=sqrt((rho#eta)[,+]);
pe2=distr_e[2,]; pc2=distr_c[2,];  pbar=(r1#pc2+r2#pe2)/(r1+r2);
z_alpha=probit(1-alpha/2);
sig2 = (rho#gamma)[,+];
ed = sig2/sig;
%mend markov;

%macro markov2;
free dstr_e dstr_c;
if (&prop = 1 | (&prop = 0 & &lag = 0))then do;
distr_e={0,0,1,0}; distr_c={0,0,0,1};
trans = i(4);
do year = 1 to nptmp2;
 ls=1-(1-loss[,year])##(1/n_intrvl);
 dro=1-(1-noncmpl[,year])##(1/n_intrvl);
 dri=1-(1-dropin[,year])##(1/n_intrvl);
 pc1=1-(1-pc[,year])##(1/n_intrvl);
 if &prop = 1 then
 pe1 = 1-exp(k#log(1-pc1));
 else
 pe1 = 1-exp(k[year,]#log(1-pc1));
 if (year = nptmp2 & nptmp2 ^= nptmp) then
  ni = mod(nn,n_intrvl);
 else  ni = n_intrvl;
 do ii = 1 to ni;
 trans[,3]=ls//pe1//(1-(ls+pe1+dro))//dro;
 trans[,4]=ls//pc1//dri//(1-(ls+pc1+dri));
 distr_e=trans*distr_e; distr_c=trans*distr_c;
*the next 6 lines adjust for staggered entry;
 temp_e=distr_e[{3 4},1]#(1-ad_cens[,ii+(year-1)#ni]);
 distr_e[1,]=distr_e[1,]+(distr_e[{3 4},]-temp_e)[+,];
 distr_e[{3 4},]=temp_e;
 temp_c=distr_c[{3 4},1]#(1-ad_cens[,ii+(year-1)#ni]);
 distr_c[1,]=distr_c[1,]+(distr_c[{3 4},]-temp_c)[+,];
 distr_c[{3 4},]=temp_c;
*end of transition matrix loop;   end;
 dstr_e=dstr_e||distr_e;  dstr_c=dstr_c||distr_c;
*end of years loop;   end;
end;

else do; /* nonproportional hazards and lag is specified */
nactv = lag#n_intrvl;
nstates=2#nactv+2;
distr_e = 0//0//1/j(nactv-1,1,0)//j(nactv,1,0);
distr_c = 0//0//j(nactv,1,0)//1//j(nactv-1,1,0);
do year=1 to nptmp2;
  ls=1-(1-loss[,year])##(1/n_intrvl)#j(1,nactv,1);
  dro=1-(1-noncmpl[,year])##(1/n_intrvl)#j(1,nactv,1);
  dri=1-(1-dropin[,year])##(1/n_intrvl)#j(1,nactv,1);
  pc1=1-(1-pc[,year])##(1/n_intrvl)#j(1,n_intrvl,1);
  pe1 = j(1,nactv,0);
  do l = 1 to lag;
    pe1[,(l-1)#n_intrvl+1:l#n_intrvl] = 1-(exp(k[l,]#log(1-pc[,year])#
     j(1,n_intrvl,1)))##(1/n_intrvl);
  end;
  actv = pc1[,1]||pe1;
  a = j(nactv,nactv,0); b=a; d=a;
  *Start transition creationn and multiplication loop;
  c = diag(dro); b= diag(dri);
  a[2:nactv,1:(nactv-1)] =
   diag(1-((ls+dro+actv[,2:(nactv+1)])[,1:(nactv-1)]));
  d[1:(nactv-1),2:nactv] =
   diag(1-(((ls+dri)[,1:(nactv-1)]+actv[,2:nactv])));
  a[nactv,nactv] = 1-(ls[,1]+dro[,1]+actv[,nactv+1]);
  d[1,1] = 1-((ls+dri)[,1]+actv[,1]);
  if (year = nptmp2 & nptmp2 ^= nptmp) then
   ni = mod(nn,n_intrvl);
  else  ni = n_intrvl;
  do ii=1 to ni;
    trans=(i(2)||((ls||ls)//(actv[,2:nactv+1]||actv[,1:(nactv)])))
          //(j(2#nactv,2,0)||((a||b)//(c||d)));
    distr_e=trans*distr_e; distr_c=trans*distr_c;
*the next 6 lines adjust for staggered entry;
    temp_e=distr_e[3:nstates,1]#(1-ad_cens[,ii+(year-1)#n_intrvl]);
    distr_e[1,]=distr_e[1,]+(distr_e[3:nstates,]-temp_e)[+,];
    distr_e[3:nstates,]=temp_e;
    temp_c=distr_c[3:nstates,1]#(1-ad_cens[,ii+(year-1)#n_intrvl]);
    distr_c[1,]=distr_c[1,]+(distr_c[3:nstates,]-temp_c)[+,];
    distr_c[3:nstates,]=temp_c;
  *end of transition matrix loop;   end;
    dstr_e=dstr_e||distr_e;  dstr_c=dstr_c||distr_c;
  *end of years loop;   end;
  collaps1=3:(nactv+2); collaps2=(nactv+3):(2#nactv+2);
  dstr_e=dstr_e[{1 2},]//dstr_e[collaps1,][+,]//dstr_e[collaps2,][+,];
  dstr_c=dstr_c[{1 2},]//dstr_c[collaps1,][+,]//dstr_c[collaps2,][+,];
end; /*end nonproportional lag model */

if nptmp2 = 1 then do;
event_c=dstr_c[2,];
event_e=dstr_e[2,];
loss_c=dstr_c[1,]-({0}||dstr_c[1,1:(ncol(dstr_c)-1)]);
loss_e=dstr_e[1,]-({0}||dstr_e[1,1:(ncol(dstr_e)-1)]);
end;
else do;
event_c=dstr_c[2,]-({0}||dstr_c[2,1:(ncol(dstr_c)-1)]);
event_e=dstr_e[2,]-({0}||dstr_e[2,1:(ncol(dstr_e)-1)]);
loss_c=dstr_c[1,]-({0}||dstr_c[1,1:(ncol(dstr_c)-1)]);
loss_e=dstr_e[1,]-({0}||dstr_e[1,1:(ncol(dstr_e)-1)]);
end;
atrisk_c=(dstr_c[{3 4},][+,]+loss_c+event_c);
atrisk_e=(dstr_e[{3 4},][+,]+loss_e+event_e);
phi=atrisk_c#r1/(atrisk_e#r2);
/*theta=(event_c/atrisk_c)/(event_e/atrisk_e);*/
theta=log(1-event_c/atrisk_c)/log(1-event_e/atrisk_e);
rho=(event_c#r1+event_e#r2)/((event_c#r1+event_e#r2)[,+]);
gamma=abs(phi#theta/(1+phi#theta)-phi/(1+phi));
eta=phi/((1+phi)##2);
sig=sqrt((rho#eta)[,+]);
pe2=distr_e[2,]; pc2=distr_c[2,];  pbar=(r1#pc2+r2#pe2)/(r1+r2);
z_alpha=probit(1-alpha/2);
sig2 = (rho#gamma)[,+];
ed = sig2/sig;

%mend markov2;