/*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;