CODES 2000 User Forum -- Data Network Note #6
SAS Procedures for Meta-Analysis of Helmet Use Effects
Applies to: SAS 8.1 & 8.2.
Last updated: Thursday January 03, 2002.
SUMMARY
The effect of helmet use on inpatient charges is not identical in all states, nor does each state have a completely unique effect. Instead, we can model the effect of helmet use as normally distributed in each state with state means also normally distributed across the country. We can estimate the parameters of this hierarchical normal model from effects reported by the data network states. The SAS procedures described here implement the meta-analysis methodology presented by Gelman et al. (1995) in their Example 5.5 and are applied to their data scaled to plausible helmet use effects. These are not real CODES data.
PROCEDURE
- Simulate posterior distributions for model parameters and assess model fit.
The following SAS data steps and SAS sort, freq, print, gplot, and gchart procedures implement Example 5.5 from Gelman et al. (1995), Bayesian Data Analysis. We simulate posterior distributions for all model parameters and summarize 5,000 random draws as histograms and quantile tables. To assess model fit, we replicate the data 5,000 times and summarize the results as histograms.
*------------ Meta-Analysis of Helmet Use Effect ------------*
| Adapted from Gelman et al. (1995). |
| 5.5 Example: Combining information from educational |
| testing experiments at eight schools (from Rubin (1981)). |
*-------------------------------------------------------*;
data SATCoach Imputed;
input yBar1-yBar8 Sigma1-Sigma8 @@;
array yBar[8] yBar1-yBar8;
array Sigma[8] Sigma1-Sigma8;
array tauValue[100];
array tauProb[100];
c = 7.8348E-21;
grid0 = 0.1;
grid1 = 3000.1;
step = 50;
trials = 5000;
tauSum = 0;
tauCDF = 0;
index = 0;
do tau=grid0 to grid1 by step;
* Compute muHat and Vmu for posterior normal distribution of mu;
Sum1 = 0;
Sum2 = 0;
do i=1 to 8;
Sum1 = Sum1 + (yBar[i]/(Sigma[i]**2 + tau**2));
Sum2 = Sum2 + (1 /(Sigma[i]**2 + tau**2));
end;
muHat = Sum1 / Sum2;
Vmu = 1 / Sum2;
Product1 = 1;
do i=1 to 8;
Product1 = Product1 * (1 / Sqrt(Sigma[i]**2 + tau**2));
Product1 = Product1 * Exp(-(yBar[i] - muHat)**2 / (2 * (Sigma[i]**2 + tau**2)));
end;
tauSum = tauSum + step * (Sqrt(Vmu) * Product1);
tauPost = (1/c) * Sqrt(Vmu) * Product1;
tauCDF = tauCDF + (step * tauPost);
index=index+1;
tauValue[index] = tau;
tauProb[index] = tauCDF;
keep tau tauPost tauCDF tauSum;
output SATCoach;
end;
do draw=1 to trials;
tauMin = 0;
tauMax = 0;
ProbMin = 0;
ProbMax = 0;
tau = 0;
index = 1;
* Draw tau from its posterior distribution;
RanDraw = ranuni(1);
do while (RanDraw > tauProb[index] and index < 100);
ProbMin = tauProb[index];
tauMin = tauValue[index];
index = index + 1;
end;
if index < 100 then do;
ProbMax = tauProb[index];
tauMax = tauValue[index];
tau = tauMin + (tauMax - tauMin) * (RanDraw - ProbMin) / (ProbMax - ProbMin);
end;
* Draw mu from its posterior normal distribution given tau;
* Compute muHat and Vmu for posterior normal distribution of mu;Sum1 = 0;
Sum2 = 0;
do i=1 to 8;
Sum1 = Sum1 + (yBar[i]/(Sigma[i]**2 + tau**2));
Sum2 = Sum2 + (1 /(Sigma[i]**2 + tau**2));
end;
muHat = Sum1 / Sum2;
Vmu = 1 / Sum2;
sigmaMu = sqrt(Vmu);
mu = muHat + sigmaMu * rannor(1);
* Draw thetaJ from its posterior normal distribution given mu and tau;
* Calculate muJ and VJ for thetaJ given mu and tau; do J=1 to 8;VJ = 1 / ((1 / Sigma[J]**2) + (1 / tau**2));
sigmaJ = sqrt(VJ);
muJ = ((yBar[J] / Sigma[J]**2) + (mu / tau**2)) * VJ;
thetaJ = muJ + sigmaJ * rannor(1);
yJ = thetaJ + Sigma[J] * rannor(1);
keep draw tau mu J thetaJ yJ;
output Imputed;
end;
end;
datalines;
-2839 -794 275 -682 064 -063 -1801 -1216
1490 1020 1630 1100 940 1140 1040 1760
;
proc print data=SATCoach (firstobs=61);
title1 'Helmet Use Effect';
title2 'Integral of tau CDF';
var tau tauPost tauCDF tauSum;
run;
proc gplot data=SATCoach;
title2 'Marginal posterior density p(tau | y)';
plot tauPost*tau;
run;
proc gplot data=SATCoach;
title2 'CDF of p(tau | y)';
plot tauCDF*tau;
run;
*proc print data=Imputed;
* var draw tau mu J sigmaJ muJ thetaJ; *run; data Impute1;Set Imputed;
if J=1;
run;
proc gchart data=Impute1;
title2 'Simulated Population Std Dev';
vbar tau;
run;
proc gchart data=Impute1;
title2 'Simulated Population Mean';
vbar mu;
run;
proc sort data=Impute1;
by mu;
run;
data Quantile;
set Impute1;
if (mod(_N_,1000) in (25,250,500,750,975)) then do;
quantile = mod(_N_,1000)/10;
keep quantile mu;
output;
end;
run;
proc print data=Quantile;
title2 'mu Quantiles';
var quantile mu;
run;
proc sort data=Impute1;
by tau;
run;
data Quantile;
set Impute1;
if (mod(_N_,1000) in (25,250,500,750,975)) then do;
quantile = mod(_N_,1000)/10;
keep quantile tau;
output;
end;
run;
proc print data=Quantile;
title2 'tau Quantiles';
var quantile tau;
run;
proc sort data=Imputed;
by J thetaJ;
run;
proc gchart data=Imputed;
by J;
title2 'Simulated Effect thetaJ';
vbar thetaJ;
run;
data Quantile;
set Imputed;
if (mod(_N_,1000) in (25,250,500,750,975)) then do;
quantile = mod(_N_,1000)/10;
keep J quantile thetaJ;
output;
end;
run;
proc print data=Quantile;
title2 'thetaJ Quantiles';
var J quantile thetaJ;
run;
proc sort data=Imputed;
by J yJ;
run;
proc gchart data=Imputed;
by J;
title2 'Simulated Future Data yJ';
vbar yJ;
run;
data Quantile;
set Imputed;
if (mod(_N_,1000) in (25,250,500,750,975)) then do;
quantile = mod(_N_,1000)/10;
keep J quantile yJ;
output;
end;
run;
proc print data=Quantile;
title2 'yJ Quantiles';
var J quantile yJ;
run;
proc sort data=Imputed;
by draw J;
run;
data compare;
set Imputed;
retain;
if J = 1 then
theta1 = thetaJ;
if J = 3 then do;
theta3 = thetaJ;
if theta1 > theta3 then
Larger = 1;
else Larger = 3;
keep Larger;
output;
end;
run;
proc freq data=compare;
title2 'Comparison of theta1 vs. theta3';
tables Larger;
run;
data predict;
set Imputed;
array y[8];
retain;
if J < 8 then
y[J] = yJ;
if J = 8 then do;
y1 = y[1];
y2 = y[2];
y3 = y[3];
y4 = y[4];
y5 = y[5];
y6 = y[6];
y7 = y[7];
y8 = yJ;
yMax = max(y1,y2,y3,y4,y5,y6,y7,y8);
yMin = min(y1,y2,y3,y4,y5,y6,y7,y8);
yMean = mean(y1,y2,y3,y4,y5,y6,y7,y8);
yStd = std(y1,y2,y3,y4,y5,y6,y7,y8);
keep yMax yMin yMean yStd;
output;
end;
run;
proc gchart data=predict;
title2 'Posterior Predictive Distribution';
vbar yMax;
vbar yMin;
vbar yMean;
vbar yStd;
run;
SAS is a registered trademark of SAS Institute, Inc.