**************************************************
*  A monte carlo program to sample studies from 
*  existing meta analyses taken from JAP, PPsych
*  and Academy of Mgmt Journal.
**************************************************;
* grab data on meta analyses from Excel
**************************************************;
PROC IMPORT OUT= WORK.meta 
            DATAFILE= "C:\Documents and Settings\pcd4100\Desktop\all_metas_11.17.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="Sheet1$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
*output average n skew and kurtosis by study;
* output datset is nbar;
proc means data= work.meta mean kurtosis skew nway;
class MA;
var N;
output out=nbar
mean=m_nbar
kurtosis=k_nbar
skew=s_nbar;

* get the mega distribution of the statistics from each M-A;
*proc univariate data=nbar plot;
*var m_nbar k_nbar s_nbar;

*get medians;
*proc means data=nbar median nway;
*var m_nbar k_nbar s_nbar;
*we want to group the meta-analyses into four groups based on average sample size and on skew of sample size;
*median splits;
Data conditions;
set nbar;
if m_nbar le 168.5689655 then nbardi=0;
if m_nbar gt 168.5689655 then nbardi=1;
*proc print;
*look at the design;
*proc freq data=conditions;
*tables skewdi*nbardi;
*Run;
*make groups;
data conditionsrevised;
set conditions;
if nbardi=0 then group=1;
if nbardi=1 then group=2;
keep ma group;
*********************************************************
*  Data on delta statistics (rather than rho)
*  for updated meta
*********************************************************
Ave(d)_L Ave(d)_M Ave(d)_H 
0.20 0.45 0.98 
SDd_L SDd_M SDd_H 
0.04 0.24 0.50 
*********************************************************;
*proc print;
*descriptives by group;
*proc means data=conditionsrevised;
*class group;
*var m_nbar k_nbar s_nbar;
*we want each meta-analysis represented by a row of sample sizes;
*dropping variables and transposing;
data dropped;
set work.meta;
keep MA n;
proc sort data=dropped;
by MA;
proc transpose data=dropped out=drop2;
by MA;
run;
*this has produced a matrix with MA number and then the number of sample sizes
* contained in the study, one study per column;
* print;
*merge all info;
data all;
merge drop2 conditionsrevised;
by MA;
*proc print data=all;
Proc sort data=all;
by group;
data final;
set all;
keep col1-col97 group;
*proc print;

data group1;
set final;
if group=1;
*proc print;
data group2;
set final;
if group=2;
*proc print;
run;

proc iml;
****************************************************************************************************
******************************START IML*************************************************************
****************************************************************************************************;
* read data into a matrix for each group;
use group1 var _ALL_;
read all var _ALL_ into x1;
numbmetas1= nrow(x1);
numbcol1= ncol(x1)-1;
read all var _ALL_ into x2;
numbmetas2= nrow(x2);
numbcol2= ncol(x2)-1;
*print x1 x2 x3 x4;
*print numbmetas1 numbcol1 numbmetas2 numbcol2 numbmetas3 numbcol3 numbmetas4 numbcol4;
* set up bootstrap estimates;
**************************************************;
*set population values for Delta and SDdelta;
delta = .20; *high infinite sample ES;
SDdelta = .04; *high infinite sample SD for ES;
**************************************************;
*Input the number of iterations or simulated meta-analyses, i;
i=5000;
print delta SDdelta i;
out10 = j(i,8,-9); *placeholder for results of a meta-analysis;
Y = j(i,16,-9); * Y is the placeholder for all the outputs from all the groups
**************************************************************************************
*  1. first, we have to pick the group we want to sample from.
*  2. then, we have to pick the study from the group.
*  3. then we have to find the number of samples in the study.
*  4. then we have to sample that number of samples of the given sizes for that study
**************************************************************************************;
*  1.  Pick the group;
**************************************************************************************;
p=2;
do group = 1 to p;* p is the number of iterations and group is its counter;
*  group is the main loop for different groups;
if group = 1 then x = x1;
if group = 2 then x = x2;
***********************************
**************************************;
* start to randomly pick meta from the chosen group;
**************************************;
do l = 1 to i;  * i is the number of iterations and l is its counter;
*  l is the main loop;
*print x;
**********************************************************************************;
*  2.  Choose the study from the group;
**********************************************************************************;
rowcount = nrow(x); * find the number of rows in X;
*print rowcount;
study = round(uniform(0)*rowcount+.5); *find a random number from 1 to the number of studies;
if study < 1 then study = 1;
if study > rowcount then study = rowcount;  *limit the values of study just in case;
*print group study;
***********************************************************************************;
* 3.  find the number of samples in the study
***********************************************************************************;
samples = x[study,];
*print samples;
Ns = 0;
do a = 1 to ncol(x) until (samples[,a]<4);
Ns = a-1; if Ns = 98 then Ns = 97;
end;
*print Ns;
holder1 = j(Ns,3,-9);  *placeholder for the d values for this meta-analysis;
***********************************************************************************
*   Ns is the number of studies in the chosen meta analysis
***********************************************************************************;
*  4.  Sample the number of studies of the given sample sizes from the population
*      of interest;
***********************************************************************************;
sizes = samples[,1:Ns]; *pick the columns with nonzero sample sizes;
sizes = sizes`; *transpose it to be a column vector;
*print sizes;
****************************************************************************************
* Now we can begin the main loop, by sampling data for a meta-analysis
****************************************************************************************;
  *Generate Data for a Meta-Analysis;
  *****************************************************;
  **Set the sample size to be the first column in the sample1 matrix
  *****************************************************;
k = Ns;       *k is the number of studies for the meta-analysis;
do h=1 to k;  *this loop generates data for 1 meta;
  ********************************************************************;
  * Sample delta from a distribution;
  * Note the distribution was set up above from the listed parameters;
  ********************************************************************;
 N = sizes[h,1];   *N is the sample size for the kth study;
 if mod(N,2)=1 then N=N+1; *Make N an even number so subgroup size is exactly half N;
 InfES=RANNOR(0)#SDdelta+delta;       *InfES is the infinite sample Effect Size for the kth study;
   Nhalf =round(N#.5);                *Number of people per group, assuming equal groups;
   sam1 = normal(j(Nhalf,1,0));       *data for control group;
   sam2 = normal(j(Nhalf,1,0))+InfES; *data for experimental group;
  ***********************************
  *  compute the sample d           *
  ***********************************;
   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;
   *print d;
   *****************************************************************;
  holder1[h,1]=d;       *d for the analysis;
  holder1[h,2]=N;       *N for the analysis;
  holder1[h,3]=InfES;   *Infinite sample d (local delta) to check on the program functioning;
end; *end the k loop to create data for one meta analysis;
******************************************************************************
* Holder1 has the data for 1 meta analysis and Infinite sample ES as check
******************************************************************************;
************************************************
* Compute Mean and SD for Infinte ES as a check;
************************************************;
MInfES = holder1[+,3]/k;
SDinfES = (holder1[,3]-MInfES)#(holder1[,3]-MInfES);
SDinfES = sqrt(SDinfES[+,]/(k-1));
out10[l,1]=MInfES;    *mean of infinite sample d values;
out10[l,2]=SDinfES;   *SD of infinite sample d values;
*print MInfES SDinfES;
*******************************************************************************
*  Compute meta analyses for 1 meta
*******************************************************************************;
***********************************
**Schmidt and Hunter Bare Bones***;
***********************************;
out1=holder1;
es = out1[,1];
w = out1[,2];
nbar=w[+,]/k;
*print k nbar;
SHdelta=(w`*es)/w[+,];
bigd=j(nrow(out1),1,SHdelta);
vard=w`*((es-bigd)#(es-bigd))/w[+,];
vare=((nbar-1)/(nbar-3))#(4/nbar)#(1+(SHdelta#SHdelta/8));
SHresd=vard-vare;
if SHresd < 0 then SHresd=0;
SHsddelta=sqrt(SHresd);
out10[l,3]=SHdelta;
out10[l,4]=SHsddelta;
*print SHdelta SHsddelta vard vare;
****************************************
Unit weight model
****************************************;
unitvec=j(k,1,1);
unaved = (unitvec`*es)/k;
bigu = j(nrow(out1),1,unaved);
varun = unitvec`*((es-bigu)#(es-bigu))/k;
sampdu=((nbar-1)/(nbar-3))#(4/nbar)#(1+(unaved#unaved/8));
unresd = varun - sampdu;
if unresd < 0 then unresd = 0;
unsddelta = unresd**.5;
out10[l,5] = unaved;
out10[l,6] = unsddelta;
*print unaved unsddelta varun sampdu;
*******************************************************************
*RE Model (Hedges & Vevea), with d as the ES           *
*Computations;                        *
*******************************************************************;
* First find initial inverse variance weights
*******************************************************************;
 do b=1 to k; * compute weights for each study using sample size and d value;
  halfN = holder1[b,2]/2; *N1 = N2 = .5*N;
  d = holder1[b,1];       *the corrected d value;
  w[b,1] = (2#halfN#halfN#(halfN+halfN))/((2#(halfN+halfN)#(halfN+halfN))+halfN#halfN#d#d); *the weight;
 end; * end computation of Hedges initial inverse variance weights;
df = k - 1;
 mes = sum(es#w)/sum(w);
 q = sum((es#es)#w) - sum(es#w)#sum(es#w)/sum(w);
 c    = (q - df)/(w[+,]-sum(w#w)/w[+,]) ;
 if c<0 then c = 0 ;
 HVSDdelt = sqrt(c); 
 ***********************************************************;
 * find the random effects weights;
 ***********************************************************;
wre    = 1/(1/w + c) ;
 ******************************************
 * Find the RE mean ES
 ******************************************;
mesre = sum(es#wre)/sum(wre) ;
out10[l,7]= mesre;
out10[l,8]= HVSDdelt;
*print mesre HVSDdelt;
end; * end main i loop;
*print out10;
*********************************************************
* Find average and RMS across iterations
*********************************************************;
popmean = j(i, 1, delta);
poptau = j(i, 1, SDdelta);
avdelta = out10[+,1]/i;
avSDdelt = out10[+,2]/i;
*schmidt hunter;
rmsrSHM = (out10[,3]-delta)#(out10[,3]-delta);
rmsrSHM = sqrt(sum(rmsrSHM)/i);
rmsrSHT = (out10[,4]-SDdelta#SDdelta)#(out10[,4]-SDdelta#SDdelta);
rmsrSHT = sqrt(sum(rmsrSHT)/i);
wtmeanSH = sum(out10[,3])/i;
wtTauSH = sum(out10[,4])/i;

*hedges & Vevea;
rmsrHVM = (out10[,7]-delta)#(out10[,7]-delta);
rmsrHVM = sqrt(sum(rmsrHVM)/i);
rmsrHVT = (out10[,8]-SDdelta#SDdelta)#(out10[,8]-SDdelta#SDdelta);
rmsrHVT = sqrt(sum(rmsrHVT)/i);
wtmeanHV = sum(out10[,7])/i;
wtTauHV = sum(out10[,8])/i;

*unit weight;
rmsrunM = (out10[,5]-delta)#(out10[,5]-delta);
rmsrunM = sqrt(sum(rmsrunM)/i);
rmsrunT = (out10[,6]-SDdelta#SDdelta)#(out10[,6]-SDdelta#SDdelta);
rmsrunT = sqrt(sum(rmsrunT)/i);
wtmeanun = sum(out10[,5])/i;
wtTauun = sum(out10[,6])/i;

*****************************************************************;
print avdelta AvSDdelt wtmeanSH wtmeanHV wtmeanun wtTauSH wtTauHV wtTauun
rmsrSHM rmsrHVM rmsrunM rmsrSHT rmsrHVT rmsrunT;
*print out10;
Y[,((group-1)*8+1)] = out10[,1];
Y[,((group-1)*8+2)] = out10[,2];
Y[,((group-1)*8+3)] = out10[,3];
Y[,((group-1)*8+4)] = out10[,4];
Y[,((group-1)*8+5)] = out10[,5];
Y[,((group-1)*8+6)] = out10[,6];
Y[,((group-1)*8+7)] = out10[,7];
Y[,((group-1)*8+8)] = out10[,8];
end; * end the p loop for the 2 groups;
*print Y;
create final from Y;
append from Y;
*print final;
quit;
data last;
set final;
MInfES_ln=col1;
SDInfES_ln=col2;
SHd_ln = col3;
SHtaud_ln = col4;
UwMd_ln = col5;
Uwtaud_ln = col6;
HVMd_ln = col7;
HVtaud_ln = col8;
MInfES_hn=col9;
SDInfES_hn=col10;
SHd_hn = col11;
SHtaud_hn = col12;
UwMd_hn = col13;
Uwtaud_hn = col14;
HVMd_hn = col15;
HVtaud_hn = col16;
Keep MInfES_ln SDInfES_ln SHd_ln SHtaud_ln UwMd_ln Uwtaud_ln HVMd_ln HVtaud_ln MInfES_hn SDInfES_hn SHd_hn SHtaud_hn UwMd_hn Uwtaud_hn HVMd_hn HVtaud_hn;
file 'C:\Documents and Settings\pcd4100\Desktop\StatMos (.20,.04).dat';
put (MInfES_ln SDInfES_ln SHd_ln SHtaud_ln UwMd_ln Uwtaud_ln HVMd_ln HVtaud_ln MInfES_hn SDInfES_hn SHd_hn SHtaud_hn UwMd_hn Uwtaud_hn HVMd_hn HVtaud_hn) (7.3);
run;
proc univariate data=last plot;
var MInfES_ln SDInfES_ln SHd_ln SHtaud_ln UwMd_ln Uwtaud_ln HVMd_ln HVtaud_ln MInfES_hn SDInfES_hn SHd_hn SHtaud_hn UwMd_hn Uwtaud_hn HVMd_hn HVtaud_hn;
proc means data=last;
var MInfES_ln SDInfES_ln SHd_ln SHtaud_ln UwMd_ln Uwtaud_ln HVMd_ln HVtaud_ln MInfES_hn SDInfES_hn SHd_hn SHtaud_hn UwMd_hn Uwtaud_hn HVMd_hn HVtaud_hn;
Run;
*proc univariate data=last plot;
*proc means data=last;
*var SHrho_g1 SHresr_g1 rM_HVr_g1 TauHVr_g1 rM_HVz_g1 TauHVz_g1 SHrho_g2 SHresr_g2 rM_HVr_g2 TauHVr_g2 rM_HVz_g2 TauHVz_g2 
SHrho_g3 SHresr_g3 rM_HVr_g3 TauHVr_g3 rM_HVz_g3 TauHVz_g3 SHrho_g4 SHresr_g4 rM_HVr_g4 TauHVr_g4 rM_HVz_g4 TauHVz_g4;
*proc plot data=last; 
*plot col2*col4 col7*col9 col12*col14 col17*col19;
*data l1;
*set last;
*type = 1;
*y = col1;
*keep type y;
*data l2;
*set last;
*type = 2;
*y = col5;
*keep type y;
*data all;
*set l1 l2;
*proc shewhart;
*boxchart y*type;
*run;