Data McLeod;
Input author$ year N age Parent$ TechP$ InformP$ Dx$ Tf$ TechD$ InformD$ rxy;
es =.5*log((1+rxy)/(1-rxy));
w = n-3;
info = 3;
if informp eq 'C' then info = 1;
if informp eq 'P' then info = 2;
cards;
Asarnow 1994	56	10.43	M	I 	P 	Y	C	S	C	0.26
Asanrow 2001	156	12.31	M	I 	P 	Y	C	S	C	0
Barber 1996	158	12	B	O	O	N	C	O	P	0.16
Baron 1989	144	15.2	B	Q	C	N	C	S	C	0.19
BARRERA 1992	94	14.6	B	Q	C	N	C	S	C	0
BIGGAM 1998	125	18.8	B	Q	C	N	CU	S	C	0.18
BRAGE 1993	156	14	B	Q	C	N	CU	S	C	0.22
CORONA 2005	111	13.2	M	O	O	N	C	S	C	0.16
ELDER 1992	76	12	B	O	O	N	C	SO	CO	0.28
FIELD 1987	38	5.1	M	QO	PO	N	C	O	O	0.1
FIELD 1995	455	16.6	B	Q	C	N	C	S	C	0.27
FOREHAND 1988	89	13.08	M	O	O	N	C	S	C	0.23
FOREHAND 1988	69	13.42	M	Q	C	N	C	SO	CO	0.18
FURUKAWA 1992	165	17.47	B	Q	C	N	U	S	C	0
FURUKAWA 1997	144	17.47	B	Q	C	N	U	S	C	0.15
GALAMBAS 1990	91	11.58	M	Q	CP	N	U	S	C	0.26
GALLIMORE 1992	35	13.63	B	Q	CP	N	C	S	C	0.22
GARBER 1997	240	11.86	M	Q	CP	N	C	SO	CP	0.19
GARBER 2001	240	11.86	M	Q	CP	N	U	S	C	0.11
GREENBERGER 1996	84	13.1	B	Q	C	N	C	S	C	0.56
GREENBERGER 1996	89	13.2	B	Q	C	N	C	S	C 	0.53
HAROLD 1997	146	12.83	B	Q	C	N	C	S	C	0.6
HAROLD 1997	380	13	B	Q	C	N	C	S	C	0.32
HEAVEN 2004	276	15.34	B	Q	CP	N	C	S	C	0.18
HUNTLEY 1990	76	7.78	B	Q	P	N	C	S	C	0.15
JACQUEZ 2004	72	15	M	Q	CP	N	C	SO	CP	0.31
KOBAK 1991	48	15.7	M	O	O	N	U	S	C	0.22
MARMORSTEIN 2004	249	17.5	B	Q	C	Y	L	S	C	0.31
MARTIN 1994	681	15	B	Q	C	N	L	S	C	0.3
MCCLELLAN 2004	493	15.5	M	Q	C	N	C	S	C	0.12
MCFARLANE 1995	801	17.1	B	Q	C	N	U	S	C	0.28
MESSER 1995	20	10.1	B	Q	CP	N	C	S	C	0.29
MILNE 2001	59	15.7	M	Q	C	N	U	S	C	0.41
PUIG-ANTICH 1985	92	9.43	B	Q	P	Y	C	S	C	0.54
PUIG-ANTICH 1993	93	15.2	B	Q	P	Y	C	S	C	0.67
ROGERS 2003	306	11.7	B	Q	C	N	C	S	C	0.34
RUDOLPH 1997	81	9.65	M	Q	C	N	C	S	C	0.37
SANDERS 1992	46	10.4	M	QO	CPO	Y	C	S	C	0.26
SHEEBER 1998	52	15.5	M	Q	CP	Y	C	S	C	0.42
SHEK 1989	2150	16	B	Q	C	N	CU	S	C	0.26
STEIN 2000	68	10.5	B	Q	C	Y	C	S	C	0.37
STOCKER 1994	85	7.92	M	Q	C	N	CU	S	C	0.25
TESSER 1991	147	13	B	Q	C	N	C	S	C	0.11
THOMPSON 1999	54	13.9	M	Q	C	N	U	S	C	0.39
WHITBECK 1992	451	12	B	QO	CPO	N	U	SO	CP	0.2
*proc print;
************************************************
*The macro that runs the regression begins here;
************************************************;
%MACRO metareg(es,w,x,dsn=_last_) ;
*********************************************************
* Much of this program was written by David Wilson;
* Parts also written by Guy Cafri and Michael Brannick;
**********************************************************;
      proc iml;
        use &dsn ;
        read all var {&es} into es  where(&es^=. & &w^=.) ;
        read all var {&w}  into w   where(&es^=. & &w^=.) ;
        read all var {&x}  into x   where(&es^=. & &w^=.) ;
*************************************************************
		You must have created a SAS dataset that contains the 
		effect sizes (ES), the inverse variance weights (w)
		and the moderator (x).  This statement reads them 
		all into IML, and deletes any cases missing on any
		of the variables.
*********************************************************************;
*********************************************************************
		First, compute the unconditional REVC,
		The random-effects variance component without
		considering moderators.  This is our estimate of the amount 
		of 'true' variablility we have to examine and 
		account for using moderators.
*********************************************************************;
 w_orig = w;
 x_orig = x;
 k = nrow(es)  ;                                           *k effect sizes (studies);
 dft = k - 1   ;                                           *total degrees of freedom;
 mes = sum(es#w)/sum(w);                                   *fixed weighted mean effect size;
 qt = sum((es#es)#w) - sum(es#w)#sum(es#w)/sum(w);         *total q;
 pt   = 1-probchi(qt,dft) ;                                *p-value for total q;
 c    = (qt - dft)/(w[+,]-sum(w#w)/w[+,]) ;                *unconditional REVC (random-effects variance component);
 if c<0 then ; do ; c = 0 ; end ;                          *set to zero if negative;
print k [label="number of studies for regression analysis"];
print mes    [label="Fixed effects Weighted Mean"];   
print c [label="unconditional random-effects variance component"];
print 'chi-square test of H0: unconditional REVC = 0 is: ' qt pt;
*RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR
* Compute Regression Analyses
*RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR;
****************************************;
*Method of moments, fixed regression;
****************************************;
es_MM = es;              *to keep the name of ES unique for each analysis;
df = k - 1   ;           * redefine total df;
dfe = nrow(es)-ncol(x) ; * regression error degrees of freedom;
j = j(k,1,1) ;           * row vector of 1;
x = j(k,1,1) || x ;      * augmented design (IV) matrix to include intercept;
p = j(1,ncol(x),1) ;     * column vector of 1;
v = 1/w ;                * estimation (sampling) variances;

				  *Estimate regression parameters and between studies variance;
				  xw = x#(w*p); 
				  xwx = T(xw)*x;
                  B = inv(xwx)*T(xw)*es_MM ;              *fixed regression weights;
                  qe = sum(es_MM#w#es_MM) - T(B)*xwx*B ;  *weighted sum of squared residuals;
                  qr = qt-qe;                             *Analog to regression SS;
                  
				  ww = w#w ; 
				  xww = x#(ww*p) ; 
                  xwwx = T(xww)*x ;
                  MM_ca = (qe-(k-ncol(x)))*inv(sum(w) - trace((xwwx)*inv(xwx))) ; *conditional variance;
                  if MM_ca<0 then MM_c = 0;
				  if MM_ca>0 then MM_c= MM_ca;            *MM_c is the conditional REVC at end of fixed reg;
				  dfr = ncol(x)-1;
				  se = sqrt(vecdiag(inv(xwx))) ;
                  zvalues = B#vecdiag(inv(diag(se))) ;
                  pvalues = (1 - probnorm(abs(zvalues)))*2 ;
                  lowerB = B - se*1.96 ;
                  upperB = B + se*1.96 ;
				  RSQ1 = (qt-qe)/qt;                  *RSQ analog;
				  RSQ_MMf = (c-MM_c)/c;         	 *RSQ change from unconditional to fixed conditional;
				  if RSQ_MMf < 0 then RSQ_MMf = 0;
                  pe = 1 - probchi(qe,dfe) ;              *Chi-square for test of Conditional REVC;
				  pr = 1 - probchi(qr,dfr);               *Chi-square for fixed regression SS;
************************************************
*  Results for MM Fixed
************************************************;	
Print / '------Fixed-effects Weighted Regression-----';

                    
              print  "------- Homogeneity Analysis  -------";
                qtab1 = qr // qe ;
                qtab2 = dfr // dfe ;
                qtab3 = pr // pe  ;
                print qtab1 [rowname={"Model","Residual"}]
                            [colname="Q"][label="Source"] [format=12.5]
                            qtab2 [label="df"]
                            qtab3 [label="p"] [format=12.5];
Print 'Conditional REVC (after removing moderator) is:' MM_c;

              print "------- Regression Coefficients -------" ;
              names = {&x} ;
              names = T(names) ;
              group = {"Constant"} // names ;
              print b [rowname=group] [colname="B"] [label="Variable"]
                          [format=10.5]
                        se [label="SE"] [format=10.5]
                        lowerB [label='-95%CI'] [format=10.5]
                        upperB [label='+95%CI'] [format=10.5]
                        zvalues  [label="z"] [format=10.5]
                        pvalues  [label="p"] [format=10.5]  ;
Print '---------R-Square Analogs-------------------';
             RSQS = RSQ1//RSQ_MMf;
             print RSQS [rowname={"Model Q/Total Q","Pct Reduction in REVC"}][colname="RSQ"]; 
****************************************;
*Method of moments, Random regression;
****************************************;
Print / '-----------Random Weighted Regression, Method of Moments---------------';
                  w = w_orig;
                  xw = x#(w*p) ; xwx = T(xw)*x ;                                   *run fixed regression;
                  B = inv(xwx)*T(xw)*es ;
                  qe = sum(es#w#es) - T(B)*xwx*B ;
                  ww = w#w ; xww = x#(ww*p) ; xwwx = T(xww)*x ;
                  cr = (qe-(k-ncol(x)))*inv(sum(w) - trace((xwwx)*inv(xwx))) ;      *random CRevc;
                  if cr<0 then cr = 0 ;
                  wr = 1/(1/w + cr) ;                            *Random Weights;
        xw = x#(wr*p) ; xwx = T(xw)*x ;                          *Recompute intermediates;
        B = inv(xwx)*T(xw)*es ;                                 *Random regression coefficients;
        meanes = sum(es#wr)/sum(wr) ;                             *Random conditional wts mean;
        q =  sum((meanes - es)#(meanes - es)#wr) ;               *Random q total (not meaningful);
        qe = sum(es#wr#es) - T(B)*xwx*B ;                        *Random q residual (not meaningful);
        qr = q - qe ;                                           *Random q for model;
        dfr = ncol(x)-1 ;                                       *Df for regression model;
        pr = 1 - probchi(qr,dfr) ;                              *p value for regression model;
        se = sqrt(vecdiag(inv(xwx))) ;                          *standard errors;
        zvalues = B#vecdiag(inv(diag(se))) ;
        pvalues = (1 - probnorm(abs(zvalues)))*2 ;
        lowerB = B - se*1.96 ;
        upperB = B + se*1.96 ;
        r2 = qr/(qr+qe) ;                                       *RSQ analog;
		RSQ_MMr = (c-cr)/c;          	                        *RSQ change from unconditional to random conditional;
        if RSQ_MMr < 0 then RSQ_MMr = 0;                        *set to zero if negative;
        RSQS = r2//RSQ_MMr;                                     *concatonate;
*********************************************************************;
*  BOOTSTRAP ESTIMATES                  *
*********************************************************************;
**************************************************;
* The bootstrap is used to estimate a Confidence Interval
* for the conditional REVC in method of moments;
**************************************************;

*************set up**************************************************;
iter=2000;                     *Set number of iterations for bootstrap (leave this alone);
bootest=j(iter,1,-9);          *Set placeholder for output (tau_2);
counter=j(iter,1,-9); 
full= es || w_orig ||x[,2];            *Full contains the matrix of all data from which to sample;
**********************************************************************;
do az = 1 to iter;  *BEGIN main bootstrap loop;                                                           
**********************************************************************;
full2 = j(k,ncol(full),-9); * Placeholder for studies selected at random with replacement;
***********************************************;
* This loop picks k studies at random with replacement from original studies;
  do i = 1 to k;         	* Start picking study estimates for bootstrap sample;                                                      
    chk:pick = round(uniform(0)*k)+1;  *pick a random integer between 1 and k;
    if pick > k then goto chk;         *if the integer is greater than k, try again;
	full2[i,] = full[pick,];           *place the row from full into the store in full2;
  end;                      *End i loop picking study estimates for bootstrap sample;  
*********************************************************;
* FIND Tau_2 (do the meta-analysis)on the new sample    *
*********************************************************;
*Method of moments;
es_BS= full2[,1];
w=full2[,2];
df = k - 1   ;
j = j(k,1,1) ;
x = j(k,1,1) || full2[,3] ;
p = j(1,ncol(x),1) ;

				  *Estimate regression parameters and between studies variance;
				  xw = x#(w*p);
				  xwx = T(xw)*x ;
                  B = inv(xwx)*T(xw)*es_BS ;
                  qe = sum(es_BS#w#es_BS) - T(B)*xwx*B ;
                  ww = w#w ; xww = x#(ww*p) ; 
				  xwwx = T(xww)*x ;
                  BS_ca = (qe-(k-ncol(x)))*inv(sum(w) - trace((xwwx)*inv(xwx))) ;
                  if BS_ca<0 then BS_c = 0;
				  if BS_ca>=0 then BS_c = BS_ca;     *BS_c is the conditional REVC for the sample;
				  if BS_ca < MM_ca then count=1;     *Count is used to adjust the bootstrap for bias;
				  if BS_ca >= MM_ca then count=0;
***********************************************
*  record estimate from the meta-analysis     *
***********************************************;
bootest[az,1]=BS_Ca;
counter[az,1]=count;
end;   * END az loop for bootstrap;
proportion=counter[+,]/iter;
*************************************************
*  Bootstrap is finished
*************************************************;
rankvar = rank(bootest[,1]);
boot_low = bootest[loc(rankvar=50),1];
boot_up = bootest[loc(rankvar=1950),1];
if boot_low<0 then boot_low=0;
if boot_up<0 then boot_up=0;
*Calculate bias correction (BC estimates);
zup=1.96;
zlow=-1.96;
z0=probit(proportion);
z_bc_low= z0+(z0+zlow);
z_bc_up= z0+(z0+zup);

*********************************;
*Makes Lower bound 1 instead of 0;
quant1=probnorm(z_bc_low)#iter;
if quant1<.5 then quant1=.5;
z_bc_low=round(quant1);

quant2=probnorm(z_bc_up)#iter;
if quant2<.5 then quant2=.5;
z_bc_up=round(quant2);
********************************;

bc_low = bootest[loc(rankvar=z_bc_low),1];
bc_up = bootest[loc(rankvar=z_bc_up),1];
if bc_low<0 then bc_low=0;
if bc_up<0 then bc_up=0;

*Adjustment to bias corrected interval;
if z_bc_up < 50 then z_bc_up_adj=50;
if z_bc_up >= 50 then z_bc_up_adj=z_bc_up;
if z_bc_low > 1950 then z_bc_low_adj=1950;
if z_bc_low <= 1950 then z_bc_low_adj=z_bc_low;


bc_low_adj = bootest[loc(rankvar=z_bc_low_adj),1];
bc_up_adj = bootest[loc(rankvar=z_bc_up_adj),1];
if bc_low_adj<0 then bc_low_adj=0;
if bc_up_adj<0 then bc_up_adj=0;
bootsl = j(2,1,-9);
bootsu = bootsl;
bootsl[1,1] = boot_low;
bootsu[1,] = boot_up;
bootsl[2,1] = bc_low;
bootsu[2,1] = bc_up;
bootlab = {'raw bootstrap','bias corrected bootstrap'};
************************************************
*  Results for MM Random
************************************************;	
print 'Chi-square test of model';
print qr [label = 'Model Q'] pr [label ='p-value'];

print cr [label= ' '] [rowname = 'Conditional REVC (moderator removed) = ']
                                [format=12.8];
              print bootsl [rowname=bootlab][label ="Conditional REVC Confidence Intervals"]
              [colname='-95_CI'][format = 10.5]
                    bootsu [label ='+95_CI'][format = 10.5];
print "------- Regression Coefficients -------" ;
              names = {&x} ;
              names = T(names) ;
              group = {"Constant"} // names ;
              print b [rowname=group] [colname="B"] [label="Variable"]
                          [format=10.5]
                        se [label="SE"] [format=10.5]
                        lowerB [label='-95_CI'] [format=10.5]
                        upperB [label='+95_CI'] [format=10.5]
                        zvalues  [label="z"] [format=10.5]
                        pvalues  [label="p"] [format=10.5]  ;
			Print '---------R-Square Analogs-------------------';
             print RSQS [rowname={"Model Q/Total Q","Pct Reduction in REVC"}][colname="RSQ"]; 

***********************************************************;
***********************************************************;
*Full information ML = iterative generalized least squares;
* (Maximum likelihood random)
***********************************************************;
			 Print / '-----------Random Weighted Regression, Maximum Likelihood---------------';

es_ML = es;
w=w_orig;
df = k - 1   ;
j = j(k,1,1) ;
x = j(k,1,1) || x_orig ;
p = j(1,ncol(x),1) ;
v = 1/w ;
				  xw = x#(w*p) ; 
                  xwx = T(xw)*x ;
                  B = inv(xwx)*T(xw)*es_ML ;
                  qe = sum(es_ML#w#es_ML) - T(B)*xwx*B ;
                  ww = w#w ; 
                  xww = x#(ww*p) ; 
                  xwwx = T(xww)*x ;
                  ML_c = (qe-(k-ncol(x)))*inv(sum(w) - trace((xwwx)*inv(xwx))) ; *staring conditional REVC;
                  if ML_c<0 then ML_c = 0 ;
                  w = 1/(v + ML_c) ;                                             *iterated solution;
						  do i=1 to 100;
                          xw = x#(w*p) ; 
                          xwx = T(xw)*x ;
                          B = inv(xwx)*T(xw)*es_ML ;
                          r = es_ML - x*B ;
                          ML_c = sum((w##2)#(r##2-v))/
                                (sum(w##2)) ;                                     *final conditional REVC;
                          if ML_c<0 then ML_c = 0 ;
                          w = 1/(v+ML_c) ;
						  end ;

		se_c = sqrt(2/sum(w#w));                      *Standard error for Conditional REVC;
		u_ci_se_c= ML_C+1.96*se_c;                    *upper bound estimate for REVC;
		l_ci_se_c= ML_C-1.96*se_c;                    *lower bound estimate for REVC;
		if l_ci_se_c< 0 then l_ci_se_c=0;

        meanes = sum(es#w)/sum(w) ;                           *Random conditional wts mean;
        q =  sum((meanes - es)#(meanes - es)#w) ;              *Random q total (not meaningful);
        qe = sum(es#w#es) - T(B)*xwx*B ;                       *Random q residual (not meaningful);
        qr = q - qe ;                                           *Random q for model;
        dfr = ncol(x)-1 ;                                       *Df for regression model;
        pr = 1 - probchi(qr,dfr) ;                              *p value for regression model;
        se = sqrt(vecdiag(inv(xwx))) ;                          *standard errors;
        zvalues = B#vecdiag(inv(diag(se))) ;
        pvalues = (1 - probnorm(abs(zvalues)))*2 ;
        lowerB = B - se*1.96 ;
        upperB = B + se*1.96 ;
        r2 = qr/(qr+qe) ;                                       *RSQ analog;
		RSQ_ML = (c-ML_c)/c;          	                        *RSQ change from unconditional to random conditional;
        if RSQ_ML < 0 then RSQ_ML = 0;                        *set to zero if negative;
        RSQS = r2//RSQ_ML;                                     *concatonate;
************************************************
*  Results for Maximum Likelihood Random
************************************************;
print 'Chi-square test of model';
print qr [label = 'Model Q'] pr [label ='p-value'];
print ML_c [label= ' '] [rowname = 'Conditional REVC (moderator removed) = ']
                                [format=12.8];
              print l_ci_se_c [label ='-95_CI'][format = 10.5]
                    u_ci_se_c [label ='+95_CI'][format = 10.5];
print "------- Regression Coefficients -------" ;
              names = {&x} ;
              names = T(names) ;
              group = {"Constant"} // names ;
              print b [rowname=group] [colname="B"] [label="Variable"]
                          [format=10.5]
                        se [label="SE"] [format=10.5]
                        lowerB [label='-95_CI'] [format=10.5]
                        upperB [label='+95_CI'] [format=10.5]
                        zvalues  [label="z"] [format=10.5]
                        pvalues  [label="p"] [format=10.5]  ;
			Print '---------R-Square Analogs-------------------';
             print RSQS [rowname={"Model Q/Total Q","Pct Reduction in REVC"}][colname="RSQ"]; 

quit;
%MEND metareg ;
%metareg(es,w,age,dsn=McLeod);
run;