%macro gamfrail(mydata,factors,initbeta,convcrit,leftend,stepsize,options,data1,data2,data3);

use &mydata;
read all var _num_ into x;

covlist=&factors;
frailty={'Frailty'};
fcovlist=frailty||covlist;
itend=100;
max=0;
mid=0;
left=0;
p=ncol(x)-3;
d=0;
countall=1;
countmax=0;
countmin=0; 
groupnum=x[1,3];
summary=j(1,4,0);
do i=1 to nrow(x);
   if x[i,3]>groupnum then groupnum=x[i,3];
   if x[i,2]=1 then summary[1,2]=summary[1,2]+1;
   else if x[i,2]=0 then summary[1,3]=summary[1,3]+1;
end;
summary[1,1]=nrow(x);
summary[1,4]=100*summary[1,3]/summary[1,1];
mattrib summary colname={"Total" "Event" "Censored" "% Censored"};
print 'Summary of Number of Event and Censored Times'
 summary;
wi=j(nrow(x),1,1);
mi=j(1,groupnum,0);
di=j(1,groupnum,0);
hi=j(1,groupnum,0);

do i=1 to nrow(x);   
   do j=i+1 to nrow(x);
      if x[i,1] <= x[j,1] then do;
         temp = j(1,ncol(x),0);
	 do k=1 to ncol(x);
	    temp[1,k]= x[i,k];
	    x[i,k]= x[j,k];
	    x[j,k]= temp[1,k];
	 end;
	 if x[i,1] = x[j,1] then do;
	    if x[i,2] > x[j,2] then do;
	       tempvar=j(1,ncol(x),0);
	       do k=2 to ncol(x);
		  tempvar[1,k]= x[i,k]; 
		  x[i,k]= x[j,k];
		  x[j,k]= tempvar[1,k];
	       end;
	    end;
	 end;
      end;
   end;
   if ((i=1) & (x[1,1]^= x[2,1])&(x[1,2]=1)) then d=1;
   if i > 1 then do;
      if x[i,1] ^= x[i-1,1] then do;
         countall=countall+1;
	 if x[i,2]=1 then d=d+1;
      end;
      else if x[i,1]=x[i-1,1] then do;
         if ((x[i,2]^=x[i-1,2]) & x[i-1,2]=0) then d=d+1;
      end;
   end;
   if x[i,2]=1 then di[,x[i,3]]=di[,x[i,3]]+1;
   mi[,x[i,3]]=mi[,x[i,3]]+1;
end;

/* In the above countall is the number of unique times in the data and 
   d is the number of unique event times 
   the above portion of the program also sorts the data so that the obs 
   with the longest time comes first and the obs with the shortest time 
   comes last in the event of a tie a censored observation will come 
   before an event.
*/
					

index=1:nrow(x);
x1=x||t(index);
ipi=j(d,ncol(x)-1,0);
times=j(d,1,0);
ip=j(countall,1,0);
k=1;
if x1[1,1]=1 then do;
   ip[1,1]=1;
   k=2;
end;
do j=2 to nrow(x1);
   if x1[j,1] < x1[j-1,1] then do;
      ip[k,1]=j-1;
      k=k+1;
   end;
   if ((j=nrow(x1)) & ((x1[j,2]=1) | (x1[j,1] ^= x1[j-1,1]))) then ip[k,1]=j;
end;
k=1;
do j=1 to nrow(ip);
   if (x1[ip[j,1],2]=1 & k<= nrow(ipi)) then do;
      ipi[k,1]=ip[j,1];
      k=k+1;
   end;
end;

/* In the above ip is a 1 column matrix that  gives the index of each unique */
/* time ie. what observation from the sorted data set.  		*/
/* ipi is a 1 column matrix that gives  the indes of each unique event time */
/* x1 is a matrix of the sorted data with the time variable in column1 and */
/* the censorvar in column2 there is room in the matrix for covariates all that */
/* needs to be changed is to read them into the original matrix x. these */
/* covariates make up the next columns and the last column of x1 is */
/* a column of the sorted observation numbers.				*/


covmat=j(nrow(x),p,0);
covmat1=j(nrow(x),p,0);
betav=j(101,p,0);
betav[1,]=&initbeta;
do q=1 to p;
  covmat[,q]=x[,q+3];
end;
betamat=j(3,ncol(betav),0);
thetamat=j(3,2,0);
mattrib thetamat format=15.7;
serrmat=j(1,ncol(betav),0);
cntiter=0;
thconmat=j(itend,4,0);

do it= 1 to itend+1;
   cntiter=cntiter+1;
   numiter=0;
   nriter=0;
   wibend=100;
   if it=1 then wibend=1;
   do thetwib = 1 to wibend;
      lambda=.25;
      nriter=0;
      numiter=numiter+1;
      if it ^= 1 then do;
         if it ^= itend then do;
         do j=1 to nrow(x1);
            wi[j,]=((1/theta)+di[,x1[j,3]])/((1/theta)+hi[,x1[j,3]]);
         end;
      end;
         betav[1,]=betacurr;
         betapast=betacurr;
      end;
      neg2llv=j(101,1,0);
      ident=i(p);
      hi=j(1,groupnum,0);
      sumcov=j(nrow(ipi),p,0);
      llam=j(nrow(ipi),1,0);
      clam=j(nrow(ipi),1,0);
      clamvec=j(nrow(x),1,0);
      llamvec=j(nrow(x),1,.000000001);
      loglin=0;
      llpast=0;
   
      do inn= 1 to 100;
         nriter=nriter+1;
         u=j(1,p,0);
         h=j(p,p,0);
         gmat=h;
         llpast=loglin;
         loglin=0;
         do j=1 to nrow(ipi);
            do g=1 to p;
	       sumcov[j,g]=x1[ipi[j,1],g+3];
	       counter=1;
	       do k= ipi[j,1]-1 to 1 by -1;
	          if x1[ipi[j,1],1]=x1[k,1] then do;
	             if x1[k,2]=1 then do;
		        counter=counter+1;
		        sumcov[j,g]=sumcov[j,g]+x1[k,g+3];
	             end;
	          end;
	       end;
	       ipi[j,2]=counter;
	       ipi[j,g+2]=sumcov[j,g];
            end;
            s0=0;
            s1=j(1,p,0);
            s2=j(p,p,0);
            tu=0;
            th=0;
            tloglin=0;
            do mn=ipi[j,1] to 1 by -1;
               currcov=covmat[mn,];
               currbeta=betav[inn,];
	       y=exp((betav[inn,]*t(covmat[mn,]))+log(wi[mn,]));
	       yy=covmat[mn,]*y;
	       yyy=t(covmat[mn,])*covmat[mn,]*y;
	       s0=s0+y;
	       s1=s1+yy;
	       s2=s2+yyy;
            end;
            tu=sumcov[j,]-(counter*s1/s0);
            th=counter*(((t(s1)*s1)-(s2*s0))/(s0*s0));
            tloglin=(betav[inn,]*t(sumcov[j,]))-(counter*log(s0));
            h=h+th;
            u=u+tu;
            loglin=loglin+tloglin;                 
            llam[j,]=counter/s0;
         end;
         if it=1 then do;
            varmat=inv(-h);
            do serrcov=1 to (p);
               serrmat[1,serrcov]=sqrt(varmat[serrcov,serrcov]);
            end;
         end;
         do len=1 to p;
      	       gmat[len,len]=(abs(h[len,len]))##-.5;
         end;
         betacurr=betav[inn,];
         neg2llv[inn,1]=-2*loglin;
         neg2curr=neg2llv[inn,];
         if inn > 1 then do;
            if (max(abs(betav[inn,]-betav[inn-1,]))>&convcrit) then do;
               betav[inn+1,]=betav[inn,]-
                 (u*(gmat*inv((gmat*h*gmat)-(lambda*ident))*gmat));
               if llpast<loglin then lambda=lambda/4;
               else if llpast>loglin then lambda=10*lambda;
            end;
            else do;
      	       inn=100;
      	       betav[inn+1,]=betav[inn,];
            end;
         end;
         else do;
            betav[inn+1,]=betav[inn,]-
               (u*(gmat*inv((gmat*h*gmat)-(lambda*ident))*gmat));
            if llpast<loglin then lambda=lambda/4;
            else if llpast>loglin then lambda=10#lambda;
         end;
      end;
      do ab=1 to nrow(ipi);
         if ab ^=1 then do bc=ipi[ab,1] to ipi[ab-1,1] by -1;
            if x1[bc,1]=x1[ipi[ab,1],1] then do;
               llamvec[bc,]=llam[ab,];
            end;
         end;
         else do bc=ipi[1,1] to 1 by -1;
            if x1[bc,1]=x1[ipi[1,1],1] then do;
               llamvec[bc,]=llam[1,];
            end;
         end;
         do de=ab to nrow(ipi);
            clam[ab,]=clam[ab,]+llam[de,];
         end;
         times[ab,]=x1[ipi[ab,1],1];
      end;
      do je=nrow(ipi) to 2 by -1;
         do ff=ipi[je,1] to ipi[je-1,1]+1 by -1;
            clamvec[ff,]=clam[je,];
         end;
      end;
      do je=ipi[1,1] to 1 by -1;
         clamvec[je,]=clam[1,];
      end;
      lthetbet=0;
      lvec=j(nrow(x),1,0);
      firstsum=j(1,groupnum,0);
      secsum=j(1,groupnum,0);
      if it > 1 then thetacr=theta;
      else theta=0;
      do j=1 to nrow(x);
         firstsum[,x1[j,3]]=firstsum[,x1[j,3]]+
         clamvec[j,]*exp(betacurr*t(covmat[j,]));
         secsum[,x1[j,3]]=secsum[,x1[j,3]]+
            x1[j,2]*((betacurr*t(covmat[j,]))+log(llamvec[j,])); 
         hi[,x1[j,3]]=hi[,x1[j,3]]+clamvec[j,]*exp(betacurr*t(covmat[j,]));
         if it=1 then do;
            lvec[j,]=x1[j,2]*((betacurr*t(covmat[j,]))+log(llamvec[j,]))-
               (clamvec[j,]*exp(betacurr*t(covmat[j,])));
            lthetbet=lthetbet+lvec[j,];
            indlikel=lthetbet;
         end; 
      end;
      if it > 1 then do;
         do i=1 to groupnum;
            lvec[i,]=di[,i]*log(thetacr)-log(gamma(1/thetacr))+
                log(gamma((1/thetacr)+di[,i]))-(((1/thetacr)+
                di[,i])*log(1+thetacr*firstsum[,i])) + secsum[,i];
            lthetbet=lthetbet+lvec[i,];
         end;
      end;
      initl=sum(lvec);
      if it > 1 then do;
         if (max(abs(betacurr-betapast))<&convcrit) then do;
            thetwib=wibend;
         end;
      end;
   end;
   if &options[6,1]={y} then do;
      if max<=1 then do;
         thconmat[it,1]=theta;
         thconmat[it,2]=lthetbet;
         thconmat[it,3]=numiter;
         thconmat[it,4]=nriter;
      end;
      if max=2 then do;
         countmin=countmin+1;
         if countmin=1 then do; 
           tconmat1=j(it-1,4,0);
           mattrib tconmat1 colname=({'Theta' 'Likelihood' 'Wi steps' 'M-C steps'})
                       label={''};
            do thconit = 1 to it-1;
               tconmat1[thconit,]=thconmat[thconit,];
            end; 
          thleft=thetamat[1,1];
          thright=thetamat[3,1];
          mattrib thleft label={' '};
          mattrib thright label={' '};
          print tconmat1;
          print 'Starting Grid Search Theta is between' thleft ' and' thright;
          print /;
         end;
      end;
   end;
   
   if it=1 then do;
      indmat=j(p,4,0);
      do indj=1 to p;
         indmat[indj,1]=betacurr[,indj];
         indmat[indj,2]=serrmat[,indj];
         indmat[indj,3]=indmat[indj,1]##2/indmat[indj,2]##2;
         indmat[indj,4]=1-probchi(indmat[indj,3],1);
      end;
      mattrib indmat colname=({'Estimate' 'S.E.' 'Wald' 'P-value'})
                              rowname=covlist
                              label={'Independence Model'}
                              format=10.5;
      theta=&leftend;
   end;
   if it >1 & it<itend then do;
      if max=0 then do;
         thetamat[1,]=thetamat[2,];
         thetamat[2,]=thetamat[3,];
         thetamat[3,1]=theta;
         thetamat[3,2]=lthetbet;
         betamat[1,]=betamat[2,];
         betamat[2,]=betamat[3,];
         betamat[3,]=betacurr;
         if thetamat[3,2] < thetamat[2,2] then do;
            if it=2 then do;
               thetanow=theta;
               theta=theta+stepsize;
            end;
            else if it=3 then do;
               print indmat;
               print ' ';
               print 'Need to decrease  initial value of theta called leftend in the macro';
               abort;
            end;
            else max=1;
         end;
         if thetamat[3,2]>thetamat[2,2] then do;
            thetanow=theta;
            theta=theta+stepsize;
         end;
      end;
      if max=2 then do;
         if right=1 then do;
            if lthetbet>thetamat[2,2] then do;
               thetamat[1,]=thetamat[2,];
               betamat[1,]=betamat[2,];
               thetamat[2,1]=theta;
               thetamat[2,2]=lthetbet;
               betamat[2,]=betacurr;
            end;
            else do;
               thetamat[3,1]=theta;
               thetamat[3,2]=lthetbet;
               betamat[3,]=betacurr;
            end;       
         end;
         else do;
            if lthetbet>thetamat[2,2] then do;
               thetamat[3,]=thetamat[2,];
               betamat[3,]=betamat[2,];
               thetamat[2,1]=theta;
               thetamat[2,2]=lthetbet;
               betamat[2,]=betacurr;
            end;
            else do;
               thetamat[1,1]=theta;
               thetamat[1,2]=lthetbet;
               betamat[1,]=betacurr;
            end;
         end;
         max=1;
         if thetamat[3,1]-thetamat[1,1]<.00005 then it=itend;
      end;
      if max=1 then do;
         if (thetamat[3,1]-thetamat[2,1]) >= (thetamat[2,1]-thetamat[1,1]) then do;
            right=1;
            left=0;
            thetanow=theta;
            theta=(thetamat[3,1]+thetamat[2,1])/2;
         end;
         if (thetamat[3,1]-thetamat[2,1]) < (thetamat[2,1]-thetamat[1,1]) then do;
            left=1;
           right=0;
            thetanow=theta;
            theta=(thetamat[1,1]+thetamat[2,1])/2;
         end;
         max=2;
      end;
   end;
      if it=itend then theta=thetamat[2,1];
      if it=itend+1 then do;
         x2=x1;
         do i=1 to nrow(x2);
            do j=i+1 to nrow(x2);
               if x2[i,3] >= x2[j,3] then do;
                  temp = j(1,ncol(x2),0);
	          do k=1 to ncol(x2);
	            temp[1,k]= x2[i,k];
	            x2[i,k]= x2[j,k];
	            x2[j,k]= temp[1,k];
	          end;
	          if x2[i,1] = x2[j,1] then do;
	             if x2[i,2] > x2[j,2] then do;
	                tempvar=j(1,ncol(x2),0);
	                do k=2 to ncol(x2);
		           tempvar[1,k]= x2[i,k]; 
		           x2[i,k]= x2[j,k];
		           x2[j,k]= tempvar[1,k];
	                end;
	             end;
	          end;
               end;
            end;
         end;
         do qa=1 to p;
            covmat1[,qa]=x2[,qa+3];
         end;
         infmat=j(1+p+nrow(ipi),1+p+nrow(ipi),0);
         indicat=0;
         do i= 1 to groupnum;
            hscale=0;
            ha=j(nrow(ipi),1,0);
            hb=j(p,1,0);
            hbbmat=j(p,p,0);
            habmat=j(nrow(ipi),p,0);
            haa=0;
            do j=1 to mi[,i];
               indicat=indicat+1;
               hscale=hscale+exp(betacurr*t(covmat1[indicat,]))
                  *clamvec[x2[indicat,ncol(x2)]];
               do m=1 to p;
                  hb[m,]=hb[m,]+covmat1[indicat,m]*exp(betacurr*t(covmat1[indicat,]))*
                     clamvec[x2[indicat,ncol(x2)]];
                  do q=1 to p;
                     hbbmat[m,q]=hbbmat[m,q]+covmat1[indicat,m]*covmat1[indicat,q]*
                        exp(betacurr*t(covmat1[indicat,]))*
                           clamvec[x2[indicat,ncol(x2)]];
                  end;
                  do h= 1 to nrow(ipi);
                     yind=0;
                     if x2[indicat,1] >= x1[ipi[h,1],1] then yind=1;
                     habmat[h,m]=habmat[h,m]+covmat1[indicat,m]*exp(betacurr*t(covmat1[indicat,]))*
                        yind;
                  end;
               end;
               do h=1 to nrow(ipi);
                  yind=0;
                  if x2[indicat,1] >= x1[ipi[h,1],1] then yind=1;
                  ha[h,]=ha[h,]+exp(betacurr*t(covmat1[indicat,]))*yind;
               end;
            end;
            do infj=1 to nrow(infmat);
                  do infk=infj to ncol(infmat);
                     if infj=1 then do;
                        if infk=1 then do;
                           infmat[infj,infk]= infmat[infj,infk]
                                - (3/(theta##3)) 
                                + ((2*log(theta))/(theta##3))
                                + ((2*digamma((1/theta)))/(theta##3))
                                +((trigamma((1/theta)))/(theta##4)) 
                                -((2*digamma((1/theta)+di[,i]))/(theta##3))
                                -((trigamma((1/theta)+di[,i]))/(theta##4))
                                -((di[,i]+(1/theta))/((theta##4)*(((1/theta)+hscale)##2)))
                                +(2/((theta##4)*((1/theta)+hscale)))
                                +((2*(di[,i]+(1/theta)))/((theta##3)*((1/theta)+hscale)))
                                +((2*log((1/theta)+hscale))/(theta##3));

                        end;

                         else if infk<=(p+1) then do;
                            infmat[infj,infk]= infmat[infj,infk]
                              +((di[,i]+(1/theta))*hb[(infk-1),1])
                                /((theta##2)*(((1/theta)+hscale)##2))
                              -((hb[(infk-1),1])/((theta##2)*((1/theta)+hscale)));
                         end;
                             else do;
                                infmat[infj,infk]= infmat[infj,infk]
                                  + (((di[,i]+(1/theta))*ha[infk-(p)-1,1])
                                       /((theta##2)*(((1/theta)+hscale)##2)))
                                  - (ha[infk-(p)-1,1]/((theta##2)*((1/theta)+hscale)));
                              end;

                     end;
                     else if infj<=(p+1) then do;
                        if infk<=(p+1) then do;
                           
                           infmat[infj,infk]= infmat[infj,infk]
                        	-(((di[,i]+(1/theta))*hb[infj-1,1]*hb[infk-1,1])
                        	   /(((1/theta)+hscale)##2))
                        	+(((di[,i]+(1/theta))*hbbmat[infj-1,infk-1])
                        	   /((1/theta)+hscale));
                           je=infmat[infj,infk];
                        end;
                          else do;
                             infmat[infj,infk]= infmat[infj,infk]
                              -(((di[,i]+(1/theta))*ha[infk-1-(p),1]*hb[infj-1,1])
                                   /(((1/theta)+hscale)##2))
                               +(((di[,i]+(1/theta))*habmat[infk-1-(p),infj-1])
                                   /((1/theta)+hscale));
                          end;
                     end;
                     else do;
                        infmat[infj,infk]= infmat[infj,infk]
                           -(((di[,i]+(1/theta))*ha[infj-1-(p),1]*ha[infk-1-(p),1])
                               /(((1/theta)+hscale)##2));
                     end;
                  end;
            end;
         end;
         do infj=1 to d;
                infmat[infj+1+p,infj+1+p]=infmat[infj+1+p,infj+1+p]+((ipi[infj,2])/((llam[infj,1])##2));
         end;
         free x x1 x2 ip covmat covmat1;
         infmat=infmat+t(infmat);
         do in=1 to nrow(infmat);
            infmat[in,in]=infmat[in,in]/2;
         end;
         covarmat=inv(infmat);
         free infmat hscale ha hb habmat;
         semat=j(1,1+p,0);
         do in=1 to (1+p);
            semat[1,in]=sqrt(covarmat[in,in]);
         end;         
         null=j(1,4,0);
         mattrib null rowname=({'-2 Log L'})
                 colname=({'No Frailty' 'Frailty' 'Chi-Square' 'P-value'})
                 label={'Criterion'}
                 format=9.4;
         null[1,1]=-2*indlikel;
         null[1,2]=-2*lthetbet;
         null[1,3]=null[1,1]-null[1,2];
         null[1,4]=1-probchi(null[1,3],1);
         
         print indmat;
         print ' ';
         print ' ';
         print '     Testing Null Hypothesis: Frailty=0'
         null;
         bhat=theta||betacurr;
         final=j(p+1,4,0);
         do infj=1 to p+1;
            final[infj,1]=bhat[,infj];
            final[infj,2]=semat[,infj];
            final[infj,3]=(bhat[,infj]##2)/(semat[,infj]##2);
            final[infj,4]=1-probchi(final[infj,3],1);            
         end;
         mattrib final rowname=fcovlist
                       colname=({'Estimate' 'S.E.' 'Wald' 'Pvalue'})
                       label={'Frailty Model'}
                       format=10.5;
         print ' ';
         print ' ';
         print final;
         kendall=theta/(theta+2);
         mattrib kendall label={''};
         print 'Kendalls Tau ='kendall;
         varcov=j(p+1,p+1,0);
         do infj=1 to p+1;
            do infk=1 to p+1;
               varcov[infj,infk]=covarmat[infj,infk];
            end;
         end;
         mattrib varcov label={"Variance-Covariance Matrix"};
         if &options[1,1]={y} then do;
            print ' ';
            print ' ';
            print varcov;
         end;
         /* Creating the data set with the smaller-variance covariance matrix
            the one with repect to theta and the covariates
            the first column will be the final estimates and the next
            p+1 columns will actually be the  covariances
          */
          if &options[2,1]={y} then do;
             newdat1=t(bhat)||varcov;
             create &data1 from newdat1;
             append from newdat1;
	     close &data1;
          end;
          
          /* Creating a data set that has the estimates of the 
             baseline hazard rates as well as the times that go with
             them
          */
          newdat2=times||llam;
          mattrib newdat2 colname={'Time' 'Baseline Hazard Rate'}
                          label={''};
          if &options[4,1]={y} then do;
             create &data2 from newdat2;
             append from newdat2;
	     close &data2;
          end;

           
          /* Creating the data set with full Covariance matrix
             the first column will be the estimates and the next 1+p+D columns
             will actually be the covariance matrix
           */
          if &options[5,1]={y} then do;
             allestm=bhat||t(llam);
             newdat3=t(allestm)||covarmat;
             create &data3 from newdat3;
             append from newdat3;
	     close &data3;
          end;
          if &options[3,1]={y} then do;
             print /;
             print 'Estimates of the Baseline Hazard Rates';
             print newdat2;
          end;
             
         

      end;
   end;  


%mend;