/*--------------------------------------------------------------------*
 *                                                                    *
 *  SAS/IML PROGRAM FOR COMPUTING TICV INDEX, IMPACT MATRIX AND       *
 *  ABSORPTION MATRIX                                                 *
 *  BASED ON DR. KEN FRANK'S FORMULA                                  *
 *                                                                    *
 *  MODULE WRITTEN BY LIXIONG GU                                      *
 *                                                                    *
 *  MODULE INPUT:                                                     *
 *      CORRELATIN MATRIX: NEWCORR2                                   *
 *      SPECIFY DEPENDENT VARIABLE NAME IN %LET NAMEY=                *
 *      SPECIFY TOTAL NUMBER OF INDEPENDENT VARIABLES IN %LET NIV=    *
 *      SPECIFY ALL VARIABLE NAMES IN %LET ALLNAME=                   *
 *      SPECIFY DATA FORMAT IN %LET DATAFOR=                          *
 *            0 -- SAS DATASET (DEFAULT)                              *
 *            1 -- EXCEL 2000/20002 FILE                              *
 *      SPECIFY DATA TYPE IN %LET DATATYPE=                           *
 *            0 -- CORRELATION MATRIX (DEFAULT)                       *
 *            1 -- RAW DATA                                           *
 *                                                                    *
 *                                                                    *
 *  MODULE OUTPUT:                                                    *
 *    beta SEbeta observed t                                          *
 *    TICV IEVR(pi=.5) IEVR(TAUxy=0) OIRC                             *
 *                                                                    *
 *--------------------------------------------------------------------*/
proc means data=in1.g9;
run;

proc corr data=in1.g9;
var mloipe interven;
partial pregpa pvt lginc pared6 female;
run;

proc reg data=in1.g9;
model mloipe = interven pregpa pvt lginc pared6 female;
run;

proc reg data=in1.g9;
model mloipe = interven pregpa pvt lginc pared6 female casp1;
run;

proc corr data=in1.g9;
var casp1 interven;
partial pregpa pvt lginc pared6 female;
run;

proc corr data=in1.g9;
var casp1 mloipe; 
partial pregpa pvt lginc pared6 female;
run;

options nodate pageno=1 ls=200 nocenter;
libname in1 "D:\AddHealth\bboard\AERA_wkshop\29feb_final";
data in1.g9;
set in1.grade911_grp1;
*if w1gr11_flg=1 then outgpa=ogpa3;
if w1gr10_flg=1 then outgpa=ogpa3;
if w1gr9_flg=1 then outgpa=ogpa2;

if w1gr10_flg=1 then pregpa=ogpa2;
if w1gr9_flg=1 then pregpa=ogpa1;

if w1gr9_flg=1;
missmlope=0;
if mlope=. then do;
missmlope=1;
mlope=0;
end;

if mttrack ge 3 then interven=1;
if 0 le mttrack < 3 then interven=0;

*outcome=outgpa;
outcome=mloipe;



/*  Specify data name for calculating, and alpha level for t test  */
%LET alpha=.05;

*Model 1.1 ;
/* Specify DEPENDENT variable */
%LET pre_namey  =  mloipe;
%LET namey=%upcase(&pre_namey);

/* Specify all variables to be included in the computation */
%LET allname=mloipe interven pregpa pvt lginc pared6 female;


/* Specify total number of independent varaibles */
%LET niv=6;

/* Specify file data name */
%LET filename = in1.g9;

/* Specify data format:
    0 -- SAS DATASET (DEFAULT);
    1 -- EXCEL 2000/20002 FILE  */
%LET datafor = 0;

/* SPECIFY DATA TYPE
    0 -- CORRELATION MATRIX (DEFAULT)
    1 -- RAW DATA           */
%LET datatype = 1;

/*
*Model 2.1 ;
%let namey=REALABOR;
%let allname=REALABOR WAGECORD WDINT2 DUMY9397 UNEMPLOY INFLATIO TRADE UNIONDEN LAGLABOR;
*Model 1.2;
%let namey=UNEMPLOY;
%let allname=UNEMPLOY WAGECORD WDINT1 ULINT2 ULINT3 UNION_PA LEFTIST DUMY9297 TRADE UNIONDEN EMPLOYRG  TAX_RATE BENIFITD LAGGED_U;
*/
%MACRO RTOC(FILENAME, DATAFOR, DATATYPE);
 %IF &DATAFOR = 1 %THEN %DO;
   PROC IMPORT OUT= WORK.rawdata
            DATAFILE= "&FILENAME"
            DBMS=EXCEL2000 REPLACE;
     GETNAMES=YES;
   RUN;
  %IF &DATATYPE = 1 %THEN %DO;
    proc corr data = rawdata out=newcorr2;
     var &allname;
    run;
   %END;
  %ELSE %DO;
    DATA newcorr2;
     set rawdata;
    run;
   %END;
 %END;
 %ELSE %IF &DATATYPE = 1 %THEN %DO;
    proc corr data = &filename out=newcorr2;
     var &allname;
    run;
   %END;
  %ELSE %DO;
    data newcorr2;
     set &filename;
    run;
  %END;
%MEND;

%RTOC(&FILENAME,&DATAFOR,&DATATYPE);

* Gett sample size;
data usen; set newcorr2;
if _type_="N";
ny=&namey.;
keep ny;
run;

data usename; set newcorr2;
 if _name_ ^= " ";
 keep _name_;
run;

data usestd; set newcorr2;
 if _type_="STD";
 drop _type_ _name_;
run;

proc transpose data=usestd out=usestd1 (drop=_name_ rename=(COL1=STD));
run;

data coronly;
 set newcorr2;
 where _type_ = "CORR";
 drop _type_ _name_;
run;

/**********************************************************************************
*                                                                                 *
*                        Begin IML session                                        *
*                                                                                 *
***********************************************************************************/
proc iml;

*Start TICV module;
start TICV (TICVx,TICVy,x_ind,sdyx,sx,ny,opi,b1,sebeta,tp,rp,risq,tresult,tname,subimpact);
   title 'TICV';
   ***********************
   *     Compute ryg     *
   ***********************;
   nx=ncol(TICVx);

   TICVg1=TICVx[2:nx,2:nx];
   TICVy1=TICVy[2:nx];
   subrii=inv(TICVg1);
   Byg=subrii*TICVy1;
   ryg=sqrt(sum(Byg#TICVy1));
   subpr=i(nx-1);
   do u=1 to nx-1;
     do w=1 to nx-1;
       if w ^= u then do;
         subpr[u,w]=-subrii[u,w]/sqrt(subrii[u,u]#subrii[w,w]);
       end;
     end;
   end;

   do subi=1 to nx-1;
    subb1=Byg[subi]#sdyx[1]/sdyx[subi+1];
    subseb=(sdyx[1]/sdyx[subi+1])#sqrt((1-ryg#ryg)/(ny-nx-2))#sqrt(subrii[subi,subi]);
    subtp=subb1/subseb;
    subrp=subtp/sqrt(ny-nx-2+subtp#subtp);
    subrxypg=subrxypg||j(nx-1,1,subrp);
    subrxypg[subi,subi]=TICVy[subi+1];
   end;

*   print subpr;
*   print subrxypg;

   subimpact=subpr#subrxypg;

  /*******************************************************
   *   compute t critical value and r#                   *
   *******************************************************/
   q=nx+1;
   df=ny-q;
   tpct=&alpha;
   *tpct=&alpha/2;
   tcritical = tinv(1-tpct, df);
   if rp < 0 then tcritical=-tcritical;
   rpound = tcritical / sqrt(df + tcritical # tcritical);

/********************
 CHECK THIS
*********************/
   rpp=tcritical/sqrt((ny-q)+tcritical**2);

  /*******************************************************
   *             Compute TICV index                      *
   ******************************************************/
   m=sqrt((1-risq)#(1-ryg#ryg));
   if rp > 0 then signrp=1;
   if rp < 0 then signrp=-1;
   if rp =0 then signrp=0;
   TICV = signrp # m # ((Abs(rp)-abs(rpound))/(1-abs(rpound)));

   rycvsq = abs(TICV) # sqrt((1 - ryg # ryg) / (1 - risq));
   rxcvsq = abs(TICV) # sqrt((1 - risq) / (1 - ryg # ryg));

   IEVRPI5=2 # rpound - rp;
   IEVRr0 = 1 - (rpound / rp);
*   print rpound rp rpp tp ny nx IEVRPI5 IEVRr0;

   star1=(ny#rp)#(ny#rp)+4#tcritical#tcritical#(tcritical#tcritical-q);
   if star1<0 then do;
    nstar=0;
    opi=0;
   end;
   else do;
    nstar=ny#(ny#rp#rp-2#tcritical#tcritical+abs(rp)#sqrt(star1))/(2#tcritical#tcritical);   IEVCr0=ny*rp/(ny+round(nstar,1));
    opi=nstar/(ny+nstar);
   end;

   IEVCpi5=2 # rpp - rp;
   setpi=.5;
   q=3+nx+1;
   m1=1.96*sqrt(1/((1-setpi)#(ny-q))+1/(setpi#(ny-q)));
   e2m=exp(2*m1);
   rxylb=(rpound-e2m#(1+rpound)-1-sqrt(4#(e2m-1)#(setpi-1)#(e2m#setpi+rpound+e2m#rpound-setpi)+(1-rpound+e2m#(1+rpound))**2))/(2#(e2m-1)#(setpi-1));
   OIRC=(b1-sebeta # tcritical) / b1;
   sypx=sdyx[1]#sqrt(1-risq);
   tsbn=tp#sypx/(b1#sqrt(ny));
   TVAI1star=1-tsbn;
   TVAI2star=1-4#tsbn/3;
   TVAI3star=1-.94#tsbn;
   ESR=1-.2#sypx/b1;
   lbb1=b1-tcritical#sebeta;
   TVAI0=sx#sqrt(1-risq)#(lbb1/b1-1)+1;
   tresult = b1 || sebeta || tp || TICV || rycvsq || rxcvsq || IEVRPI5 || IEVRr0 || OIRC || nstar || IEVCPI5 || rxylb || TVAI1star || TVAI2star || TVAI3star || ESR || TVAI0;
   tname = {"beta" "SEbeta" "observed t" "TICV" "rycv square" "rxcv square" "IEVR(pi=.5)" "IEVR(rxy=0)" "OIRC" "n*" "IEVC(pi=.5)" "rxylb" "TVAI*(Beta)" "TVAI**(Beta)" "TVAI***(Beta)" "ESR(Beta,.2)" "TVAI*(beta0)"};
finish;


/* read data for computing TICV */

* Read variable names;
namey={&namey.};
varname=t({&allname.});

*Get the position of y variable in the complete correlation matrix;
nvar=nrow(varname);
do pos=1 to nvar;
 if varname[pos]="&namey" then y_pos=pos;
end;
*print varname nvar namey;
print varname nvar Y_POS;
*Read standard deviation into matrix;
use usestd1;
read all var ({std}) into stdall;

*Read correlation matrix with all variables included;
*Transform the whole correlation matrix so that the first variable would be y;
use coronly;
read all var ({&allname.}) into allmat;

vname=varname;
if y_pos ^=1 then do;
  sw=i(nvar);
  sw[1,1]=0;
  sw[1,y_pos]=1;
  sw[y_pos,1]=1;
  sw[y_pos,y_pos]=0;
  allv=sw * allmat * sw;
  stdall=sw * stdall;
  vname[1]=varname[y_pos];
  vname[y_pos]=varname[1];
end;
else allv=allmat;
  namex=vname[2:nvar];
  tnamex=t(namex);

print stdall;

/*************************************************
* Make matrix x and y so that                    *
*  y={ryx1,ryx2,...,ryxn}.                       *
*  x={1,rx1x2,...rx1xn,                          *
*     rx2x1,1,...rx2xn,                          *
*     ...                                        *
*     rxnx1,rxnx2,...1};                         *
**************************************************/
y=allv[2:nvar,1];
x=allv[2:nvar,2:nvar];
print y;
print x;
*Read sample size into matrix n;
use usen;
read all var ({ny}) into n;

*Get number of independent variables;
nx=ncol(x);

*Get parameter estimates for the full regression model;
invrx=inv(x);
b=invrx*y;
rsquared=sum(b#y);
print b rsquared;
/*****************************************************
*    rii is the diagonal elements of the inverse     *
*    of correlation matrix X. rii=1/(1-Ri_Squared)   *
******************************************************/
rii=invrx;
pr=i(nx);
do u=1 to nx;
 do w=1 to nx;
  if w ^= u then do;
    pr[u,w]=-rii[u,w]/sqrt(rii[u,u]#rii[w,w]);
  end;
 end;
end;

* print pr [rowname=namex colname=tnamex format=8.3];

result_init =1;
subimp=j((nx-1)*nx,nx-1,0);

do x_index=1 to nx;
  stdyx=stdall[1];
  stdyx=stdyx//stdall[x_index+1];
  if x_index<nx then
    sdyx=stdall[1:x_index]//stdall[x_index+2:nx+1];
  else sdyx=stdall[1:x_index];
*  print stdall stdyx sdyx;
  sw=i(nx);
  rname = namex;
  if x_index >1 then do;
    sw[1,1]=0;
    sw[1,x_index]=1;
    sw[x_index,1]=1;
    sw[x_index,x_index]=0;
    rname[1]=namex[x_index];
    rname[x_index]=namex[1];
  end;

  TICVx=sw * x * sw;
  TICVy=sw * y;
  b1=b[x_index]#stdyx[1]/stdyx[2];  /* parameter estimate   */
  risq=1-1/rii[x_index,x_index];
  sebeta=(stdyx[1]/stdyx[2])#sqrt((1-rsquared)/(n-nx-1))#sqrt(rii[x_index,x_index]);
  tp=b1/sebeta;
  rp = tp / sqrt(n-nx-1 + tp # tp);  /*convert observed t-ratio to a correlation for each predictor*/
  rxypg=rxypg||j(nx,1,rp);
  rxypg[x_index,x_index]=TICVy[1];
*  print rxypg TICVy;
  run TICV(TICVx,TICVy,x_index,sdyx,stdyx[2],n,0.5,b1,sebeta,tp,rp,risq,result,cname,subimpact);
*  print subimpact;
*  print rname;
  do w=1 to nx-1;
   do u=1 to nx-1;
     wi=(nx-1)*(x_index-1)+w;
     subimp[wi,u]=subimpact[w,u];
   end;
  end;
  if result_init=1 then do;
      ticvout=result;
      result_init=0;
    end;
  else ticvout = t(t(ticvout) || t(result));
end;

* print subimp;
* print rxypg [rowname=namex colname=tnamex format=8.3];

/***************************
*     Impact Matrix        *
***************************/
impact=pr#rxypg;

/***************************
*     Absorption Matrix    *
***************************/
impall=j((nx-1)*nx,nx-1,0);
do x_index=1 to nx;
  sw=i(nx);
  rname = namex;
  if x_index >1 then do;
  sw[1,1]=0;
  sw[1,x_index]=1;
  sw[x_index,1]=1;
  sw[x_index,x_index]=0;
  rname[1]=namex[x_index];
  rname[x_index]=namex[1];
  end;

  impact1=sw * impact * sw;
  woname=rname[1];
  rname=rname[2:nx];
  tcname=t(rname);
  impact1=impact1[2:nx,2:nx];
  do w=1 to nx-1;
   do u=1 to nx-1;
     wi=(nx-1)*(x_index-1)+w;
     impall[wi,u]=impact1[w,u];
   end;
  end;
  if x_index =1 then do;
    rrname=rname;
    wwoname=woname;
  end;
  else do;
    rrname=rrname//rname;
    wwoname = wwoname//woname;
  end;
end;

  woxi=impall/subimp;

* print woname "absorb colname with respect to rowname",
  woxi [rowname=rname colname=tcname format=8.3];
  subimp1=rrname||char(woxi,8,3);
*  subname={name}||t(rrname);
*  create subimpdata from subimp1[colname=subname];
  create subimpdata from subimp1;
  append from subimp1;

  create wond from wwoname;
  append from wwoname;

*print impact [rowname=namex colname=tnamex format=8.3];
*print ticvout [rowname=namex colname=cname format=8.3];
imp1=namex||char(impact,8,3);
tn1={name}||tnamex;
create impdata from imp1[colname=tn1];
append from imp1;
ticv1=namex||char(ticvout,8,3);
cn1={name}||cname;
create ticvdata from ticv1[colname=cn1];
append from ticv1;

quit;

*title "Absorption Matrix: how much x1 absorbs x2 with respect to x3";
*proc print data=subimpdata noobs;
*run;
title "Impact Matrix: impact of column on row, correlation with outcome on diagonal";
proc print data=impdata noobs;
run;
title "TICV Indices";
proc print data=ticvdata noobs;
run;

%macro abmat(dataf=subimpdata, namef=wwoname, nvar=%eval(&niv-1));

%let nv = %eval(&nvar+1);

proc sql noprint;
  select col1 into :v1-
  :v&nv
    from work.wond;
  quit;
run;

%do nmat = 1 %to &nv;
 data subimp&nmat;
   set &dataf (firstobs=%eval((&nmat-1)*&nvar+1) obs=%eval(&nmat*&nvar));
 run;

proc sql noprint;
  select col1 into :var2-
  :var&nv
    from work.subimp&nmat;
  quit;
run;

proc datasets library=work;
  modify subimp&nmat;
  rename col1=var_name
    %do i=2 %to &nv;
             col&i=&&var&i
    %end;
  ;
  quit;

title "Absorption Matrix &&v&nmat.";
title2 "How much predictor &&v&nmat. absorbs the impact of colname on rowname";
proc print data=subimp&nmat noobs;
run;

%end;

%mend;

*%abmat();