/*iml_LEAAF_for_Bootstrap.sas*/

/*Step 1-2: Read SAS dataset in the following order:
Put dichotomous risk factors that you want to calculate the
LE-AAFs for First. In this example it is arthritis to highschool. Next
put the dichotomous covariates (in this case age_dichot and female). 
Next put indicator variables for each time point. In this example
there are 9 time periods, so the indicators are time1 to time9.
Finally, the outcome (hospitalization in our case), and then a
continuous version of the time variable*/
/*User input is highlighted in Yellow*/

proc iml;
varnames={arthritis asthma afib kidney copd dementia depress diabetes hf 
hyperlip htn ischemic ami osteo stroke_tia pro_medicaid_status highschool
age_dichot female time1 time2 time3 time4 time5 time6 time7 time8 time9
hosp_outcome time};
use work.prepare_&i;
read all var _num_ into originaldata[colname=varnames];
close work.prepare_&i;

/*Step 3: Read the design dataset into matrix desmtx*/
varnames2={arthritis asthma afib kidney copd dementia depress diabetes hf 
hyperlip htn ischemic ami osteo stroke_tia pro_medicaid_status highschool age_dichot female};
use bootdoc.design_whites_&i;
read all var _num_ into desmtx[colname=varnames2];
close bootdoc.design_whites_&i;

/* Step 4: Enter the number of risk factors you have read into matrices.
   Risk factors only not covariates*/
nrisk=17; 

/* Step 5: Enter the number of covariates you have read into matrices.*/
ncov=2;

/* Step 6: Enter the maximum number of time points.*/
maxyear=9;


nx=nrisk+ncov; 
ninter=0; 
nmandc=nrisk;
x=originaldata[, (1:(nx+maxyear+1))] ;  
*x=originaldata[, (1:(nx+2))] ;
xsub=originaldata[, (1:(nx))] ;
year=originaldata[, ((nx+1):(nx+maxyear))];
*year=originaldata[, (1:(nx+2))];

outcomea=originaldata[,(nx+maxyear+1)];

N=nrow(x);

package load listutil;

lendesmtx=nrow(desmtx);
/*Select the columns for the risk factors 
  from the design matrix*/
sumnonzero1=desmtx[,1:Nrisk];
*Sum the cols of the design matrix risk factors only;
sumnonzeros=sumnonzero1[,+] ;
*Get maximum # of risk factors; 
maxnonzeros=max(sumnonzeros);

/*Function that creates factorial matrix of 2 ## n risk factors*/
start FullFactorial( n );
  r = 2 ## n;
  xf = j(r, n);
  pot = 2 ## ( (n-1):0 );
  do i = 1 to r;
  xf[i, ] = band( i-1, pot ) > 0;
end;
  return( xf );
finish;

nz2enum=listcreate(maxnonzeros);
subsets=listcreate(maxnonzeros);
loczero=listcreate(maxnonzeros);

do i=1 to maxnonzeros;
valf=0:(2##i-1);
call ListSetItem (nz2enum,i,valf);
call ListSetItem (subsets,i,{0});
call ListSetItem (loczero,i,{0});
a = FullFactorial( i );
call ListSetItem (subsets,i,a);
call ListSetItem (loczero, i, loc(a)); /*Locations of ones*/ 
end;

/*Step 7 Read SAS dataset of vector of logistic regression coefficients into 
SAS Order should be intercept, time, risk factors, covariates.
Should only be one variable in dataset*/

varnames3={betacoeff};
use bootdoc.outbeta_hosp_whites_&i;
read all var _num_ into outbeta_condIntercova[colname=varnames3];
close bootdoc.outbeta_hosp_whites_&i;

betayeara=outbeta_condIntercova[1];
betayrona=outbeta_condIntercova[2];
rowbetaa=nrow(outbeta_condIntercova);
betaa=outbeta_condIntercova[(3):rowbetaa];
betaconda=betaa[1:nrisk];
betaintera=betaa[(nrisk+1):(nrisk+ninter)];
if NCov ^= 0 then do;
betacova=outbeta_condIntercova[(2+nrisk+ninter+1):rowbetaa];
end;

AAFbyYearMtxa=j(maxyear,nmandc,0);

NskByYearMtx=j(maxyear,lendesmtx,0);

N1a=j(maxyear,1,0);

Nyr=j(maxyear,1,0);

do yr=1 to maxyear;

N1a[yr]=sum( outcomea[loc(year[,yr]=1),] );

beta0a=betayeara[1]+(betayrona[1]*yr);

xmat=xsub[loc(year[,yr]=1),];

Nyr[yr]=nrow(xmat);

xPs0va = j(NMandC,1,0);

sumbetacova=j(1,1,0);

/*Below looking for row of xmat that matches row i of
desgnmatrix. If sum of logical elements are equal then it is a match*/

do r=1 to lendesmtx; 
   xmat2=repeat(desmtx[r,],Nyr[yr],1);
   xmat3=(xmat[,]=xmat2[,]);
   temv=xmat3[,+];
   temv2=(temv[,]=nx);
   NskByYearMtx[yr,r]=temv2[+,];

   row=DesMtx[r,];
   rowMandC=row[,1:NMandC];
   numNonZeros=row[,+];

*C Loop;
if Ncov > 0  then do; 
rowcov=row[,(1+NMandC):NX];
sumBetaCova=rowcov*betacova; 
end; /*C Loop*/

rowcond=row[,1:NMandC];
If Ncov=0 then sumbetaCova=0;

/*First get indicators for columns*/
/*Find indices for nonzero elements in design matrix row, for meds and conditions*/
nzsloc=loc(rowMandC);  
numNzs=ncol(nzsloc);   /*Count how many not zero*/

*D Loop;
if numNzs > 0 then do;

subs = FullFactorial( numNzs );
numSubs=nrow(subs);
ps0vectora=j(numsubs,1,0);
total=j(numsubs,1,0);

*E Loop;
do j=1 to numsubs;
rowsubMandC=j(1,NMandC,0);
mat_row=j(j,numNzs,0);

*F Loop;
do k=1 to numnzs;
mat_row[j,k]=subs[j,k]*nzsloc[1,k];
scale=mat_row[j,k];
if scale ^=0 then rowsubMandC[1,scale]=1;
end;/*F Loop*/

rowsubcond=rowsubMandC[1,1:NMandC];
condsuma=rowsubcond*betaconda;

sumcoeffa=((-1*beta0a) - condsuma -sumBetaCova);

denoma=1+exp(sumcoeffa);

ps0vectora[j]=1 / denoma;

end; /*E loop*/

/*DO Loop G*/
*********************************************;
*Need to generate a nz2enum list with 3 rows.;
*These are cells with 0,1  0,1,2,3 0,1,2,3,4,5,6,7;
************************************************; 
enum=(listgetitem(nz2enum,numnzs))`;
do g=1 to numnzs;
xk=nzsloc[g];
p=2##(numnzs-g);
x = (band(enum,p)=p);
y1a=sum(Ps0vectora[loc(x[,1]=1),]);

y2a=sum(Ps0vectora[loc(x[,1]=0),]);

combinea=(y1a-y2a)/(numsubs/2);

combinea=combinea/n1a[yr];

combinea=NskByYearMtx[yr,r]*combinea;

xPs0va[xk]=xPs0va[xk]+combinea;

end; /*G Loop*/

end; /*D Loop*/
end; /*A Loop*/
AAFbyYearMtxa[yr,]=xPs0va`;

AAFbyYearMtx2a=AAFbyYearMtxa#Nyr;

end;

allaveragea=(AAFbyYearMtx2a[+,])/N;

sumLEAAFa=sum(allaveragea);

create outa from allaveragea;  /* create data set and variables */
append from allaveragea;              /* write all rows */
close;

/*Step 8-Specify SAS dataset name to output LE-AAFs*/
create LEAAF_all_white_&i var {allaveragea sumLEAAFa};
append;
close LEAAF_all_white_&i;
run;
quit;
/*Step 9-Specify SAS dataset in the set statement below*/
data bootdoc.LEAAF_all_white_&i;
set LEAAF_all_white_&i;
LE_AAFA=ALLAVERAGEA*100;
LE_AAFA=round(LE_AAFA,.01);
run;
proc datasets nolist noprint;
delete outa LEAAF_all_white_&i;
run;
quit;