SNAFTM.sas

/*code to simulate data according to a SNAFTM*/

options nomprint nosymbolgen nomlogic nonotes;


%macro simulate(
subjects=2500,            /* number of subjects in the simulated data sets */
    tpoints= 10,            /* maximum number of measurements for each subject in the simulated data */
    
           
   psi1= 0.3,      /* true psi*/
	
  gamma0 = log(3/7), gamma1 = 2, gamma2 = log(0.5), gamma3 = log(1.5),
                              /* gamma =parameters for L*/
	
   alpha0= log(2/7) , alpha1=0.5 , alpha2 = 0.5, alpha3 = log(4),                      
                              /* alpha = parameters for A*/
	

	cut=30,

	lam=0.01,/*scale parameter for weibull*/

	shape=1,/*shape parameter for weibull*/

	dataout= Simdata_all
);


     data &dataout;
       length  A Am1 L    Am1L  censor1 censor2 Y Ym eligible  maxT T 3 ;

         

        	 sT0 = %sysevalf(&sT0) ;
             sL = %sysevalf(&sL) ;
             sA = %sysevalf(&sA) ;
             sY = %sysevalf(&sY);
                 
        
          
                              
               do id=1 to %sysevalf(&subjects); 
					call ranexp(sT0,expvar); 
                    T0=(expvar/(&lam**&shape))**(1/&shape); 

					
					if T0<&cut then IT0=1;
					else IT0=0;    
					
					maxT=&tpoints;
                    initA = .;
                    C1time = .; 
                    C2time = . ;
                    array Larray L1-L&tpoints;
                    array Aarray A1-A&tpoints;
                    array Yarray Y1-Y&tpoints;
					array Ymarray Ym1-Ym&tpoints;
                    array pAarray pA1-pA&tpoints;
                    array pLarray pL1-pL&tpoints;
					array chkarray chk1-chk&tpoints;
					array tarray t1-t&tpoints;
                    
                   
                    
                    array Censearray Cense1-Cense&tpoints;
                    array Censarray  Cens1-Cens&tpoints;

                    array pcensarray pCens1-pCens&tpoints;
                    array pcensearray pCense1-pCense&tpoints;

                    time = 1.0;
                    j = 1;

					Ym1=0;
					
                    logitL = &gamma0 + &gamma1 * IT0 ;
					
					
                    pL = 1/(1+(1/exp(logitL)));
                    pLarray(1) = pL;
                    
                    if pL=1 then do; 
                         L1 = 1;
                    end;
                    else do;                      
                         call ranbin(sL,1,pL,L1);                 
                    end;                             
                                        
                    logitA = &alpha0 + &alpha1 * L1 ;                    
                    pA1 = 1/(1+(1/exp(logitA)));
					
                    
                                    
                         call ranbin(sA,1,pA1,A1);                      
                  
                
                    
					/*generate Y's*/
					chk1 = exp(&psi1*A1);

					if T0 > chk1 then do;
						Y1=0;
						t1=.;
					end;
					else if T0 <= chk1 then do;
						Y1=1;
						t1=T0*exp(-&psi1*A1);
					end;
					
      
                    cens_1 = 0; 
                    cense = 0;
                    j = 1;                    
                   cens = 0; ;
                    pcens = 0;
                    Censearray(1) = cense;
                    Censarray(1) = cens; 
                    cens_1 = cens;
                    pcensarray(1) = pCens; 
                   
           
           
                    do j=2 to &tpoints;
                         time = j;
                         if j = C2time + 1 then leave ; 
                         if j = C1time + 1 then leave;  
						

						 Ymarray(j)=Yarray(j-1);
						 

						
                         logitL = &gamma0 + &gamma1 * IT0 + &gamma2 * Aarray(j-1) + &gamma3 * Larray(j-1);
						
						pL=exp(logitL)/(1+exp(logitL));
                         pLarray(j) = pL;
                                          
                         if pL = 1 then do; 
                             Larray(j) = 1;
                         end;
                         else do;                          
                              call ranbin(sL,1,pL,Larray(j));                            
                         end;
                     
 						logitA = &alpha0 + &alpha1 * Larray(j) + &alpha2 * Larray(j-1) + &alpha3 * Aarray(j-1);
                        
                         pAarray(j)=exp(logitA)/(1+exp(logitA));
						
                         pp= pAarray(j);                                                
                                                     
                              call ranbin(sA,1,pAarray(j),Aarray(j));                                 
                    
                         if Aarray(j)=1 then timeA=j;
                   
						 if Ymarray(j)=1 then do;
							Aarray(j)=0;
							Larray(j)=0;
						 end;
						
						/*generate Y's*/
						chkarray(j) = chkarray(j-1) + exp(&psi1*Aarray(j));

						if T0 > chkarray(j) then do;
							Yarray(j)=0;
							tarray(j)=.;
							
						end;
						else if T0 <= chkarray(j) then do;
							Yarray(j)=1;
							if Ymarray(j)=1 then do;
							tarray(j)=tarray(j-1);
							end;
							else do;
							tarray(j)=(j-1) + ((T0-chkarray(j-1))*exp(-&psi1*Aarray(j)));
							end;
							
						end;
						
                        cens = 0; ;
                        pcens = 0;
                         
                         Censearray(j) = cense; 
                         Censarray(j) = cens; 
                         pCensarray(j ) = pCens;
                         pCensearray(j) = pCense ;
                         cens_1 = cens;

					
                    end; /* end j = 2 to &tpoints */
                    if C1time = . then C1time = &tpoints ; 
                    if C2time = . then C2time = &tpoints ;
                     
                    time_used = min(C1time,C2time);
                    do tpoint=1 to  time_used; 
						  
                         tpoint2 = tpoint-1;                         
                         L=Larray(tpoint);
						 chk=chkarray(tpoint); 
                                             pL = pLarray(tpoint);                          
                                             A=Aarray(tpoint);
                         
                         _pA_t=pAarray(tpoint);                                                      
                         if tpoint=1 then Am1=0; else Am1=Aarray(tpoint-1);
                         if tpoint=1 then Lm1=0; else Lm1=Larray(tpoint-1);
                         Am1L = Am1*L;
                         
                    
                         censor1 = Censarray(tpoint);
                         pC_t = 1 - pCensarray(tpoint) ; /* want p(cens = 0 ) = 1 - p(cens = 1) */

                         censor2 = Censearray(tpoint);
                         pC2_t = 1 - pCensearray(tpoint);

                         
						 	  Ym = Ymarray(tpoint);
                              if time_used < &tpoints then do;
                                   if tpoint < time_used then Y = Yarray(tpoint);
                                   else Y = . ;
                              end;
                              else if time_used = &tpoints then Y = Yarray(tpoint) ;
							  T=tarray(maxT);
                    
                             
                        
                         output;

                    end; /* end of tpoint loop */
               end;  /* end of id loop */
          

          keep id  tpoint  tpoint2 T0 IT0 chk Y Ym  A Am1 L Lm1 _pA_t T maxT 
               pL  ;

		
               call symput('sT0',sT0);
               call symput('sA',sA);
               call symput('sL',sL);
              
               call symput('sY',sY);
               
        
 
           run;


			

		   data &dataout;
				set &dataout;
				if Ym=0;
			run;

  			



%mend simulate;



         %let sT0 = 123;
          %let sL =  1217607105;
          %let sA =  1920525422;
          %let sY = 2;
          


          %simulate();