**************************************************************************;
proc iml;
options nodate pagesize=60 linesize=80;

*Study characteristics;
revs = {0, .04, .10, .19, .35, .52};*if revc is 0 then we get problems, instead set to a a small nonzero number;  
stud = {13, 22, 30, 69, 112, 234, 503};
samps = {53, 231, 730};
effect={0, .1, .3, .5};
******************************;
allout = j(12,27,-9); *placeholder for main output;
index = 1;             *output condition index;
********************************************************
*  loops for the study-defining characteristics
********************************************************;
*do l1 = 1 to nrow(revs);
*do l2 = 1 to nrow(stud);
do l3 = 1 to nrow(samps);
do l4 = 1 to nrow(effect);
studies = stud[1, ];
revc = revs[2, ];
nbar=samps[l3, ];
bs=effect[l4, ];
sdn=round(.2#nbar);
sdrho= sqrt(revc);
*print revc studies nbar bs sdn sdrho;
****************************************************************************;
a = .50; *intercept d=.44 median from Brannick survey;
*correlation between moderators;
M= {1 .0 ,
   .0  1}; *corrs between IVs: 0; 
Ry= j(2,1,bs);
   *Ry = {bs,  
      bs}; *corrs of each IV with the DV: 0, .10, .30, .50, we found smallest-largest: .05 .06	.13	.17	.23;
Big = i(3);
big[2,1]= Ry[1,1];
big[1,2]=big[2,1];
big[3,1]=Ry[2,1];
big[1,3]=big[3,1];
big[3,2]=M[2,1];
big[2,3]=big[3,2];
beta = inv(M)*Ry;*beta is the expected value of the vector of regression coefficient relating ES to moderators;
beta = beta`;
E_B1=remove (beta,2);
E_B2=remove (beta,1);
*SDrho = .001; * set the standard deviation of true effect sizes-unconditional sd in d from Brannick survey smallest, 10th, 25th, 50th, 75th, 90th, largest: 0, 0, .20, .32, .44, .59, .72;
*REVC= convert to variance- unconditional variance in d from Brannick survey smallest, 10th, 25th, 50th, 75th, 90th, largest: 0, 0, .04, .10, .19, .35, .52;
Rsq = beta*Ry;*expected value of r_squared;
TarCV = (1-Rsq)#REVC;
******************************************************************************************;
*print M Ry Big beta SDrho REVC Rsq TarCV E_b1 E_b2;
*********************************************************************
choose the number of studies in your meta-analysis;
********************************************************************;
*Nbar = 230;  * smallest-largest: 34 53	105	231	358	730 2425;
*SDN = 47;    * Standard deviation of sample sizes 10, 47;
*studies = 13; *smallest-largest:  13 22	30	69	112	234 503; 
k = studies;
*******************************************************************;
*******************************************************************;
z=10000;  * z is number of iterations;
*******************************************************************;
*******************************************************************;
*print 'input data';
*print 'Corr between moderators' M;
*print 'Corrs between moderators and ES' Ry;
*print 'betas' beta;
*print 'ES variability' SDrho REVC TarCV;
*print 'rho is' M;
*print 'number of studies' studies; 
*placeholders for results, Out1 for individual studies
out2 for meta-analysis;
out1=j(z,2,-9);
out2=j(z,4,-9);
out3=j(z,10,-9);
out4=j(z,2,-9);
out5=j(z,4,-9);


do y=1 to z;
************************************************;
*find data for k studies at population es level;
************************************************;

mods = normal(j(k,3,0));
int = j(k,1,a);               *intercept defined at top;
mods = mods*half(big);
mods[,1]=mods[,1]#SDrho+int;* mods has the infinite sample ES adjusted to have the standard deviation SDrho;
*mods = mods*half(big)#SDrho+int;

*print mods int revc nbar;
*************************************************************;
* draw observed samples from underlying effect sizes
*************************************************************;
*observed = j(k,5,-9);
es= j(k,1,-9);
x=j(k,2,-9);
n_=j(k,1,-9);

do draw = 1 to k;
*PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPP
* pick a sample size;
*******************************************************;
*d1: N=round(normal(0)*SDN+Nbar);
*if N < 10 then goto d1;  *pick the sample size for 1 study;
X1 = RANNOR (0);
X1 = -0.409035830 + 0.905290969*X1 + 0.409035830*X1**2 + -0.023674572*X1**3;*skew=2.2 kurtosis=5.3;
d1: N = round(X1*SDN + Nbar);
if N < 20 then goto d1;  *pick the sample size for 1 study;
************************************************************
*  Pick infinite sample ES
***********************************************************;
ESrho = mods[draw,1];  * godown the rows of MODS to pick ES;
****************************************************
* find observed study by sampling N cases 
*  from infinite ES
****************************************************;
Nhalf =round(N#.5); 
  N= Nhalf#2;
  sam1 = normal(j(Nhalf,1,0));      
  sam2 = normal(j(Nhalf,1,0))+ESrho; 

   m1 = sam1[+,]/Nhalf;     *mean for control group (group1);
   m2 = sam2[+,]/Nhalf;     *mean for experimental group (group2);
   V1 =(sam1-m1)#(sam1-m1); *squared deviations group 1;
   V1 = V1[+,]/(Nhalf-1);   *variance group 1;
   V2 =(sam2-m2)#(sam2-m2); *squared deviations group 2;
   V2 = V2[+,]/(Nhalf-1);   *variance group 2;
   Spool = sqrt(((Nhalf-1)#V1+(Nhalf-1)#V2)/(Nhalf+Nhalf-2)); *pooled SD;
   d = (m2-m1)/Spool;       * raw d;
   d = (1-3/(4#N-9))*d;     *Hedges correction to debias d for small samples;

es[draw,1]=d;
x[draw,1]=mods[draw,2];
x[draw,2]=mods[draw,3];
n_[draw,1]=N;

*print es x n_;

end; *end draw loop;
************************************************
* Compute some preliminaries useful in metas
************************************************;
es2_=es#es;
n2_=n_#n_; 
n_2=n_/2;
nn_=n_2#n_2;
*RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR
* Compute Regression Analyses
*RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR;
****************************************;
*Method of moments;
****************************************;
es_MM = es;
w=(2#nn_#n_)/(2#n2_+nn_#es2_);
w_f=w;****new****;
df = k - 1   ;
j = j(k,1,1) ;
x = j(k,1,1) || x ;
p = j(1,ncol(x),1) ;
v = 1/w ;
				  *Estimate regression parameters and between studies variance;
				  xw = x#(w*p); 
				  xwx = T(xw)*x;
				  xwx_f=xwx;****new****;
                  B = inv(xwx)*T(xw)*es_MM ;
                  qe = sum(es_MM#w#es_MM) - T(B)*xwx*B ;
				  qe_b_fixed=qe;****new****;
                  ww = w#w ; 
				  xww = x#(ww*p) ; 
                  xwwx = T(xww)*x ;
                  MM_ca = (qe-(k-ncol(x)))*inv(sum(w) - trace((xwwx)*inv(xwx))) ;
                  if MM_ca<0 then MM_c = 0;
				  if MM_ca>0 then MM_c= MM_ca;
                  w = 1/(v + MM_c) ;
				  xw = x#(w*p) ; 
                  xwx = T(xw)*x ; *Use between studies variance to re-estimate regression parameters;
				  B = inv(xwx)*T(xw)*es_MM ; *Use between studies variance to re-estimate regression parameters;

		qe = sum(es_MM#w#es_MM) - T(B)*xwx*B ;
		qe_b_random = sum(es_MM#w_f#es_MM) - T(B)*xwx_f*B ;****new****;
        if qe < 0 then qe = 0;
		dfe = nrow(es)-ncol(x) ;
		pe = 1 - probchi(qe,dfe) ;
		alpha=.05;
		***********new***********;
		if qe_b_fixed < 0 then qe_b_fixed = 0;
		pe_b_fixed = 1 - probchi(qe_b_fixed,dfe);

		if qe_b_random < 0 then qe_b_random = 0;
		pe_b_random = 1 - probchi(qe_b_random,dfe);
		*************************;

		*Primary Results;
		if pe < alpha then QE_random_mm_power=1;
		if pe > alpha then QE_random_mm_power=0;
		MM_bias_CV=MM_c-TarCV;
		MM_MSE_CV=(TarCV-MM_c)**2;

		***********new***********;
		if pe_b_fixed < alpha then QE_b_fixed_power=1;
		if pe_b_fixed > alpha then QE_b_fixed_power=0;

		if pe_b_random < alpha then QE_b_random_power=1;
		if pe_b_random > alpha then QE_b_random_power=0;
		*************************;
*print  qe_b_fixed qe_b_random qe pe pe_b_fixed pe_b_random QE_b_fixed_power QE_b_random_power ;

***********************************************************;
*Full information ML = iterative generalized least squares;
***********************************************************;
es_ML = es;
w=(2#nn_#n_)/(2#n2_+nn_#es2_);
df = k - 1   ;
j = j(k,1,1) ;
*x = j(k,1,1) || x ;
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))) ;
                  if ML_c<0 then ML_c = 0 ;
                  w = 1/(v + ML_c) ;
						  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)) ;
                          if ML_c<0 then ML_c = 0 ;
                          w = 1/(v+ML_c) ;
						  end ;

		se_c = sqrt(2/sum(w#w));
		u_ci_se_c= ML_C+1.96*se_c;
		l_ci_se_c= ML_C-1.96*se_c;
		if l_ci_se_c< 0 then l_ci_se_c=0;
		qe = sum(es_ML#w#es_ML) - T(B)*xwx*B ;
		if qe < 0 then qe = 0;
        dfe = nrow(es)-ncol(x) ;
        pe = 1 - probchi(qe,dfe) ;	
		alpha=.05;

		*Primary Results;
		if pe < alpha then QE_random_ml_power=1;
		if pe > alpha then QE_random_ml_power=0;
		if l_ci_se_c = 0 then ML_tau_CI_power=0; 
		if l_ci_se_c > 0 then ML_tau_CI_power=1;
		ML_bias_CV=ML_c-TarCV;
		ML_MSE_CV=((TarCV-ML_c)**2);
		if TarCV >= l_ci_se_c & TarCV <= u_ci_se_c then ML_PC_SE=1;
		else ML_PC_SE=0;
		ML_width=u_ci_se_c-l_ci_se_c;



*print mL_width;

**************************************************;
*Bootstrap;
**************************************************;

w=(2#nn_#n_)/(2#n2_+nn_#es2_);
xa= x[,2];
xb= x[,3];
xc= xa||xb;
*print w es x x1;
full= es || w ||xc ;
*effect=full[,1];
*weight=full[,2];
*mod=full[,3];
k = nrow(full);
*print full;

*********************************************************************;
*  BOOTSTRAP ESTIMATES                  *
*********************************************************************;
*************set up**************************************************;
iter=2000;                     *Set number of iterations for bootstrap;
bootest=j(iter,1,-9); *Set placeholder for output (tau_2);
counter=j(iter,1,-9); 


**********************************************************************;
do az = 1 to iter;  *BEGIN main loop;                                                           
**********************************************************************;
full2 = j(k,ncol(full),-9);
  do i = 1 to k;         	*Start picking study estimates for bootstrap sample;                                                      
    chk:pick = round(uniform(0)*k)+1;
    if pick > k then goto chk;
	full2[i,] = full[pick,];
  end;                      *End i loop picking study estimates for bootstrap sample;  
	  
 *print w2 es2 x2;
***********************************************;
* FIND Tau_2 (do the meta-analysis)    *
***********************************************;
*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:4] ;
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;
				  if BS_ca < MM_ca then count=1;
				  if BS_ca >= MM_ca then count=0;
*print BS_C;


***********************************************
*  record estimate from the meta-analysis     *
***********************************************;
bootest[az,1]=BS_Ca;
counter[az,1]=count;
proportion=counter[+,]/iter;


end;   * END az loop for bootstrap;
*************************************************
*  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;
Boot_width= boot_up-boot_low;


*print boot_width; 

*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;
BC_width=bc_up-bc_low;


***********************************;
*count intervals that need adjusting;
if z_bc_up < 50 then up_adjust_count=1;
if z_bc_up >= 50 then up_adjust_count=0;
if z_bc_low > 1950 then low_adjust_count=1; 
if z_bc_low <= 1950 then low_adjust_count=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;

BC_adj_width= bc_up_adj-bc_low_adj;

*print ML_width Boot_width BC_width BC_adj_width;


************************************;

*Primary Results;
if boot_low = 0 then Boot_CI_power=0; 
if boot_low > 0 then Boot_CI_power=1;
if TarCV >= boot_low & TarCV <= boot_up then Boot_PC_SE=1;
else Boot_PC_SE=0;

if bc_low = 0 then BC_CI_power=0; 
if bc_low > 0 then BC_CI_power=1;
if TarCV >= bc_low & TarCV <= bc_up then BC_PC_SE=1;
else BC_PC_SE=0;

if bc_low_adj = 0 then BC_adj_CI_power=0; 
if bc_low_adj > 0 then BC_adj_CI_power=1;
if TarCV >= bc_low_adj & TarCV <= bc_up_adj then BC_adj_PC_SE=1;
else BC_adj_PC_SE=0;


*create out1 from bootest;
*append from bootest;
*print 'Bootstrap Confidence Interval for the random-effects variance component';
*print MM_C mm_ca /*bootest counter*/ boot_low boot_up z_bc_low z_bc_up bc_low bc_up z_bc_low_adj z_bc_up_adj bc_low_adj bc_up_adj up_adjust_count low_adjust_count  ;




******************storing values;
*Primary Results;

out1[y,1]=MM_bias_CV;
out1[y,2]=ML_bias_CV;

out2[y,1]=ML_PC_SE;
out2[y,2]=Boot_PC_SE;
out2[y,3]=BC_PC_SE;
out2[y,4]=BC_adj_PC_SE;

out3[y,1]=QE_random_mm_power;
out3[y,2]=QE_random_ml_power;
out3[y,3]=ML_tau_CI_power;
out3[y,4]=Boot_CI_power; 
out3[y,5]=BC_CI_power;
out3[y,6]=BC_adj_CI_power;
out3[y,7]=up_adjust_count;
out3[y,8]=low_adjust_count;
out3[y,9]=QE_b_fixed_power;
out3[y,10]=QE_b_random_power;

out4[y,1]=MM_MSE_CV;
out4[y,2]=ML_MSE_CV;

out5[y,1]=ML_width;
out5[y,2]=Boot_width;
out5[y,3]=BC_width;
out5[y,4]=BC_adj_width;

bias= out1[+,]/z;
MSE= out4[+,]/z;
PC= out2[+,]/z;
Power= out3[+,]/z;
Int_width= out5[+,]/z;
*********************************************
*  End Main Study Loop for one condition
*********************************************;
end; * End y loop;
rootMSE=sqrt(MSE);

dat=studies||nbar||rsq||revc||Tarcv||bias||rootMSE||PC||power||Int_width;
allout[index,]=dat;
index = index+1;

*print bias [colname={"MM_bias_CV","ML_bias_CV"}] rootMSE [colname={"MM_MSE_CV","ML_MSE_CV"}] PC [colname={"ML_PC_SE","Boot_PC_SE","BC_PC_SE","BC_adj_PC_SE"}] Power[colname={"QE_random_mm_power","ML_Hedges_Error","ML_tau_CI_Error","Boot_CI_Error","BC_CI_Error","BC_adj_CI_Error"}];
*print dat;*print dat [colname={"nbar" "studies" "rsq" "revc" "Tarcv" "MM_bias_CV" "ML_bias_CV" "MM_MSE_CV" "ML_MSE_CV" "ML_PC_SE" "Boot_PC_SE" "BC_PC_SE""BC_adj_PC_SE" "QE_random_ml_power" "ML_Hedges_Error" "ML_tau_CI_Error" "Boot_CI_Error" "BC_CI_Error" "BC_adj_CI_Error" "up_adjust_count" "low_adjust_count"}];
***********************************
* End condition loops
***********************************;
*end; * l1;
*end; * l2;
end; * l4;
end; * l3;
cname = {"studies" "nbar" "rsq" "revc" "Tarcv" "MM_bias_CV" 
"ML_bias_CV" "MM_MSE_CV" "ML_MSE_CV" "ML_PC_SE" "Boot_PC_SE" 
"BC_PC_SE" "BC_adj_PC_SE" "QE_random_mm_power" "QE_random_ml_power" 
"ML_tau_CI_power" "Boot_CI_power" "BC_CI_power" "BC_adj_CI_power" "up_adjust_count" "low_adjust_count" "QE_b_fixed_power" 
"QE_b_random_power" "ML_width" "Boot_width" "BC_width" "BC_adj_width"};

create allout from allout [colname=cname];
append from allout;
*print allout;
quit;
data new;
set allout;
proc print;
run;


PROC EXPORT DATA= WORK.new 
            OUTFILE= "C:\Documents and Settings\Guy Cafri\Desktop\trial2.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="fd"; 
RUN;