/***************************************************************************************/
/* PROGRAM: schacht.SAS                                                                 */
/*                                                                                     */
/* DESCRIPTION : ANALYZE DATA ON A NONPARAMETRIC MANNER                                */
/*                                                                                     */
/* INPUT : DATA=  DATA SET WITH VARIABLE OF INTEREST                                   */
/*                REQURIED ARGUMENT                                                    */
/*                                                                                     */
/*         OUT=   OUTPUT DATA SET                                                      */
/*                OPTIONAL ARGUMENT                                                    */
/*                                                                                     */
/*         GROUP= CATEGORICAL VARIABLE INDICATING THE DIFFERENT GROUPS                 */
/*                REQURIED ARGUMENT                                                    */
/*                                                                                     */
/*         RESP=  RESPONSE VARIABLE                                                    */
/*                REQURIED ARGUMENT                                                    */
/*                                                                                     */
/*         COV=   COVARIATES FOR WHICH A CORRECTION WILL BE MADE                       */
/*                REQURIED ARGUMENT                                                    */
/*                                                                                     */
/*         CONTRAST= CONTRAST STATEMENT FOR TESTING HYPOTHESIS                         */
/*                   REQURIED ARGUMENT                                                 */
/*                   BY DEFAULT=%STR(1 -1)                                             */                                                                   */
/*                                                                                     */
/*         ALPHA=  ALPHA LEVEL FOR CONFIDENCE INTERVAL FOR 2 GROUPS                    */
/*                 BY DEFAULT=0.05                                                     */
/*                 OPTIONAL ARGUMENT                                                   */
/*                                                                                     */
/*         WHERE= WHERE STATEMENT APPLIED TO DATA                                      */
/*                OPTIONAL ARGUMENT                                                    */
/*                                                                                     */
/*  EXAMPLE:                                                                           */
/*    %NONPAR(DATA=mydata, GROUP=tx, RESP=D30, COV=age sex,CONTRAST=%str(1 -1));       */
/*                                                                                     */
/*  NOTE: 1. THE FOLLOWING TEMPORAY DATA SETS WILL BE CREATED AND DELETED              */
/*           _DATA,                                                                    */
/*           _GROUP.                                                                   */
/*         2. THE MACRO WILL WORK FOR SAS VERSION 8.2                                  */
/*                                                                                     */
/*                                                                                     */
/* OUTPUT : SAS DATA SET NAMED BY THE OUT= ARGUMENT                                    */
/*                                                                                     */
/* CREATED BY: KRIS BOGAERTS                                                           */
/*                                                                                     */
/* CREATION DATE : 28JAN2004                                                           */
/*                                                                                     */
/* VALIDATION BY                                                                       */
/*                                                                                     */
/* ADJUSTMENT: BY K. BOGAERTS ON 03MAR04: Add confidence intervals for 2 groups        */
/* ADJUSTMENT: BY K. BOGAERTS ON 09JAN06: new calculation for gamma                    */
/* ADJUSTMENT: BY K. BOGAERTS ON 11JAN06: build in protection for overflow by          */
/*                                        calculation of CI for t-approximation        */                                                                              */
/*                                                                                     */
/***************************************************************************************/;
%macro schacht(data=,
               out=,
               group=,
               resp=,
               cov=,
               contrast=%str(1 -1),
               alpha=0.05,
               where=);
/* checking for where criterium */
/*------------------------------*/
%local s;
%if %length(&where)>0 %then
  %do;
   %let s=WHERE &where;
  %end;
%else %do; 
   %let s=;
  %end;

/* sort data by group and   */
/* create temporary dataset */
/*--------------------------*/
proc sort data=&data out=_data;
by &group;
run; 

/* remove missing observations */
data _data;
 set _data;
if nmiss(of &group &resp &cov)>0 then delete;
&s;
run;

/* Counting the number of groups */
/*-------------------------------*/
proc sql noprint;
%local ng;
select count(distinct &group) as ng into: ngroups from _data;

/* Get all values of grouping variable*/
/*------------------------------------*/
create table _group as
select distinct &group 
from _data;
quit;

/* Recoding the grouping values to 1,2,... */
/*-----------------------------------------*/
data _group;
 set _group;
newgroup=_N_;
run;

data _data;
 merge _data _group;
by &group;
run;

/* Counting the numbers of observations per group */
/*------------------------------------------------*/
proc sql noprint;
%do i=1 %to &ngroups;
%local n&i;
select count(newgroup) into: n&i from _data where newgroup=&i;
%end;
%local n;
select count(newgroup) into: n from _data;
quit;

proc iml;

/* function VAR */
/*--------------*/
START var(x);
   n=nrow(x);
   return((x-REPEAT(x[:,],n,1))[##,]/(n-1));
FINISH;

/* function COVAR */
/*----------------*/
START covar(x,y);  
   C=NCOL(x);
   N=NROW(x); 
   free covariance;
   do i=1 to C;
    covardata=x[,i]||y[,i];
    covardata=covardata-repeat(covardata[:,],N,1); 
    XPX=covardata`*covardata/(N-1); 
    covariance=covariance||XPX[1,2];
   end;
   return(covariance);                                                                                          
FINISH;  

/* read in data */
/*--------------*/
use _data;
read all var {newgroup} into group;
read all var {&resp} into resp;
read all var {&cov} into cov;
close _data;

/*Number of covariates*/
/*--------------------*/
ncov=ncol(cov);
call symput('ncov',char(ncov));

/*Total Number of observations*/
/*----------------------------*/
nobs=nrow(resp);

/*Calculating ranks */
/*------------------*/;
/* Overall Ranks for Response*/
RR=ranktie(resp);
/* Overall Ranks for Covariates*/
free RC;
%do i=1 %to &ncov;
RC&i=ranktie(cov[,&i]);
RC=RC||RC&i;
%end;

/*all except one group ranks combined with within group ranks for response*/
%do i=1 %to &ngroups;
  RR_&i=j(nobs,1,.);
  RR_&i[loc(group=&i)]=ranktie(resp[loc(group=&i),]);
  RR_&i[loc(group^=&i)]=ranktie(resp[loc(group^=&i),]);
%end; /* loop i */
/*all except one group ranks combined with within group ranks for covariates*/
  %do i=1 %to &ngroups;
   free RC_&i;
   %do j=1 %to &ncov;
    RC&j._&i=j(nobs,1,.);
    RC&j._&i[loc(group=&i)]=ranktie(cov[loc(group=&i),&j]);
    RC&j._&i[loc(group^=&i)]=ranktie(cov[loc(group^=&i),&j]);
    RC_&i=RC_&i.||RC&j._&i;
  %end; /* loop j */
%end; /* loop i */


/* Calculating F_i'(Y_ik) */
%do i=1 %to &ngroups;
     F_&i=(RR-RR_&i)/&&n&i;
%end;

/* Calculating H_i(Y_ik) */
%do i=1 %to &ngroups;
 H&i=0;
 %do j=1 %to &ngroups;
  %if &i^=&j %then %do;
    H&i=H&i + F_&j[loc(group=&i),];
  %end;
  H&i=H&i./(&ngroups-1);
 %end;
%end;

/* Calculating G_c,i'(Y_ik) */
%do i=1 %to &ngroups;
     G_&i=(RC-RC_&i)/&&n&i;
%end;

/* Calculating K_c,i(Y_ik) */
%do i=1 %to &ngroups;
 K&i=0;
 %do j=1 %to &ngroups;
  %if &i^=&j %then %do;
    K&i=K&i + G_&j[loc(group=&i),];
  %end;
   K&i=K&i./(&ngroups-1);
 %end;
%end;


/* Calculating p's    */
/*--------------------*/
%do i=1 %to &ngroups;
p&i=H&i[:,]; 
*print p&i;
%end;

/* Calculating q's    */
/*--------------------*/
%do i=1 %to &ngroups;
q&i=K&i[:,]; 
*print q&i;
%end;

/* calculation of gamma's */
/*------------------------*/
/* Covariance H and K */
covHK=0;
%do i=1 %to &ngroups;
covHK=covHK+covar(repeat(H&i,1,&ncov),K&i)/&&n&i;
%end;

/* Covariance F and G */
covFG=0;
%do i=1 %to &ngroups;
  %do j=1 %to &ngroups;
   %if &i^=&j %then %do;
     covFG=covFG+ covar(repeat(F_&j[loc(group=&i),],1,&ncov),G_&j[loc(group=&i),])/&&n&i;
  %end;
 %end;
%end;
covFG=covFG/((&ngroups-1)**2);

/* Variance for K */
VarK=0;
%do i=1 %to &ngroups;
VarK=varK+var(K&i)/&&n&i;
%end;

/* Variance for G */
VarG=0;
%do i=1 %to &ngroups;
 %do j=1 %to &ngroups;
   %if &i^=&j %then %do;
     VarG=varG+var(G_&j[loc(group=&i),])/&&n&i;
   %end;
 %end;
%end;
VarG=VarG/((&ngroups-1)**2);


* START ADDED BY KB ON 09JAN06;
/* Covariances */
 /* Covariance K and K */
free covKK;
%do q=1 %to &ncov;
 covKqK=0;
 %do p=1 %to &ngroups;
  covKqK=covKqK+covar(repeat(K&p[,&q],1,&ncov),K&p)/&&n&p;
 %end;
covKK=covKK//covKqK;
%end;
print covKK;

/* Covariance G and G */
free covGG;
%do q=1 %to &ncov;
 covGqG=0;
 %do i=1 %to &ngroups;
   %do j=1 %to &ngroups;
    %if &i^=&j %then %do;
      covGqG=covGqG+ covar(repeat(G_&j[loc(group=&i),&q],1,&ncov),G_&j[loc(group=&i),])/&&n&i;
   %end;
  %end;
 %end;
 covGqG=covGqG/((&ngroups-1)**2);
 covGG=covGG//covGqG;
%end;


/*create variance matrix V22 */
V22=covKK+covGG;
/*create variance matrix V12 */
V12=covHK+covFG;

/* calculation of gamma */ 
gamma=solve(V22,V12`)`;
*gamma=gamma#0; /* for unadjusted rel eff. */
* END ADDED BY KB ON 09JAN06;
%do j=1 %to &ncov;
gamma&j=gamma[&j];
%end;
print gamma;

/*Relative effects*/
free releff;
%do i=1 %to &ngroups;
releff&i=p&i - gamma*(q&i-0.5)`;
releff=releff||releff&i;
print releff&i;
%end;


/* Variance matrix */
/*-----------------*/
/* Variances*/
variances=j(&ngroups,1,.);
%do i=1 %to &ngroups;
variances[&i]=(&n./&&n&i.)#var(H&i - (gamma#K&i)[,+]);
 %do j=1 %to &ngroups;
  %if &i^=&j %then %do;
variances[&i]=variances[&i] + (&n./(&&n&j.*(&ngroups-1)##2))*var(F_&i[loc(group=&j),]-(gamma#G_&i[loc(group=&j),])[,+]);
 %end;
%end;
%end;/* end i loop */

/* Covariances */
covariances=j(&ngroups,&ngroups,0);
%do i=1 %to &ngroups;
 %do j=(&i+1) %to &ngroups;
covariances[&i,&j]=-&n./(&&n&i.#(&ngroups-1))*covar(H&i-(gamma#K&i)[,+],F_&j[loc(group=&i),]-(gamma#G_&j[loc(group=&i),])[,+]) -
                   &n./(&&n&j.#(&ngroups-1))*covar(H&j-(gamma#K&j)[,+],F_&i[loc(group=&j),]-(gamma#G_&i[loc(group=&j),])[,+]); 
    %do k=1 %to &ngroups;
      %if (&k^=&i) & (&k^=&j) %then %do;
covariances[&i,&j]=covariances[&i,&j] + 
                   (&n./(&&n&k.*(&ngroups-1)##2))#covar(F_&i[loc(group=&k),]-(gamma#G_&i[loc(group=&k),])[,+],F_&j[loc(group=&k),]-(gamma#G_&j[loc(group=&k),])[,+]);
	 %end; /* end if k */
    %end; /* end loop k*/ 
 %end; /* end loop j*/
%end;/*end loop i */

/*create variance matrix */
V=diag(variances) + covariances + covariances`;

/* calculate statistics */
/*----------------------*/
c={&contrast};

/* Chi-square test */
Xstat=&n.#(c*releff`)*solve(c*v*c`,(c*releff`)`);
df_x=round(trace(ginv(c*v)*c*v));
p_X=1-probchi(Xstat,df_x);
print Xstat p_X;

/* for 2 groups: t-approximation */
%if &ngroups=2 %then %do;
tau=j(1,2,.);
%do i=1 %to &ngroups;
 %let j=%eval(3-&i);
 tau[&i]=(1/&&n&i)*(var(F_&j[loc(group=&i),]) + ((gamma##2)#var(G_&j[loc(group=&i),]))[,+]
        -2#(gamma#covar(F_&j[loc(group=&i),],G_&j[loc(group=&i),]))[,+]);
%do k=1 %to &ncov;
 %do m=(&k+1) %to &ncov;
    tau[&i]=tau[&i]+(2/&&n&i)#gamma[&k]#gamma[&m]#covar(G_&j[loc(group=&i),][,&k],G_&j[loc(group=&i),][,&m]);
  %end;
 %end;
%end;
print tau;
df_t=tau[+]##2/((tau[1]##2)/(&n1-1) + (tau[2]##2)/(&n2-1));
Tstat=sqrt(Xstat);
p_t=2*(1-probt(Tstat,df_t));
print df_t p_t;
/* Add CI*/
* START CHANGE BY KB ON 11JAN06;
%let cialpha=%sysevalf((1-&alpha)*100);
%do i=1 %to &ngroups;
LLn&i=releff&i-probit(1-&alpha./2)*sqrt(variances[&i]/&n.);
ULn&i=releff&i+probit(1-&alpha./2)*sqrt(variances[&i]/&n.);
print "standard error" (sqrt(variances[&i]/&n.));
print "&cialpha% CI using normal approximation" LLn&i ULn&i;
if df_t>0.009 then do; * otherwise overflow error - value is based for 95% CI;
LLt&i=releff&i-tinv(1-&alpha./2,df_t)*sqrt(variances[&i]/&n.);
ULt&i=releff&i+tinv(1-&alpha./2,df_t)*sqrt(variances[&i]/&n.);
print "standard error" (sqrt(variances[&i]/&n.));
print "&cialpha% CI using t-approximation" LLt&i ULt&i;
end;
* END CHANGE BY KB ON 11JAN06;
%end;
%end;

/* Output */
/*--------*/

%if %length(&out)^=0 %then %do;
  create &out var{
%do i=1 %to &ngroups;
releff&i
%end;
%do j=1 %to &ncov;
gamma&j
%end;
Xstat df_X P_X 
%if &ngroups=2 %then %do; 
Tstat df_t P_t
/* Added CI for releffs*/
%do i=1 %to &ngroups;
LLn&i ULn&i LLt&i ULt&i
%end;
%end;
};
  append;
  close &out;
%end;

/* END IML */
quit;

/* Remove temporary datasets */
/*---------------------------*/
proc sql;
drop table _group, _data;
quit;

%mend schacht;