I'm trying to come up with funnel plot from meta-analysis using macro recommended by Rendina-Gobioff et al., PUB_BIAS: A SAS® Macro for Detecting, Publication Bias in Meta-Analysis.
Below macro runs fine with no error in log. However, funnel plot is not produced. Any hints or suggestions appreciated.
%macro symsize;
data _null_;
set metadat;
retain sizeh1-sizeh%eval(&n_stud+2)
fontv1-fontv%eval(&n_stud+2);
length fontv1-fontv%eval(&n_stud+2) $ 20;
array sizes sizeh1-sizeh%eval(&n_stud+2);
array fvs $ fontv1-fontv%eval(&n_stud+2);
do i=1 to (&n_stud)+2;
if i=y then do;
sizes{i}=size*11/&maxsize;
fvs{i}='font=specialu v=K';
if i=&n_stud+2 then output;
end;
end;
%do i=1 %to (&n_stud)+2;
call symput("sh&i",trim(left(put(sizeh&i,6.2))));
call symput("fv&i",trim(left(put(fontv&i,20.))));
%global sh&i fv&i;
%end;
run;
%mend symsize;
%macro doall (study);
%do %while (&study le (&n_stud+2));
if y=&study then do;
xx&study=lower95;
yy&study=y+0.2; output;
yy&study=y-0.2; output;
yy&study=y;
output;
xx&study=upper95; output;
yy&study=y+0.2;output;
yy&study=y-0.2; output;
end;
%let study=%eval(&study+1);
%end;
%mend doall;
%macro syms ;
%do i=1 %to (&n_stud)+2;
symbol&i &&fv&i l=1 interpol=none h=&&sh&i color=black;
%end;
%mend syms;
%macro plotpts ;
%do i=1 %to (&n_stud + 2);
yy&i*xx&i=%eval(&n_stud+3)
%end;
%do i=1 %to (&n_stud + 2);
y&i*x&i=&i
%end;
%mend plotpts;
%Macro PubBias(N1=size1,N2=size2,di=effsize,studyID=study,ModelType=1,Dataset=MetaPub);
proc iml;
* +----------------------------------------------------------------------------------------+
Read data from regular SAS into IML
+----------------------------------------------------------------------------------------+;
use &Dataset;
read all var{&N1} into n1_vec;
read all var{&N2} into n2_vec;
read all var{&di} into di_vec;
read all var{&n1 &n2} into n_vec;
k= nrow(di_vec);
* +----------------------------------------------------------------------------------------+
Subroutine to calculate weighted mean effect size,
standard error, and confidence interval for mean.
Inputs to the subroutine are
di_vec - column vector of effect sizes (d)
var_di - column vector of estimation errors
tau2 - scalar estimate of RANDOM EFFECTS variance (set to 0 if Modeltype-fixed(0))
Outputs are
d_mean = weighted mean d value
resum_wt = scalar, sum of the weights
5
vi_star = column vector of total variance for each study
d_SE = standard error of d
upper95, lower95 = endpoints of 95% CI
+----------------------------------------------------------------------------------------+;
start mean_d(di_vec,var_di,tau2,d_mean,resum_wt,vi_star,d_SE,upper95,lower95);
k = nrow(di_vec);
d_mean = 0;
resum_wt = 0;
vi_star = J(k,1,0);
do i = 1 to k;
d_mean = d_mean + di_vec[i,1]/(var_di[i,1]+tau2);
resum_wt = resum_wt + (var_di[i,1]+tau2)##-1;
vi_star[i,1] = var_di[i,1]+tau2;
end;
d_mean = d_mean/resum_wt;
d_SE = SQRT(resum_wt##-1);
upper95 = d_mean + 1.96#d_SE;
lower95 = d_mean - 1.96#d_SE;
finish;
* +----------------------------------------------------------------------------------------+
Subroutine to calculate the Q test
of homogeneity.
Inputs to the subroutine are
di_vec - column vector of effect sizes (d)
n_vec - matrix (k X 2) of sample sizes
corresponding to each effect size
Outputs are
QQ = the obtained value of Q
prob_qq1 = chi-square probability associated with QQ
var_di = column vector of variances of effect sizes
+----------------------------------------------------------------------------------------+;
start calcq(di_vec,n_vec,qq,prob_qq1,var_di);
k = nrow(di_vec);
var_di=J(k,1,0);
do i = 1 to k;
var_di[i,1] = ((n_vec[i,1]+n_vec[i,2])/(n_vec[i,1]#n_vec[i,2])) +
((di_vec[i,1]##2)/(2#(n_vec[i,1]+n_vec[i,2])));
end;
d_plus = 0;
fesum_wt = 0;
do i = 1 to k;
d_plus = d_plus + di_vec[i,1]/var_di[i,1];
fesum_wt = fesum_wt + var_di[i,1]##-1;
end;
d_plus = d_plus/fesum_wt;
QQ = 0;
do i = 1 to k;
QQ = QQ + ((di_vec[i,1] - d_plus)##2/var_di[i,1]);
end;
prob_qq1 = 1 - PROBCHI(QQ,k-1);
finish;
*+----------------------------------------------------------------------------------------+
Subroutine to calculate OLS and WLS tests of models.
6
Inputs to the subroutine are
di_vec - column vector of effect sizes (d)
n_vec - matrix (k X 2) of sample sizes
corresponding to each effect size
X_Matrix - Matrix of potential moderator variables
For tests of publication bias, vi are predictors
vi - reciprocals of variances
Outputs are
B_wls - regression weights for WLS
SE_B - Standard errors of the WLS weights
B_ols - regression weights for OLS
SE_B_ols - Standard errors of the OLS weights
+----------------------------------------------------------------------------------------+;
start calcreg(di_vec,n_vec,X_Matrix,vi,B_wls,SE_B,B_ols,SE_B_ols);
k = nrow(di_vec);
X = J(k,1,1)||X_Matrix;
B_wls = INV(X`*DIAG(vi)*X)*X`*DIAG(vi)*di_vec;
cov_b = INV(X`*DIAG(vi)*X);
SE_B = SQRT(vecdiag(cov_b));
B_ols =INV(X`*X)*X`*di_vec;
cov_b = INV(X`*X);
SE_B_ols = SQRT(vecdiag(cov_b));
finish;
* +----------------------------------------------------------------------------------------+
Subroutine Kendall
Computes the Kendall Tau for Norman Cliff ordinal level analyses.
Arguments to the subroutines are:
A B = vectors of observed data for the two variables
(A will be the observed tx effects
B will be the variance of the tx effects)
N = sample size
Returned are:
T_AB = Kendall Tau Coefficient Y and X1
UNTIE_A = proportion of scores that are not tied on A
UNTIE_B = proportion of scores that are not tied on B
VART_AB = VARIANCE of Y AND X1
Z_TEST = obtained value of Z for test of Kendall Tau Coefficient
+----------------------------------------------------------------------------------------+;
START KENDALL(A,B,N,T_AB,untie_A,untie_B,VART_AB,Z_TEST);
DOM_MTXA = J(N,N,0);
ties_A = 0;
counts_A = 0;
do i = 1 to N;
do j = 1 to N;
if A[i,1] > A[j,1] then do;
DOM_MTXA[i,j] = 1;
end;
if A[i,1] < A[j,1] then do;
DOM_MTXA[i,j] = -1;
end;
if A[i,1] = A[j,1] then do;
ties_A = ties_A + 1;
end;
counts_A = counts_A + 1;
end;
end;
untie_A = 1 - (ties_A - N)/(counts_A - N);
DOM_MTXB = J(N,N,0);
ties_B = 0;
counts_B = 0;
do i = 1 to N;
do j = 1 to N;
if B[i,1] > B[j,1] then do;
DOM_MTXB[i,j] = 1;
end;
if B[i,1] < B[j,1] then do;
DOM_MTXB[i,j] = -1;
end;
if B[i,1] = B[j,1] then do;
ties_B = ties_B + 1;
end;
counts_B = counts_B + 1;
end;
end;
untie_B = 1 - (ties_B - N)/(counts_B - N);
DOM_MTXD = J(N,N,0);
do i = 1 to N;
do j = 1 to N;
DOM_MTXD[i,j] = DOM_MTXA[i,j]#DOM_MTXB[i,j];
end;
end;
MTXD_sum = DOM_MTXD[,+];
T_AB = MTXD_sum[+,] #(1/(n#(n-1)));
MTX_F = J(N,1,0);
do i = 1 to N;
MTX_F[i,1] = (MTXD_sum [i,1]/(n-1));
end;
MTX_G = J(N,1,0);
do i = 1 to N;
MTX_G[i,1] = (MTX_F[i,1] - T_AB)##2;
end;
MTXG_sum = 1/(n-1)#(MTX_G[+,]);
MTX_H = J(N,N,0);
do i = 1 to N;
do j = 1 to N;
MTX_H[i,j] = DOM_MTXD[i,j]#DOM_MTXD[i,j];
end;
end;
MTXH_sum = MTX_H[+,+];
NUMER_AB = MTXH_sum - ((n) # (n-1) # (T_AB#T_AB));
DENOM_AB = n#(n-1) -1;
VART_AB = NUMER_AB/DENOM_AB;
vart_ab = ((4#(n-2)#(MTXG_sum))+(2#VART_AB)) / (n#(n-1));
IF vart_ab > 0 THEN DO ;
Z_TEST = (T_AB /SQRT(vart_ab));
end;
if vart_ab =0 then do;
Z_TEST = 5.00;
END;
FINISH;
* +----------------------------------------------------------------------------------------+
Main program
+----------------------------------------------------------------------------------------+;
* +----------------------------------------------------------------------------------------+
Calculate Q
+----------------------------------------------------------------------------------------+;
run calcq(di_vec,n_vec,qq,prob_qq1,var_di);
CC = (J(1,K,1)*var_di##-1) - ((J(1,K,1)*var_di##-2) /
(J(1,K,1)*var_di##-1));
Tau2 = (QQ - (K - 1)) / cc;
if tau2 < 0 then tau2 = 0;
if &modeltype=0 then tau2=0;
* +----------------------------------------------------------------------------------------+
Calculate weighted mean effect size and confidence interval
+----------------------------------------------------------------------------------------+;
run mean_d(di_vec,var_di,tau2,d_mean,resum_wt,vi_star,d_SE,upper95,lower95);
vi = J(k,1,0);
vi_inv = j(k,1,0);
mean = 0;
weight = 0;
t_star = J(k,1,0);
Vi_stand = J(k,1,0);
dev_di = J(k,1,0);
Root_Vi=J(k,1,0);
Egger_z = J(k,1,0);
REJ_TRIM = J(3,1,0);
do i = 1 to k;
vi[i,1] = vi_star[i,1];
vi_inv[i,1] = vi[i,1]##-1;
mean = d_mean;
weight = resum_wt;
dev_di[i,1] = ABS(di_vec[i,1] - mean);
Vi_stand[i,1] = vi[i,1] - (weight##-1);
t_star[i,1] = (di_vec[i,1] - mean)/SQRT(Vi_stand[i,1]);
Root_Vi[i,1] = vi[i,1]##-0.5;
Egger_z[i,1] = di_vec[i,1]#Root_Vi[i,1];
end;
* +----------------------------------------------------------------------------------------+
Calculate Egger OLS Regression (Precision as predictor)
+----------------------------------------------------------------------------------------+;
run calcreg(Egger_z,n_vec,Root_Vi,vi_inv,B_wls,SE_B,B_ols,SE_B_ols);
eggt_ols = B_ols[1,1]/SE_B_ols[1,1];
eggPR_tOLS =2#(1-probt(abs(eggt_ols),k-2));
* +----------------------------------------------------------------------------------------+
Calculate Begg Rank Correlation(Variance as predictor)
+----------------------------------------------------------------------------------------+;
run KENDALL(t_star,vi,k,T_X1Y,UT_A,UT_B,VART_X1Y,Z_TEST);
BeggV_z = Z_TEST;
BeggVpr_z =2#(1-probnorm(abs(Z_test)));
* +----------------------------------------------------------------------------------------+
Calculate Begg Rank Correlation (Sample size as predictor)
+----------------------------------------------------------------------------------------+;
total_n = n_vec * J(2,1,1);
run KENDALL(t_star,total_n,k,T_X1Y,UT_A,UT_B,VART_X1Y,Z_TEST);
BeggN_z = Z_TEST;
BeggNpr_z =2#(1-probnorm(abs(Z_test)));
* +----------------------------------------------------------------------------------------+
Calculate Funnel Plot WLS Regression (Sample size as predictor)
+----------------------------------------------------------------------------------------+;
run calcreg(di_vec,n_vec,total_n,vi_inv,B_wls,SE_B,B_ols,SE_B_ols);
Funt_WLS = B_wls[2,1]/SE_B[2,1];
FunPR_tWLS =2#(1-probt(abs(Funt_WLS),k-2));
* +----------------------------------------------------------------------------------------+
Calculate Trim and Fill
+----------------------------------------------------------------------------------------+;
dev_rank = rank(dev_di);
do i = 1 to k;
if di_vec[i,1] < mean then dev_rank[i,1] = -1#dev_rank[i,1];
end;
r=-1#MIN(dev_rank);
gamma=k-r;
ro=gamma-1;
r2=MAX(dev_rank);
gamma2=k-r2;
ro2=gamma2-1;
if ro > 3 then REJ_TRIM[1,1] = REJ_TRIM[1,1] +1;
if ro2 > 3 then REJ_TRIM[2,1] = REJ_TRIM[2,1] +1;
* +------------------------------------------------------------------+
Either tail (note: alpha is .10 for this
+------------------------------------------------------------------+;
if (ro > 3 | ro2 > 3) then REJ_TRIM[3,1] = REJ_TRIM[3,1] +1;
Right=REJ_TRIM[1,1]; If Right=1 then Right='Yes'; Else Right='No';
Left=REJ_TRIM[2,1];If Left=1 then Left='Yes'; Else Left='No';
Both=REJ_TRIM[3,1];If Both=1 then Both='Yes'; Else Both='No';
* +----------------------------------------------------------------------------------------+
Print output results
+----------------------------------------------------------------------------------------+;
file print;
put
@1'-------------------------------------------------------------------------------------'//
@1 'Meta-Analysis: Descriptive Information'/
@1 '--------------------------------------'//
@5'Number of Studies'@60 k//
@5'Mean Effect Size' @60 d_mean 10.4/
@7'Confidence Band'/
@10 '95% Lower Limit' @60lower95 10.4/
@10'95% Upper Limit'@60 upper95 10.4//
@5 'Test for Homogeneity'@60'Chi Square'@75'Probability'/
@5 '--------------------' @60'----------'@75'-----------'/
10
@5'Q test'@62 QQ 8.4 @75 prob_qq1 10.4//
@5'Random Effects Variance Component'@62 Tau2 8.4///
@1'--------------------------------------------------------------------------------------'//
@1 'Meta-Analysis: Tests of Publication Bias'/
@1 '----------------------------------------'//
@5 'Egger Regression' @60't value'@75'Probability'/
@5 '----------------' @60'-------'@75'-----------'/
@5 'Egger OLS Regression' @60 eggt_OLS 7.4 @75 eggPR_tOLS 10.4///
@5 'Begg Rank Correlation' @60'z value'@75'Probability'/
@5 '---------------------' @60'-------'@75'-----------'/
@5 'Begg Rank Correlation (Predictor=Variance)' @60 beggV_z 7.4 @75 beggVpr_z 10.4/
@5 'Begg Rank Correlation (Predictor=Sample Size)' @60 beggN_z 7.4 @75 beggNpr_z 10.4///
@5 'Funnel Plot Regression' @60't value'@75'Probability'/
@5 '----------------------' @60'-------'@75'-----------'/
@5 'Funnel Plot WLS Regression' @60 Funt_WLS 7.4 @75 FunPR_tWLS 10.4///
@5 'Trim and Fill'@60'Publication Bias Present'/
@5 '-------------'@60'------------------------'/
@7 'Right Tail' @70 Right/
@7 'Left Tail' @70 Left/
@7 'Both Tails' @70 Both/
@1'--------------------------------------------------------------------------------------'//;
* +----------------------------------------------------------------------------------------+
Send summary information to regular SAS for plots
+----------------------------------------------------------------------------------------+;
sum_n = n_vec*J(2,1,1);
total_n = sum(sum_n);
outvector = d_mean||lower95||upper95||total_n;
create j1 from outvector;
append from outvector;
quit;
data meta2;
set &Dataset;
SE = SQRT((&n1 + &n2) / (&n1 * &n2) + &di**2 / (2*(&n1 + &n2)));
lower95 = &di - 1.96*SE;
upper95 = &di + 1.96*SE;
size = &n1 + &n2;
y = _n_;
data total1;
set j1;
rename col1 = &di col2 = lower95 col3 = upper95 col4 = size;
length &StudyID $ 12;
y = 99;
&StudyID = 'Total';
data total2;
y = 98;
data total;
set total1 total2;
axis1 label=(height=2.9
font='Zapf' 'Effect Size')
minor=none
value=(height=2.5 font='Zapf');
axis2 label=(height=2.9 angle = 90
font='Zapf' 'Sample Size')
minor=none
value=(height=2.5 font='Zapf');
title1 font = 'Zapf' height = 2.9 'Funnel Plot' ;
symbol1 interpol=none
value=circle
height=3
cv=black
ci=black
co=black
width=2;
proc gplot data=meta2;
plot size*&di /haxis=axis1 vaxis=axis2 frame ;
run;
proc means data=meta2 noprint;
var size;
output out=meanout sum=sum max=max;
data _null_;
set meanout end=eof;
if eof then do;
call symput("n_stud",trim(left(put(_freq_,8.))));
call symput("n_subs",trim(left(put(sum,8.))));
call symput("maxsize",trim(left(put(max,8.))));
end;
proc sort data=meta2; by &di upper95;
run;
data metadat;
set total(in=a) meta2(in=b);
length metastr $ 200;
retain metastr studcnt;
metanum=y;
y=_n_;
if _n_=1 then do;
metastr="0=' ' %eval(&n_stud+3)=' '";
studcnt=&n_stud+2;
end;
else studcnt=studcnt-1;
do i=1 to %eval(&n_stud+2);
if i=studcnt then metastr=" " || trim(left(metastr)) || ' ' || trim(left(y))
||"='" || trim(left(&StudyID)) || "' ";
end;
if &StudyID=' ' then y=.;
yo=lower95; xo=upper95;
output metadat;
if _n_=%eval(&n_stud+2) then do;
call symput("metastr",trim(metastr));
end;
run;
proc format;
value stfmt &metastr;
run;
data alldata;
set metadat;
%doall (1)
%symsize
%syms
data final(drop=i);
set alldata;
array xarray{*} x1-x%eval(&n_stud+2);
array yarray{*} y1-y%eval(&n_stud+2);
do i=1 to (&n_stud)+2;
if i=y then do;
xarray{i}=&di;
yarray{i}=y;
end;
end;
format yy1-yy%eval(&n_stud+2) stfmt.;
run;
axis1 label=(height=2.9
font='Zapf' 'Effect Size')
minor=none
value=(height=2.5 font='Zapf');
axis2 label=(height=2.9 angle = 90
font='Zapf' 'Study')
minor=none
order=( 0 to %eval(&n_stud+3) by 1)
value=(height=2.5 font='Zapf');
title1 font = 'Zapf' height = 2.9 'Forest Plot' ;
proc gplot data=final;
plot
%plotpts / overlay haxis=axis1 vaxis=axis2 frame href=0;
symbol%eval(&n_stud+3) f=marker v=none l=1 w=1 i=join;
symbol1 font=marker v=P l=1 h=2.7 interpol=none;
run;
%Mend PubBias;
quit;
Data MetaPub;
Input size1 1-3 size2 5-7 effsize 9-12 study $ 14-25;
Datalines;
25 30 0.75 Able 1986
100 125 0.45 Bonk 1994
250 250 0.25 Carson 1990
90 150 0.35 Diddle 1993
50 60 0.70 Efron 1991
180 130 0.39 Flipper 1989
68 82 0.50 Goober 1975
170 200 0.30 Halcyon 2002
110 90 0.42 Illy 2004
45 45 0.65 Jersey 2002
;
%PubBias (N1=size1,N2=size2,di=effsize,studyID=study,ModelType=1,Dataset=MetaPub);
run;
Click image to register for webinar
