FOCUS AREAS

Bayesian Hierarchical Modeling for Meta-Analysis


Contents | SAS Program | PDF

data multistudies;
   input studyid ctrl ctrlN trt trtN;
   datalines;
   1 3 39 3 38
   2 14 116 7 114
   3 11 93 5 69
   4 127 1520 102 1533
   5 27 365 28 355
   6 6 52 4 59
   7 152 939 98 945
   8 48 471 60 632
   9 37 282 25 278
   10 188 1921 138 1916
   11 52 583 64 873
   12 47 266 45 263
   13 16 293 9 291
   14 45 883 57 858
   15 31 147 25 154
   16 38 213 33 207
   17 12 122 28 251
   18 6 154 8 151
   19 3 134 6 174
   20 40 218 32 209
   21 43 364 27 391
   22 39 674 22 680
   ;
data multistudies; 
   set multistudies;
   logy=log(trt/(trtN-trt)/(ctrl/(ctrlN-ctrl)));
   sigma2=1/trt+1/(trtN-trt)+1/ctrl+1/(ctrlN-ctrl);
run;
proc mcmc data=multistudies outpost=nlout seed=276 nmc=50000 thin=5 
   monitor=(OR Pooled);
   array OR[22];
   parms  mu tau2;
   prior mu ~ normal(0, sd=3);
   prior tau2~ igamma(0.01,s=0.01);
   random theta ~n(mu, var=tau2) subject=studyid;
   OR[studyid]=exp(theta);
   Pooled=exp(mu);
   model logy ~ n(theta, var=sigma2);     
run;
%CATER(data=nlout, var=OR: Pooled);
proc mcmc data=multistudies outpost=blout seed=126 nmc=50000 thin=5 
   monitor=(OR Pooled);
   array OR[22]; 
   parms  mu_theta mu_phi s_theta s_phi;
   prior mu:~ normal(0,sd=3);
   prior s:~ igamma(0.01,s=0.01);
   random theta ~n(mu_theta, var=s_theta) subject=studyid;
   random phi ~n(mu_phi, var=s_phi) subject=studyid;
   p = logistic(theta + phi);
   q = logistic(phi);
   model trt ~ binomial(trtN, p);
   model ctrl ~ binomial(ctrlN, q);
   OR[studyid]=exp(theta);
   Pooled=exp(mu_theta);
run;

%CATER(data=blout, var=OR: Pooled);

data blout; 
   set blout;
   array OR {22} OR1-OR22;
   array y_{22} y_1-y_22;
   keep y_1-y_22;
   do i=1 to 22;
      y_{i}=OR{i};
   end;
run;
data all;
   set nlout blout;
   keep OR13 OR19  y_13 y_19;
run;
proc template;
   define statgraph OneByTwo;
   begingraph;
   layout lattice / rows=1 columns=2 rowgutter=10 columngutter=10;

   layout overlay;
      densityplot OR13 / kernel() lineattrs=GRAPHDATA2 name='n' legendlabel='normal';
      densityplot y_13 / kernel() lineattrs=GRAPHDATA1 name='b' legendlabel='binomial';
   endlayout;

   layout overlay;
      densityplot OR19 / kernel()lineattrs=GRAPHDATA2 name='n' legendlabel='normal';
      densityplot y_19 / kernel() lineattrs=GRAPHDATA1 name='b' legendlabel='binomial';
   endlayout;

   sidebar/ align=bottom;
      discretelegend "n" "b"/ border=1 ;  
   endsidebar;
      
   endlayout;
   endgraph;
   end;
run;
proc sgrender data=all template=OneByTwo;
run;