* READ from a permanent SAS dataset;
LIBNAME DS 'c:\SAS\perm\AnhBin';  
DATA mydata; SET DS.AnhBin; RUN;

TITLE ' pooled logistic BY time';
PROC SORT DATA=mydata OUT=mydatasorted; 
BY time; RUN;
PROC logistic DATA=mydatasorted NOPRINT; 
BY time;
MODEL biRadmit (event='1') =NDX NPR LOS DX101 / aggregate scale=none;* lackfit;
OUTPUT OUT=outpool3 P=mu;
RUN;
DATA outpool3; SET outpool3;
wt = mu*(1-mu);
rsdraw = biRadmit-mu; 

PROC SORT DATA=outpool3 OUT=outpool3 ;
  BY PNUM_R time; RUN;

QUIT;
proc iml;
libname reflib "C:\Users\Justin\Documents\My SAS Files\reflib";
RESET STORAGE=reflib.MVIntegration;
LOAD;

use outpool3;                                                            * ## change ####;
read all VARIABLES {NDX NPR LOS DX101} into Zmat;   * ## change ####;
read all var {wt} into wt;
read all var {rsdraw} into rsd;
close outpool3;

corrind = 1; * if corrind =1, then Identity correlation matrix, o/w all elemants =1;
N=1625;T =3; Np = 4;* Np does not include the intercept;
alpha = 0.05;
 
cutoff2 = 0.05; cutoff3 = 0.01;
level1 = 25000;level2 = 250000;level3 = 2500000; 
ABSEPS = 0.000001;
RELEPS = 0;

start rho(a,rsd) global(N,T);
abm = j(N,2*T,.);
abm[,1:T] = shape(rsd,N);	* N x T;
abm[,T+1:2*T] = shape(a,N);
corr = corr(abm);  
rho = corr[1:T,T+1:2*T];    * T x T;
return(rho);
finish rho;

start stddev(a,rsd) global(N,T);
bm = shape(rsd,N);    		 * N x T;
bdev = bm-j(N,1,1)*bm[:,];   * bdev N x T,   bm[:,] Col Mean is a row vector   1 x T;
bdev2 = bdev#bdev;      
am = shape(a,N);   
adev = am-j(N,1,1)*am[:,];  
adev2 = adev#adev;      
stddev = sqrt( (1/N)*t(bdev2)*adev2 );   * T x T;
return(stddev);
finish stddev;

* corrected standardization;
start stdzn(x) global(N,T);
xrows = shape(x,N);   *N x T;
y = xrows - xrows[:,];  *N x T;
vcv = (1/(N-1))*t(y)*y; * T x T;
v = diag(vcv);
sinv = sqrt(inv(v));
x2 = y*sinv;   *N xT;
x2  = shape(x2,N*T,1);
return(x2);
finish stdzn;


pvec = j(Np*T*T,1,.);sevec = j(Np*T*T,1,.);  * pvec   (T*T) x Np;
r4out = j(T,T*Np,.); se4out = j(T,T*Np,.); z4out = j(T,T*Np,.); p4out = j(T,T*Np,.); 

y = rsd;
y_std = stdzn(y);
DO i=1 TO Np;
x = wt#Zmat[,i]; 			 * (N*T) x 1;
x_std = stdzn(x);
r = rho(x_std,y_std);		 * T x T;
se = stddev(x_std,y_std);	 * T x T;
z = sqrt(N)*(r/se);
p = 2*(1-cdf('normal',abs(z)));  * T x T;

r4out[,T*(i-1)+1:T*i] = r;
se4out[,T*(i-1)+1:T*i] = se;
z4out[,T*(i-1)+1:T*i] = z;
p4out[,T*(i-1)+1:T*i] = p;

DO j = 1 TO T;
p[j,j] = 1;  * Not going to test diagonal elements;
END;
pvec[T*T*(i-1)+1:T*T*i,1] = shape(p,T*T,1);    *(T*T) x 1;
sevec[T*T*(i-1)+1:T*T*i,1] = shape(se,T*T,1);   *(T*T) x 1;
END;

pse = j(Np*T*T,2,.); pse[,1] = pvec; pse[,2] = sevec; 
call sort(pse,{1});  * sort by 1st column in ascending order (by default);

stop = 0;
rnking = 1;

DO WHILE (stop<1);

pmin = pse[rnking,1];

se4test = pse[rnking:Np*T*T,2]; * row vector;

L = Np*T*T - rnking +1;


PRINT rnking;
PRINT pmin;

if pmin = 0 then DO; pact = 0 ; print 'Integration Skipped'; print pact; END;
if pmin >=0.5 then DO; pact = 1; print 'Integration Skipped'; print pact; END;

if pmin>0 & pmin <0.5 then DO;

 LOWER = probit(pmin*0.5)*J(1,L,1);
 UPPER = probit(1-pmin*0.5)*J(1,L,1);
 INFIN = J(1,L,-1);
 Corr1 = I(L); Corr2 = j(L,L,1); if corrind = 1 then Corr = Corr1; else Corr = Coor2;
 COVAR = diag(se4test)*Corr*diag(se4test); * variance matrix for test;

/* MAXPTS = level1;*/
/* RUN MVN_DIST( L, LOWER, UPPER, INFIN, COVAR, MAXPTS, ABSEPS, RELEPS,  ERROR, VALUE, NEVALS, INFORM );*/
/* PRINT ERROR VALUE NEVALS INFORM;*/
/* pact = 1- VALUE;*/
/* print 'Level1';  print pact;*/
/* if pact<cutoff2 then DO;*/
/*	MAXPTS = level2;*/
/*	RUN MVN_DIST( L, LOWER, UPPER, INFIN, COVAR, MAXPTS, ABSEPS, RELEPS,  ERROR, VALUE, NEVALS, INFORM );*/
/*	pact = 1- VALUE;*/
/*	print 'Level 2 Integreation'; print pact;*/
/* END;*/
/* if pact<cutoff3 then DO;*/
		MAXPTS = level3;
		RUN MVN_DIST( L, LOWER, UPPER, INFIN, COVAR, MAXPTS, ABSEPS, RELEPS,  ERROR, VALUE, NEVALS, INFORM );
	pact = 1- VALUE;
	print 'Level 3 Integreation';	print pact;
/* END; */

END;   * to close    if pmin>0 & pmin <0.5 then DO;

 if pact >= alpha | L=1 then stop = 1;
 if pact < alpha & L>1 then DO;	rnking = rnking +1; END;

END;  * to close   DO WHILE (stop<1);


print 'final pvalue';
print pact, L;


* by Multiple Test;
TypeVec = (pvec >= pmin*j(Np*T*T,1,1) );
Type = shape(TypeVec, Np, T*T);

TypeMtx = j(Np+T,T*T,.);
Tint = I(T);
TypeMtx[1,]= shape(Tint,1,T*T);
TypeMtx[2:1+Np,] = Type;
Tt2 = j(T,T,0);Tt2[,2] = Tint[,2];
Tt3 = j(T,T,0);Tt3[,3] = Tint[,3];
TypeMtx[Np+2,] = shape(Tt2,1,T*T);
TypeMtx[Np+3,] = shape(Tt3,1,T*T);

CREATE TypeMtx from TypeMtx;
APPEND from TypeMtx;
CLOSE TypeMtx;
print 'Each row of TypeMtx is the type vector for each of the predictors, by multiple test';
print 'But the intercept admitd only T valid equations';
print TypeMtx;


* by Individual Tests;
TypeVec2 = (pvec >= alpha*j(Np*T*T,1,1) );
Type2 = shape(TypeVec2, Np, T*T);
TypeMtx2 = TypeMtx;
TypeMtx2[2:1+Np,] = Type2;
CREATE TypeMtx2 from TypeMtx2;
APPEND from TypeMtx2;
CLOSE TypeMtx2;

print 'Each row of TypeMtx2 is the type vector for each of the predictors, by individual tests';
print 'But the intercept admitd only T valid equations';

print TypeMtx2;


Type4out = j(T,T*(Np+T),.);
DO i=1 to Np+T;
Type4out[,T*(i-1)+1:T*i] = shape(TypeMtx[i,],T,T);
END;


print '===================================================================================';
print 'Note: r4out, se4out, z4out, p4out consider only non-intercept,non-time-dummy variables';
print '===================================================================================';
print r4out;
print se4out;
print z4out;
print p4out;
print 'Type4out is TypeMtx rearranged, for Excel output, each TxT block is the type vector for a predictor';
print 'i.e. by multiple test and the intercept admitd only T valid equations';
print Type4out; 

x = wt; 			 * (N*T) x 1;
x_std = stdzn(x);
r = rho(x_std,y_std);		 * T x T;
se = stddev(x_std,y_std);	 * T x T;
z = sqrt(N)*(r/se);
p = 2*(1-cdf('normal',abs(z)));  * T x T;
print r, se, z, p;
T_wt = p>0.05;
T_wtT2 = j(T,T,0);
T_wtT2[,2] = T_wt[,2];
T_wtT3 = j(T,T,0);
T_wtT3[,3] = T_wt[,3];
print T_wt;
T_wt = shape(T_wt,1); 
T_wtT2 = shape(T_wtT2,1); 
T_wtT3 = shape(T_wtT3,1); 
print T_wt, T_wtT2, T_wtT3;
TypeMtx3 = j(Np+T,T*T,.);
TypeMtx3[1,] = T_wt;
TypeMtx3[2:Np+1,] = Type2;
TypeMtx3[Np+2,] = T_wtT2;
TypeMtx3[Np+3,] = T_wtT3;

CREATE TypeMtx3 from TypeMtx3;
APPEND from TypeMtx3;
CLOSE TypeMtx3;
print 'Each row of TypeMtx3 is the type vector for each of the predictors, by individual tests';
print 'But the intercept and time dummies may admit more than T valid equations';
print TypeMtx3;
Quit;