/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~* * MACRO 'ICSTEST' ------------------------------------------------------------------------------------------------------------------------------------------ The %ICSTEST macro computes the generalized log-rank tests of Zhao & Sun (2004) and Sun, Zhao, and Zhao (2005). The following are arguments listed (separated by commas) enclosed in parantheses. All arguments are optional, except LEFT=, RIGHT=, and GROUP=. DATA= SAS data set to be analyzed. Default is DATA=_last_. LEFT= A numeric variable representing the left endpoint of the time interval. Left-censored observations have a missing value or a value of 0. RIGHT= A numeric variable representing the right endpoint of the time interval. Right-censored observations have a missing value. FREQ= A single numeric variable whose values represent the frequency of occurrence of the observations. GROUP= A variable identifying different treatment groups. ERRORTYPE= Convergence criterion to be used. 1 -- The maximum of the closeness of consecutive estimates, 2 -- The closeness of the log likelihood function, 3 -- The gradient of the log likelihood function, 4 -- The maximum measures of ERRORTYPE=1, ERRORTYPE=2, and ERRORTYPE=3. Default is ERRORTYPE=1. RATECONV= Rate of convergence for the selected ERRORTYPE. Default is RATECONV=1e-7. mRS= Number of resampling for the generalized Greenwood formula. Default is mRS=50. SEED= Random seed for resampling for the generalized Greenwood formula. Default is SEED=8375. *~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/ %macro ICSTEST (version, data=_last_, left= , right= , group= , freq=, errortype=1, rateconv=1e-7, mRS=50, seed= 8375, rho=0, gamma=0) ; %let _version=2.1; %if &version ne %then %put ICSTEST macro Version &_version; %let _opts = %sysfunc(getoption(notes)); options nonotes; /* Check for newer version */ %if %sysevalf(&sysver >= 8.2) %then %do; %let _notfound=0; filename _ver url 'http://ftp.sas.com/techsup/download/stat/versions.dat' termstr=crlf; data _null_; infile _ver end=_eof; input name:$15. ver; if upcase(name)="&sysmacroname" then do; call symput("_newver",ver); stop; end; if _eof then call symput("_notfound",1); run; %if &syserr ne 0 or &_notfound=1 %then %put &sysmacroname: Unable to check for newer version; %else %if %sysevalf(&_newver > &_version) %then %do; %put &sysmacroname: A newer version of the &sysmacroname macro is available.; %put %str( ) You can get the newer version at this location:; %put %str( ) http://support.sas.com/ctx/samples/index.jsp; %end; %end; options &_opts; %if %INDEX(%str(&data), .) = 0 %then %let dataname= WORK.&data ; %else %let dataname = &data ; ************** EXISTENCE for the required variables in LEFT=, RIGHT=, GROUP= ***************; %if %sysfunc(EXIST(&data))=0 %then %do; %put ERROR: The data set %UPCASE(&dataname) does not exist. ; %goto EXIT; %end; %if (%bquote(&left)=) %then %do; %put ERROR: Expecting a name in 'LEFT= '. ; %goto EXIT; %end; %if (%bquote(&right)=) %then %do; %put ERROR: Expecting a name in 'RIGHT= '. ; %goto EXIT; %end; %if (%bquote(&group)=) %then %do; %put ERROR: Expecting a name in 'GROUP= '. ; %goto EXIT; %end; %local savenote; %let savenote=%sysfunc(getoption(notes)); Options nonotes; data _datacheck; set &data; run; Options &savenote; ************** Check if LEFT=, RIGHT=, FREQ=, and GROUP= are in &DATA ***************; %let dsid = %sysfunc(open(%str(_datacheck))); %let xid = %sysfunc(varnum(&dsid, &left)); %if &xid <= 0 %then %do; %put ERROR: Variable %upcase(&left) does not exist.; %let rc1 =%sysfunc( close(&dsid) ); %goto EXIT2; %end; %let _vartype = %sysfunc( vartype(&dsid, &xid) ); %if %upcase(&_vartype)^= N %then %do; %put ERROR: Variable %upcase(&left) is non-numeric.; %let rc1 =%sysfunc( close(&dsid) ); %goto EXIT2; %end; %let xid = %sysfunc(varnum(&dsid, &right)); %if &xid <= 0 %then %do; %put ERROR: Variable %upcase(&right) does not exist.; %let rc1 =%sysfunc( close(&dsid) ); %goto EXIT2; %end; %let _vartype = %sysfunc( vartype(&dsid, &xid) ); %if %upcase(&_vartype)^= N %then %do; %put ERROR: Variable %upcase(&right) is non-numeric.; %let rc1 =%sysfunc( close(&dsid) ); %goto EXIT2; %end; %if (%bquote(&group)^=) %then %do; %let xid = %sysfunc(varnum(&dsid, &group)); %if &xid <= 0 %then %do; %put ERROR: Variable %upcase(&group) does not exist.; %let rc1 =%sysfunc( close(&dsid) ); %goto EXIT2; %end; %end; %if (%bquote(&freq)^=) %then %do; %let xid = %sysfunc(varnum(&dsid, &freq)); %if &xid <= 0 %then %do; %put ERROR: Variable %upcase(&freq) does not exist.; %let rc1 =%sysfunc( close(&dsid) ); %goto EXIT2; %end; %let _vartype = %sysfunc( vartype(&dsid, &xid) ); %if %upcase(&_vartype)^= N %then %do; %put ERROR: Variable %upcase(&freq) is non-numeric.; %let rc1 =%sysfunc( close(&dsid) ); %goto EXIT2; %end; %end; %let _nobs = %sysfunc( attrn(&dsid, NOBS) ); %let _nvars = %sysfunc( attrn(&dsid, NVARS) ); %put NOTE: There were &_nobs observations read from the data set %UPCASE(&dataname).; %let rc1 =%sysfunc( close(&dsid) ); ************** Check if mRS > 1 ***************; %if &mRS<=1 %then %do; %put ERROR:The value for 'mRS=' should be greater than one. ; %goto EXIT; %end; Options nonotes; proc sort data=&data out=indata; by &group; run; /* Delete obs with invalid time values (valid if L=0, R>=0) or invalid FREQ values */ data indata; *drop _rm _zF _zT _bad; retain _rm _zF _zT 0; set indata END=eof; _bad=0; %if %bquote(&freq) ^= %then %do; if &freq < 1 then do; _bad= 1; _zF= 1; end; %end; if (missing(&left)) then &left=0; if (&left < 0) then do; _bad=1; _zT= 2; end; else if missing(&right)=FALSE and &right < &left then do; _bad=1; _zT= 2; end; _rm= _rm + _bad; if eof then do; call symput('_rmN', trim(left(put(_rm,12.)))) ; call symput('_missType', put(_zF+_zT,1.)); end; if _bad=1 then delete; %if %bquote(&freq) ^= %then %do; do i=1 to &freq; output; end; %end; %else output;; run; %let dsid = %sysfunc(open(%str(indata))); %let _nobs = %sysfunc( attrn(&dsid, NOBS) ); %let rc1 =%sysfunc( close(&dsid) ); %if &_nobs=0 %then %do; Options &savenote; %put ERROR: There were no valid observations in the data set %UPCASE(&data).; %goto EXIT; Options nonotes; %end; %else %if &_missType=3 %then %do; Options &savenote; %put NOTE: &_rmN observations deleted due to invalid values in the LEFT=, RIGHT=, or FREQ= variables.; Options nonotes; %end; %else %if &_missType=2 %then %do; Options &savenote; %put NOTE: &_rmN observations deleted due to invalid values in the LEFT= or Right= variables.; Options nonotes; %end; %else %if &_missType=1 %then %do; Options &savenote; %put NOTE: &_rmN observations deleted due to invalid values in the FREQ= variable.; Options nonotes; %end; data indata _GrpName(keep= _Grp &group N) ; set indata END=eof; by &group; retain _Grp N 0; if first.&group=1 then do; N = 1; _Grp+1; end; else do; N + 1; end; output indata; if last.&group=1 then output _GrpName; if eof then call symput('_nGrp', trim(left(put(_Grp,12.)))) ; run; proc print data=_GrpName noobs ; var &group N ; title 'Number of Observations by Group'; run; goption reset=title; %if &_nGrp=1 %then %do; %put ERROR: The GROUP variable should have more than one level.; proc datasets NOLIST NOWARN NODETAILS; delete _datacheck indata _grpname; quit; run; %goto EXIT; %end; proc iml; use indata; read all var{&Left &Right _Grp &group}; L= &Left; R=&Right; maxR = max(R); * handle exact event times ; ExactT= sum(L=R); if ExactT>0 then do; ExactTime = loc(L=R); L[ExactTime, ] = L[ExactTime, ] - 0.00001; end; * assign max(R)*2 for . right-censoring ; MissR = sum(R <= . ); if MissR> 0 then do; locR = loc(R <= .); maxR= max(R); maxL= max(L); if maxR > maxL then maxR2= maxR*2; else maxR2= maxL*2; R[locR, ] = maxR2 ; end; nobs=nrow(L); *Grp: Group Index ; mGrp= &_nGrp; GrpName= &group[1]; do ii= 2 to nobs; if _Grp[ii] ^= _Grp[ii-1] then GrpName = GrpName // &group[ii] ; end; /********************************* construct the non-overlapping intervals (q, p] and determine the number of pseudo-parameters m **********************************/ pp= unique(r); npp= ncol(pp); qq= unique(l); nqq= ncol(qq); q= J(1, npp, .); do; do i=1 to npp; do j=1 to nqq; if (qq[j] < pp[i]) then q[i]= qq[j]; end; if q[i]= qq[nqq] then goto lab1; end; lab1: end; if i > npp then nq= npp; else nq= i; q= unique(q[1:nq]); m= ncol(q); p= J(1, m, .); do i=1 to m; do j= npp to 1 by -1; if (pp[j] > q[i]) then p[i] = pp[j]; end; end; sj = unique(p); m=ncol(sj); alpha= J(nobs, m, 0); do i=1 to nobs; alpha[i, ]= (sj > L[i,1])#(sj <= R[i,1]); end; /******************************************************/ /* Module to calculate isotonic estimator of x /******************************************************/ start isotonic(x,weight); n_X = nrow(x); phat=x; wt=weight; plevel= 1:n_X ; dummy=max(phat)*1000; is_vlt = ((phat//dummy) <(-dummy//phat)) # t(1:(n_X +1)); vlt_loc =min( (^is_vlt )*10*n_X + is_vlt); do while (vlt_loc < 10*n_X); idx= t(loc(phat[1:vlt_loc, ]=phat[vlt_loc-1]) || vlt_loc) ; phat[idx, ] = sum(wt[idx] #phat[idx]) / sum(wt[idx]) ; is_vlt = ((phat//dummy) < (-dummy//phat)) # t(1:(n_X +1)); vlt_loc =min( (^is_vlt )*10*n_X + is_vlt); end; return(phat); finish isotonic; if m > 1 then do; error&errortype=1; epsilon= 0.2; * STEP 1: set up the initials for beta; iter=1; beta_old=J(m-1,1,1/m)#t(1:m-1); phat_old= (beta_old//1) - (0//beta_old); loglik= sum(log(alpha*phat_old)) ; alpha_dff = (alpha[ ,1:m-1]-alpha[ ,2:m]) ; /******************************************************/ /* Module for EM-Iterative Convex Minorant Algorithm /******************************************************/ do while (error&errortype > &rateconv); * STEP 2: Apply ICM ; ap= alpha*phat_old ; if sum(ap =0) > 0 then ap [loc(ap =0)] =0.00001; d1l_n= alpha_dff / repeat(ap, 1, m-1); d1l= t(d1l_n[+, ]); d2l_n= - (d1l_n##2); d2l=t(d2l_n[+, ]); W=diag(-d2l); x_new= beta_old + d1l/vecdiag(W); if sum(x_new<0) then x_new[loc(x_new<0)]=0.00001; beta_up = isotonic(x_new, vecdiag(W)) ; phat_up = (beta_up//1) - (0//beta_up); ap_up = alpha*phat_up; if sum(ap_up=0) > 0 then ap_up[loc(ap_up=0)] =0.00001; loglik_up = sum(log(ap_up)) ; /*** line search ***/ if loglik_up >= loglik then beta_new= beta_up; if loglik_up < loglik then do; lambda=1; gamma= 0.5; zz= beta_up; phat_old= (beta_old//1) - (0//beta_old) ; ap_old = alpha*phat_old ; if sum(ap_old=0) > 0 then ap_old[loc(ap_old=0)] =0.00001; d1l_n_old= alpha_dff / repeat(ap_old , 1, m-1); d1l_old= t(d1l_n_old[+, ]); loglik_old= sum(log(alpha*(phat_old))) ; loglik_zz= loglik_up ; lb = loglik_old + epsilon*t(d1l_old) *(zz-beta_old); ub= loglik_old + (1-epsilon)*t(d1l_old) *(zz-beta_old) ; do while ( (loglik_zz < lb) | (loglik_zz > ub) ); if loglik_zz < lb then lambda= lambda-gamma; if loglik_zz > ub then lambda= lambda+gamma; zz= beta_old + lambda*(beta_up - beta_old); gamma = gamma/2 ; ap_zz = alpha*( (zz//1) - (0//zz) ); if sum(ap_zz=0) > 0 then ap_zz[loc(ap_zz=0)] =0.00001; loglik_zz= sum( log(ap_zz) ) ; lbtemp = loglik_old + epsilon*t(d1l_old) *(zz-beta_old); ubtemp= loglik_old + (1-epsilon)*t(d1l_old) *(zz-beta_old) ; lbtemp = loglik_old + epsilon*t(d1l_old) *(zz-beta_old); ubtemp= loglik_old + (1-epsilon)*t(d1l_old) *(zz-beta_old) ; lb= min(lbtemp, ubtemp); ub=max(lbtemp, ubtemp); end; beta_new = zz; end; beta_up = beta_new; * STEP 3: Apply Self-consistency ; phat_old= (beta_up //1) - (0//beta_up ); ap = alpha*phat_old; if sum(ap =0) > 0 then ap [loc(ap =0)] =0.00001; p_temp= (alpha#t(phat_old))/repeat(ap , 1, m); p_new = t(p_temp[+, ]/nobs); beta_new=cusum(p_new)[1:m-1]; * STEP 4: Convergence checking ; p_new = (beta_new //1) - (0//beta_new); ap_new = alpha*p_new ; if sum(ap_new =0) > 0 then ap_new [loc(ap_new =0)] =0.00001; loglik_new= sum( log(ap_new) ) ; d1l_n= alpha_dff / repeat(ap_new, 1, m-1); d1l= t(d1l_n[+, ]); error2= max(abs(loglik_new - loglik)); error1= max(abs(beta_new - beta_old)); error3 = sum(beta_new#d1l); error4 = max(error1, error2, abs(error3)); summary_itertp = iter || loglik_new || error1 || error2 || error3 ||error4 ; if iter=1 then summary_iter = summary_itertp ; else summary_iter = summary_iter // summary_itertp; if error&errortype > &rateconv then do; beta_old = beta_new; iter= iter +1; phat_old = p_new; loglik=loglik_new; end; end; beta_new = beta_new // 1; surv = 1-beta_new ; phat= beta_new - (0//beta_new[1:m-1]); end; if m=1 then do; phat=1; beta_new=1; end; /******************************************************/ /* Generalized Log-rank Test I (Zhao and Sun, 2004) /******************************************************/ deltai =repeat(R <= maxR, 1, m) ; rhoij =(1-deltai) # (repeat(L, 1, m) >= repeat(sj[1, 1:m], nobs,1) ); ap = alpha*phat; if sum(ap<=0) >0 then do; dj=J(nobs, m, 0); nonzeroi = loc(ap>0); dj[nonzeroi, ] = deltai[nonzeroi, ] #alpha[nonzeroi, ] #repeat(t(phat), nobs,1)[nonzeroi, ] / repeat(ap[nonzeroi, ] , 1, m); end; else do; dj = deltai#alpha#repeat(t(phat), nobs,1)/ repeat(ap, 1, m); end; dj = dj[+, ]; nj = J(1, m, 0); do j=1 to m; nj[ ,j]= sum(dj[ ,j:m]) + rhoij [+, j]; end; dj = dj[1, 1:m]; djk = J(mGrp, m, 0); njk = J(mGrp, m, 0); do k=1 to mGrp; indGrpk = loc(_Grp =k); alphak= alpha[indGrpk, ]; deltaik= deltai[indGrpk, ]; apk =alphak*phat; if sum(apk<=0) >0 then do; dj_temp =J(ncol(indGrpk), m, 0); nonzeroi = loc(apk>0); dj_temp[nonzeroi, ] = deltaik[nonzeroi, ] #alphak[nonzeroi, ] #repeat(t(phat), ncol(indGrpk),1)[nonzeroi, ] / repeat(apk[nonzeroi, ] , 1, m); end; else do; dj_temp = deltaik#alphak#repeat(t(phat), ncol(indGrpk),1)/ repeat(apk, 1, m); end; djk[k, ] = dj_temp[+, ]; do j=1 to m; njk[k, j] = sum(djk[k ,j:m]) + sum(rhoij[indGrpk, j]); end; end; Uk = J(mGrp, m, 0); nonzero = loc(nj > 0); if sum(nj>0) > 0 then Uk= djk[ ,nonzero] - njk[ ,nonzero]#repeat(dj[ ,nonzero], mGrp,1)/repeat(nj[ ,nonzero], mGrp,1); Uk = Uk[ ,+] ; /*** Generalized Greenwood formula based on resampling ***/ mRS= &mRS; Uk_RS = J(mGrp, mRS, 0); V_RS1 = J(mGrp, mGrp, 0); do b=1 to mRS; Tib= J(nobs, 1, 0); deltaib = deltai[ , 1]; do ii=1 to nobs; if deltai[ii, 1]=0 then do; Tib[ii, ]=L[ii , ]; end; else if deltai[ii, 1]=1 then do; Li= sum(sj <= L[ii ,1]) ; if Li=0 then SLi=1; else SLi= 1 - beta_new[Li, ]; Ri= sum(sj <= R[ii ,1]) ; SRi= 1 - beta_new[Ri, ]; fi = (sj>L[ii ,1])#(sj<= R[ii ,1])#t(phat) / (SLi - SRi); positivef = loc(fi>0); fi= fi[ ,positivef] / sum(fi[ ,positivef]); * to make sure sum(fi)=1::: required to run "RANDMULTINOMIAL" ; call randseed(&seed); Tib[ii, ] = sj[ ,sum(RANDMULTINOMIAL(1, 1, fi) # positivef)]; end; end; db = J(1, m, 0); *dr[m]=0 always; nb= J(1, m, 0); *nr[m]=0 always; dkb=J(mGrp, m, 0); nkb= J(mGrp, m, 0); do jj=1 to m; db[ ,jj] = sum((Tib#deltaib) = sj[ ,jj]); nb[ ,jj] = sum(Tib >= sj[ ,jj]); do k=1 to mGrp; indGrpk = loc(_Grp=k); dkb[k, jj] = sum((Tib[indGrpk, ]#deltaib[indGrpk, ]) = sj[ ,jj]); nkb[k, jj] = sum(Tib[indGrpk, ] >= sj[ ,jj]); end; end; Vb= J(mGrp, mGrp, 0); *** if nj<=1 then Vbj=0 ; if sum(nb<= 1) > 0 then do; Vb_temp= J(1, m, 0); posloc= loc(nb> 1) ; Vb_temp[ ,posloc]= db[ ,posloc] # (nb[ ,posloc]-db[ ,posloc]) / ((nb[ ,posloc]##2) # (nb[ ,posloc]-1)) ; end; else do; Vb_temp = db # (nb-db) / ((nb##2) # (nb-1)) ; end; do k1=1 to mGrp; do k2=k1 to mGrp; if k2= k1 then Vb[k1, k2] = sum(nkb[k2, ] # (nb-nkb[k2, ]) # Vb_temp); if k2> k1 then do; Vb[k1, k2] = sum(-nkb[k1, ] # nkb[k2, ] # Vb_temp) ; Vb[k2, k1] = Vb[k1, k2] ; end; end; end; V_RS1 = V_RS1 + Vb ; if sum(nb< 1) > 0 then do; Ukb = J(mGrp, m, 0); posloc0= loc(nb> 0) ; Ukb[ ,posloc0] = dkb[ ,posloc0] - nkb[ ,posloc0] #repeat(db[ ,posloc0] , mGrp,1) / repeat(nb[ ,posloc0] , mGrp,1) ; end; else do; Ukb = dkb - nkb#repeat(db, mGrp,1) / repeat(nb, mGrp,1) ; end; Uk_RS[ ,b] = Ukb[ ,+] ; end; *end for b=1 to mRS; V_RS1 = V_RS1/mRS; Uk_RSm= Uk_RS[ ,:]; V_RS2 = (1+1/mRS)*(Uk_RS-repeat(Uk_RSm, 1,mRS))*t(Uk_RS-repeat(Uk_RSm, 1,mRS)) / (mRS-1); V_RS = V_RS1 + V_RS2; ****** chisquare test *******; *chisq= t(uvec) * ginv(sigma) * uvec; x= (V_RS || Uk) // (t(Uk)||{0}); a=sweep(x, 1:mGrp); df=0; do i=1 to mGrp; if a[i,i] > 0 then df =df+1; end; if df >0 then do; chisq= -a[mGrp+1,mGrp+1]; pvalue= 1.- probchi(chisq, df); end; else do; chisq= .; pvalue= .; end; /******************************************************/ /* Generalized Log-rank Test II (Sun, Zhao and Zhao, 2005) /******************************************************/ /* Sun and Zhao(2005) p 51 */ fl= j(nobs, 1, 0); /* \hat{F}_n(L_i) */ fr= j(nobs, 1, 0); /* \hat{F}_n(R_i) */ etafl=j(nobs, 1, 0); /* \eta(\hat{F}_n(L_i) */ etafr=j(nobs, 1, 0); /* \eta(\hat{F}_n(R_i) */ k_n= j(nobs, 1, 0); /* K_n(U_i,V_i,\Delta_i,\Gamma_i) */ do i= 1 to nobs; do j=1 to m; if sj[j] <= L[i] then fl[i]=beta_new[j]; end; do j=1 to m; if sj[j] <= R[i] then fr[i]=beta_new[j]; end; den= fr[i] - fl[i]; if abs(den) < 1e-8 then K_n[i]= 0; else do; a= fl[i]; if a<1e-8 then etafl[i]=1; else do; if (&rho>0 | &gamma>0) then etafl[i]= 1 - ((1-a)**(1+&rho)) * log(1-a) * (a**&gamma); else do; a=1-a; etafl[i]= 1 - a*log(a); end; end; b= 1-fr[i]; if b<1e-8 then etafr[i]= 1; else do; if (&rho>0 | &gamma>0) then etafr[i]= 1- (b**(1+&rho)) * log(b) * ((1-b)**&gamma); else etafr[i]= 1- b*log(b); end; k_n[i]= (etafr[i] - etafl[i]) / den; end; end; *************** Test Statistic uvec ****************; uvec= j(mGrp, 1, 0); nvec= j(mGrp, 1, 0); do k=1 to mGrp; indGrp = loc(_Grp=k); nvec[k]= ncol(indGrp); uvec[k]= sum(k_n[indGrp]); end; **************** Covariance Matrix of uvec ***********; /* Q_0(K_0^2) S, J & J(2005) Theorem 1 p52 */ sigma0= sum(k_n#k_n)/nobs; sigma= sigma0#(diag(nvec) - nvec*t(nvec)/nobs); ****** chisquare test *******; *chisq= t(uvec) * ginv(sigma) * uvec; x2= (sigma || uvec) // (t(uvec)||{0}); a2=sweep(x2, 1:mGrp); df2=0; do i=1 to mGrp; if a2[i,i] > 0 then df2=df2+1; end; if df2>0 then do; chisq2= -a2[mGrp+1,mGrp+1]; pvalue2= 1.-probchi(chisq2, df2); end; else do; chisq2=.; pvalue2=.; end; /******************************************************/ /* output /******************************************************/ reset noname; print ,,, "Generalized Log-rank Test I (Zhao & Sun, 2004)", ,, "Test Statistic and Covariance Matrix", , GrpName[label="&group"] Uk[label='U' format=12.4] V_RS[label='cov(U)' format=12.4] ,,chisq[label= 'Chi-Square' format=12.4] df[label= 'DF' format=5.0] pvalue[label= 'Pr > Chi-Square' format=pvalue6.4] ; print ,,,"Generalized Log-Rank Test II (Sun, Zhao & Zhao, 2005)", %if &rho=0 %then %do; %if &gamma=0 %then "xi(x)=xlog(x)"; %else "xi(x)=xlog(x)*(1-x)**&gamma"; %end; %else %do; %if &gamma=0 %then "xi(x)=xlog(x)*x**&rho"; %else "xi(x)=xlog(x)*x**&rho*(1-x)**&gamma"; %end; ,, "Test Statistic and Covariance Matrix",, GrpName[colname={"&group"}] uvec[colname={'U'} format=10.4] sigma[colname={'cov(U)'} format=10.4]; print , chisq2[colname={'ChiSquare'} format=12.4] df2[colname={'DF'} format=5.0] pvalue2[colname={"Pr>ChiSquare"} format=pvalue6.4]; if (m=1) then print ,,, "WARNING: The test results shown may be unreliable."; quit; proc datasets NOLIST NOWARN NODETAILS; delete _datacheck indata _grpname; quit; run; Options &savenote; %put NOTE: ICSTEST macro completed.; %EXIT2: Options nonotes; proc datasets NOLIST NOWARN NODETAILS; delete _datacheck; quit; run; Options &savenote; %EXIT: %mend ICSTEST;