OPTIONS FORMCHAR="|----|+|---+=|-/\<>*";

ods rtf file='P:\meta_macros\sgs.rtf' style=Journal;
/*
*************************************************************
ods statement is optional.

This ods statement creates a rich text file of the output 
and places it in the folder meta_macros on the P-disk drive.
You will need to change this for your own computing 
environment.  

The program runs properly even if you do not have a valid 
drive name. SAS does give you an error warning, but it can 
be ignored. (Try running it as is.)

The output file name is   sgs.rtf in the folder P:\meta_macros.

Note that repeated use will overwrite this, so open the output 
file in word, and save it under another name when done.

Note: The SAS log will complain about Division by Zero in some
applications. These are not errors.

***********************************************************************
Reference for MACRO:  

Shuster JJ, Guo JD, Skyler JS. Meta-analysis of safety for 
low event-rate binomial trials. Res Synth Methods. 2012 
Mar;3(1). doi: 10.1002/jrsm.1039. PubMed PMID: 24339834; 
PubMed Central PMCID: PMC3856441.

***********************************************************************

*/

%macro sgs (ddset,ft,nt,fc,nc);
/*

***********************************************************************
This Macro will produce point and 95% confidence limits 
for relative risk based, on the Shuster-Guo-Skylar random effects 
meta-analysis method for low event-rate two-sample binomial comparisons.  
Note that all studies are weighted in proportion to the total sample size 
for the trial, and the method does not exclude studies where zero failures 
occur on one or both treatments. (No continuity corrections are made and 
they are strongly discouraged.) 

We shall apply this to the Nissen-Wolski 2007 NEJM paper (referenced below)
for rosiglitazone risk of Myocardial Infarction in Type II Diabetes.  We 
retrieved data for the six studies with no events on both arms, raising 
the number of studies from 42 to 48 (NW discarded 6 for that reason) 

Nissen SE, Wolski K.  Effect of rosiglitazone on the risk of myocardial 
infarction and death from cardiovascular causes.  New England Journal
of Medicine 2007; 356(24): 2457-2471.
***********************************************************************

The input for this Macro, called as below, is defined as follows:

%sgs (tab3,mi1,n1,mi2,n2);

ddset=tab3 tells the macro the dataset to be analyzed is tab3, 
a SAS dataset. This contains data on each study to be 
analyzed.  You can use any dataset name instead of tab3.

ft=mi1  The treatment failure number in dataset tab3 is mi1
(Myocardial infarctions)

nt=n1  The treatment sample size in dataset tab3 is n1

fc=mi2  The control failure number in dataset tab3 is mi2

nc=n2  The control sample size in dataset tab3 is n2

***********************************************************************
Output:

Printouts of (1) the raw data, and 
(2) study by study point and interval estimates for the 
individual study relative risk (risk ratio) (done by 
natual logs then exponentiating the result back to 
the original scale [natural antilog]) are provided.

The tabulated summary data provide the  number of studies combined, 
the point estimate for summary relative risk, the 95% confidence 
interval for the summary relative risk , and the two-sided P-value.  

***********************************************************************
Notes:

Instead of a normal approximation, we use a more conservative but much 
more accurate t-approximation with M-2 degrees of freedom,(M=number 
of studies) based on our massive simulation study of nearly 40,000 
scenarios for low event-rates and low number of trials being combined 
(5-20), each replicated 100,000 times.

We do not provide forest plots, Cochran's Q, or I-squared.   The study 
specific relative risk and their error properties are highly 
inaccurate in low event-rate situations (not just no event 
rate situations). 

Our methods are based on the ratio of two summary proportions,to
create an overall proportion for each treatment, and then take 
the ratio of these. Study-specific proportions can be estimated 
without bias even when event rates are low. See our reference 
for details.  Logs of risk ratios and odds ratios  have poor low event 
rate approximations both bias and assessment of sampling error.
***********************************************************************
*/
data a;set &ddset;
ft=&ft;nt=&nt;fc=&fc;nc=&nc;
	 rr=(ft/nt)/(fc/nc);
	 logrr=log(rr);
	 pt=ft/nt;pc=fc/nc;xxx=1;
	 if nt=. and ft=0 then pt=0;
	 if nc=. and fc=0 then pc=0;
	 se=((1-pt)/(pt*nt))+ ((1-pc)/(pc*nc));
	 se=sqrt(se);
	 ll=exp(logrr-1.96*se);
      ul=exp(logrr+1.96*se);
	  k=_n_;y=k-30*int(k/29.99);
	 file print;
	 if y=1 then do;
	 put _page_;
	 put 'se=' se;
	 put  // @25 'Printout of Raw Data'; 
	 put // @10 'Study' @22 'Treatment_Fail/N'  @44 'Control_Fail/N'  //;
     end; 
	 put @10 k @22 ft  ' / ' nt   @44 
fc  ' / ' nc ;
data h;set a;n1=nt;n2=nc;f1=ft;f2=fc;study=_n_;
keep study f1 n1 f2 n2;
	data a;set a;
		 file print;
	 if y=1 then do;
	 put _page_;
	 put  // @25 'Printout of Study-specific Sumary Relative Risk Data'; 
	 put // @10 'Study' @22 'Est Relative Risk'  @44 '95% Confidence Limits' //;
	 end;
	     put @10 k @22 rr @44 '(' ll ' to ' ul ')'  ;  

		




Data a;set h;
/*
***********************Ingredients*******************************************
Here we obtain the weight W as the total sample size for the trial

The e1 and e2 are the observed failure probability multiplied by the weight

e1e2 is e1*e2

the RR estimate is be the ratio of the mean of e1 to the mean of e2

Note that e1 and e2 are the projected number of failures on the study if all 
subjects had received treatment 1 or all subjects received treatment 2 respectively.

The total of all the e1(e2) over all studies in the meta-analysis provides the number
of subjects projected to fail if all subjects received treatment 1(2) respectively
over all trials in the meta analysis.  The ratio of these two sums (or equivalently 
the ratio of the of the two means, is therefore a completely logical estimate of 
relative risk.

Random effects allows the presumption that these studies form a random sample of studies from 
a conceptual universe of studies.
 
Using the delta method, the square of the SE of log(e1/e2) is  V/M
with M=number of studies and 

V=[var(e1)/(mean(e1)**2)]+[var(e2)/{mean(e2**2)}]
-[2*cov(e1,e2)/[{mean(e1)*mean(e2)}]

SE{log(RR)}=sqrt(V/M)

****************************************************************************

*/
W=n1+n2;
e1=f1*w/n1;
e2=f2*w/n2;
e1e2=e1*e2;
proc means mean var n noprint;
var e1 e2 e1e2 ;
output out=new mean= e1 e2 e1e2
var= ve1 ve2 ve1e2 

n=m  me2 me1e2;
data new;set new;
data new;set new;
cove1e2=(e1e2-e1*e2)*(m/(m-1));

point_rr=e1/e2;
log_point_rr=log(e1/e2);
se2_logrr=(ve1/(e1*e1))+(ve2/(e2*e2))-2*(cove1e2/(e1*e2));
se_logrr=sqrt((se2_logrr)/m);
ll_rr=exp(tinv(.025,m-2)*se_logrr)*point_rr;
ul_rr=exp(tinv(.975,m-2)*se_logrr)*point_rr;
z_rr=log_point_rr/se_logrr;
p_rr=2*probt(-abs(z_rr),m-2);

data new;set new;

df=m-2;
file print;
put @10 'Study Size Weighted Low-Event Rate Binomial Meta Analysis';
put // @10 'Relative Risk Computation';
put // @10 'Values above 1.00 mean higher Failure Rate for Treatment Group';
put // @10 'Point Estimate=' point_rr best6.  @35 '95% CI= (' ll_rr best6.
'  to ' ul_rr best6. ' )';
if p_rr>0.0001 then
put // @10 'Two-sided P-Value=' p_rr z6.4 ' based on T with df=#studies-2=' 
df;
if p_rr<0.0001 then
put // @10 'Two-sided P-Value<0.0001 based on T with df=#studies-2=' df;


run;


;
*/
run;
%mend;

RUN;

/*
If you are creating an automatic .rtf file you need this statement to close  
it out.


*/
run;

data tab3;
length study $20.;
input study $ n1 mi1 cd1 n2 mi2 cd2;
cards;
49653/011 357 2 1 176 0 0 
49653/020 391 2 0 207 1 0 
49653/024 774 1 0 185 1 0
49653/093 213 0 0 109 1 0 
49653/094 232 1 1 116 0 0
100684 43 0 0 47 1 0 
49653/143 121 1 0 124 0 0 
49653/211 110 5 3 114 2 2 
49653/284 382 1 0 384 0 0 
712753/008 284 1 0 135 0 0 
AVM100264 294 0 2 302 1 1 
BRL49653C/185 563 2 0 142 0 0 
BRL49653/334 278 2 0 279 1 1 
BRL49653/347 418 2 0 212 0 0 
49653/015 395 2 2 198 1 0 
49653/079 203 1 1 106 1 1 
49653/080 104 1 0 99 2 0
49653/082 212 2 1 107 0 0 
49653/085 138 3 1 139 1 0
49653/095 196 0 1 96 0 0
49653/097 122 0 0 120 1 0
49653/125 175 0 0 173 1 0
49653/127 56 1 0 58 0 0
49653/128 39 1 0 38 0 0
49653/134 561 0 1 276 2 0
49653/135 116 2 2 111 3 1
49653/136 148 1 2 143 0 0
49653/145 231 1 1 242 0 0
49653/147 89 1 0 88 0 0
49653/162 168 1 1 172 0 0
49653/234 116 0 0 61 0 0
49653/330 1172 1 1 377 0 0
49653/331 706 0 1 325 0 0
49653/137 204 1 0 185 2 1
SB-712753/002 288 1 1 280 0 0 
SB-712753/003 254 1 0 272 0 0
SB-712753/007 314 1 0 154 0 0 
SB-712753/009 162 0 0 160 0 0
49653/132 442 1 1 112 0 0
AVA 394 1 1 124 0 0 
DREAM 2635 15 12 2634 9 10 
ADOPT 1456 27 2 2895 41 5
49653/044 101 0 0 51 0 0
49653/096 232 0 0 115 0 0
49653/109 52 0 0 25 0 0
49653/325 196 0 0 195 0 0
49653/282 70 0 0 75 0 0
49653/351 28 0 0 29 0 0
						

proc print;title 'Printout of 48 Studies';run;
;;;;

title 'Reanalysis of Nissen_Wolski 2007 rosiglitazone data: Myocardial Infarction';
%sgs (tab3,mi1,n1,mi2,n2);



run;
/*
This last piece of code will close out the rtf file you created
*/
run;
ods rtf close;
run;