/*Create Data Set in SAS
   VD= Number Vaccinated With Disease
   VWD= Number Vaccinated Without Disease
   NVD = Number Not Vaccinated With Disease
   NVWD = Number Not Vaccinated Without Disease
*/
OPTIONS FORMCHAR="|----|+|---+=|-/\<>*";

data BCG;
   input TRIAL VD VWD NVD NVWD LATITUDE YEAR ALLOC;
   LN_OR = log((vd/vwd)/(nvd/nvwd));
   EST   = 1/vd + 1/vwd +  1/nvd + 1/nvwd;
cards;
1	4	  119		11	   128	44	48	1
2	6	  300		29	   274	55	49	1
3	3	  228		11	   209	42	60	1
4	62 	13536	     248	 12619  52	77	1
5	33	 5036	      47	  5761  13	73	2
6	180	 1361	     372	  1079  44	53	2
7	8	 2537	      10	   619	19	73	1
8	505	87886	     499	 87892	13	80	1
9	29	 7470	      45	  7232	27	68	1
10	17	 1699	      65	  1600	42	61	3
11	186	50448	     141	 27197	18	74	3
12	5	2493	       3	  2338	33	69	3
13	27	16886	      29	 17825	33	76	3
;

proc print;run;

/*Mantel-Haenszel method for combining odds ratios*/
data MH;
  set BCG;
  treat=1; response="Disease"; count=vd;   output;
  treat=1; response="NoDis";   count=vwd;  output;
  treat=2; response="Disease"; count=nvd;  output;
  treat=2; response="NoDis";   count=nvwd; output;
  keep trial treat response count;run;

proc print;run;

proc freq data=mh;
weight count; 
tables trial*treat*response/chisq relrisk   cmh;
run;

/*Vote-Counting*/
data pvalues;
input trial pvalue positive;
cards;
1   .1030 1
2   .0001 1
3   .0235 1
4   .0001 1
5   .3354 1
6   .0001 1
7   .0001 1
8   .8494 2
9   .0461 1
10  .0001 1
11  .0022 1
12  .5378 2
13  .9483 1
;

proc freq data=pvalues;
  tables positive / binomial(p=.5);
  exact binomial;
run;

data fisher;
  set pvalues;
  P=-2*log(pvalue);
run;
proc means data=fisher Sum;
var P;
output out=new sum(P) = TotalP;
run;

options nocenter;
data fisher;
  TotalP = 126.4;
  df = 2*13;
  FisherP = 1-probchi(TotalP,df);
run;
proc print data=fisher2;
var TotalP df FisherP;
format FisherP 20.19;
run;

/*Fixed Effects Analysis, Combined Log Odds Ratio*/
ods output fitstatistics = fit1;
proc mixed data=bcg;
class trial;
model ln_or= / s cl;
repeated / group=trial;
parms  / parmsdata=bcg eqcons=1 to 13;
run;

/*Random Effects Analysis, Combined Log Odds Ratio*/
ods output fitstatistics = fit2;
proc mixed data=bcg;
class trial;
model ln_or= / s cl;
random int / subject=trial s;
repeated / group=trial r;
parms (.5) (0.35712)(0.20813)
           (0.43341)(0.02031)
           (0.05195)(0.00991)
           (0.22701)(0.00401)
           (0.05698)(0.07542)
           (0.01253)(0.53416)
           (0.07164) /eqcons=2 to 14;run;


/*Meta Regression using Latitude as a Predictor*/
proc mixed data=bcg itdetails;
class trial;
model ln_or= latitude/ s cl ddf=1000;
random int / subject=trial s;
repeated / group=trial;
parms (.5) 
           (0.35712)(0.20813)
           (0.43341)(0.02031)
           (0.05195)(0.00991)
           (0.22701)(0.00401)
           (0.05698)(0.07542)
           (0.01253)(0.53416)
           (0.07164) /eqcons=2 to 14;
run;

/*Check for Publication Bias*/

data BCG;
input TRIAL VD VWD NVD NVWD LATITUDE YEAR ALLOC;
LN_OR = log((vd/vwd)/(nvd/nvwd));
est   = 1/vd + 1/vwd +  1/nvd + 1/nvwd;
wi = 1/est;
wt_ln_or = wi*ln_or;
group=1;
n=sum(vd,vwd,nvd,nvwd);
root_n = sqrt(n);
cards;
1	4	  119		11	   128	44	48	1
2	6	  300		29	   274	55	49	1
3	3	  228		11	   209	42	60	1
4	62 	13536	     248	 12619      52	77	1
5	33	 5036	      47	  5761      13	73	2
6	180	 1361	     372	  1079      44	53	2
7	8	 2537	      10	   619	19	73	1
8	505	87886	     499	 87892	13	80	1
9	29	 7470	      45	  7232	27	68	1
10	17	 1699	      65	  1600	42	61	3
11	186	50448	     141	 27197	18	74	3
12	5	2493	       3	  2338	33	69	3
13	27	16886	      29	 17825	33	76	3
;

proc means data=bcg noprint;
  by group;
  output out=meandat sum(wt_ln_or)=wt_sum_ln_or
                     sum(wi) = sum_wt;
run;
  
proc print data=meandat;
run;
data combine;
 merge bcg meandat;
 by group;

 tbar = wt_sum_ln_or / sum_wt;
 wt_sq_ln_or_diff = wi*(ln_or-tbar)**2;
 vi_star = est - (1/sum_wt);
 Ti_star = (ln_or - tbar)/(sqrt(est - (1/sum_wt)));
 Zi = (ln_or - Tbar)/sqrt(est);
 se = sqrt(est);
run;

proc print data=combine;
run;
proc univariate plot;
  var zi;
  histogram;
  qqplot / normal(mu=est sigma=est);
run;

/*Examples of Tests for Publication Bias*/

options nodate pageno=1;
proc reg data=combine;
  model zi = se;
  title "Publication Bias Regression Model";
  title2 "Based on Precision";
run;
proc reg data=combine;
  model zi = root_n;
  title "Publication Bias Regression Model";
  title2 "Based on Square Root of Sample Size";
run;
proc corr data=combine kendall;
  var vi_star ti_star;
  title "Publication Bias Kendall's Tau";
  title2 "Between Vi* and Ti*";
run;