%meta_rr
%meta_or
%meta_rd
%my_meta_norm
%ForestMacro
%PubBias
%macro pseudo_rsq(log_lh=, log_lh_i=, k=, obs=);
proc iml;
value = 1-((&log_lh.)/&log_lh_i.);
/* cox_snell_rsq = 1 - (&log_lh_i./&log_lh.)##(2/&obs.);*/
title "McFadden's R-Squared";
print value;
title;
quit;
%mend;
%macro rsq(data = , p=);
%let open_d = %sysfunc (open(&data.));
%let obs = %sysfunc(attrn(&open_d., nobs));
%let close_d = %sysfunc(close(&open_d.));
proc iml;
use &data.;
read all var _all_;
close &data.;
SSres = sum((rri-pred)##2);
/* SSt = sum(rri##2) - &obs.#(mean(rri)##2);*/
sst = sum((rri-mean(rri))##2);
ADJ_R_SQ = 1 - ((ssres/(&obs.-&p.)) / (sst/(&obs.-1)));
print ADJ_R_SQ;
quit;
%mend;
%macro meta_rr(data = , ai = , bi = , ci = , di = ,references=);
data new_data dorp_data;
set &data.;
if &ai. = 0 or &bi. = 0 or &ci. = 0 or &di. = 0 then do;
output dorp_data;
delete;
end;
if ai = bi and ci = di then do;
output dorp_data;
delete;
end;
output new_data;
run;
proc iml;
use new_data;
read all var _all_;
close new_data;
%let open_d = %sysfunc (open(new_data));
%let obs_ = %sysfunc(attrn(&open_d., nobs));
%let close_d = %sysfunc(close(&open_d.));
%let obs_1 = %sysevalf(&obs_ - 1);
n1 = &ai.+&ci.;
n2 = &bi.+&di.;
m1 = &ai.+&bi.;
m2 = &ci.+&di.;
N_total = &ai. +&bi.+&ci.+&di.;
/*99% 2.576*/
eer = &ai. / n1; /*98% 2.326*/
cer = &bi. / n2;
nnt = 1/abs(eer-cer);
nnt_ci_l = 1/((eer-cer)+1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
nnt_ci_u = 1/((eer-cer)-1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
nnh = 1/abs(eer-cer);
nnh_ci_l = 1/((eer-cer)+1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
nnh_ci_u = 1/((eer-cer)-1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
/*RR*/
pti = &ai./(&ai.+&ci.);
pci = &bi./(&bi.+&di.);
rri = pti/pci;
ln_rri = log(rri);
var_ln_rri = (1-pti)/(n1#pti) + (1-pci)/(n2#pci);
wi = 1/var_ln_rri;
rr_pooled = exp(sum(wi#log(rri))/sum(wi));
se_ln_rr_pooled = sqrt(1/(sum(wi)));
se_ln_rr=((1-pti)/(&ai.+&ci.)/pti+(1-pci)/(&bi.+&di.)/pci)##0.5;
rri_ci_l = exp(log(rri)-(1.96#se_ln_rr));
rri_ci_u = exp(log(rri)+(1.96#se_ln_rr));
rr_pooled_ci_l = exp(log(rr_pooled)-(1.96#se_ln_rr_pooled));
rr_pooled_ci_u = exp(log(rr_pooled)+(1.96#se_ln_rr_pooled));
print rri rri_ci_l rri_ci_u rr_pooled rr_pooled_ci_l rr_pooled_ci_u se_ln_rr_pooled ;
Zrr = log(rr_pooled)/se_ln_rr_pooled; /*ho: rr=1 reject h0 if |zrr|>=1.96*/
/* zrr_pv2 = 1 - cdf("Normal", Zrr); */
zrr_pv = (1-probnorm(abs(Zrr)))#2;
print zrr zrr_pv;
Qrr = sum(wi#((log(rri)-log(rr_pooled))##2));/*ho: rr1=rr2....=rr reject h0 qrr_pv>qrr_pv_005*/
/* qrr_pv_005 = cinv(0.95, 8);*/
qrr_pv = (1-probchi(Qrr, &obs_1.));
print Qrr qrr_pv;
i2 = (Qrr - &obs_1.)/ qrr;
print i2;
/* Random Effect */
t_2 = max(0, (Qrr-(&obs_1.)) / (sum(wi)- (sum(wi##2)) /sum(wi))) ;
wi_ = 1/(var_ln_rri+t_2);
theta_ = sum(wi_#rri) / sum(wi_);
se_theta_ = 1 / sqrt(sum(wi_));
theta__ci_l = theta_ - 1.96 # se_theta_;
theta__ci_u = theta_ + 1.96 # se_theta_;
r_pooled_rr = exp(sum(wi_#log(rri))/sum(wi_));
se_r_pooled_rr = 1 / sqrt(sum(wi_));
r_pooled_rr_ci_l = r_pooled_rr - 1.96#se_r_pooled_rr;
r_pooled_rr_ci_u = r_pooled_rr + 1.96#se_r_pooled_rr;
print t_2 wi_ theta_ se_theta_ theta__ci_l theta__ci_u ;
print r_pooled_rr se_r_pooled_rr r_pooled_rr_ci_l r_pooled_rr_ci_u;
/*theta_的有效性檢驗， h0 : theta_ =0 */
utheta_ = ((sum(wi_ # rri))##2 ) / sum(wi_);
utheta__pv = (1-probchi(utheta_, 1));
zutheta_ = theta_ / se_theta_;
zutheta__pv = (1-probnorm(abs(zutheta_)))#2;
print utheta_ utheta__pv zutheta_ zutheta__pv;
median_ln_rri = median(ln_rri);
wi_pct = wi/sum(wi);
sum_w = sum(wi_pct);
create output var {rri rri_ci_l rri_ci_u rr_pooled rr_pooled_ci_l rr_pooled_ci_u se_ln_rr zrr Qrr
qrr_pv ln_rri wi n1 n2 median_ln_rri i2 wi_pct wi_};
append;
close output;
GroupId =3;
create pooled var {rr_pooled rr_pooled_ci_l rr_pooled_ci_u GroupId};
append;
close output;
create ran_pooled var {r_pooled_rr r_pooled_rr_ci_l r_pooled_rr_ci_u GroupId};
append;
close output;
quit;title;
data output2;
merge output new_data;
obs = _n_;
GroupId = 1;
run;
data output3;
set output2
pooled (rename=(RR_POOLED=rri rr_pooled_ci_l=rri_ci_l rr_pooled_ci_u=rri_ci_u))
ran_pooled (rename=(r_pooled_rr=rri r_pooled_rr_ci_l=rri_ci_l r_pooled_rr_ci_u=rri_ci_u));
if _n_ = &obs_.+1 then &references. = "Pooled_F";
if _n_ = &obs_.+2 then &references. = "Pooled_R";
run;
proc sgplot data = output2;
HIGHLOW y=obs HIGH=rri_ci_u LOW=rri_ci_l / close = rri ;
REFLINE 1 / label = "1" AXIS= X;
yaxis REVERSE;
quit;
proc print data = output2(keep = &references. rri rri_ci_u rri_ci_l);
run;
/*
proc sgplot data = output2;
STEP x = obs y= wi / DATALABEL = wi;
quit;
*/
proc sgplot data = output2;
scatter y = se_ln_rr x = ln_rri / datalabel = obs;
yaxis REVERSE;
refline median_ln_rri / axis = x;
quit;
%mend;
%macro meta_or(data = , ai = , bi = , ci = , di = ,references=);
data new_data dorp_data;
set &data.;
if &ai. = 0 or &bi. = 0 or &ci. = 0 or &di. = 0 then do;
output dorp_data;
delete;
end;
if ai = bi and ci = di then do;
output dorp_data;
delete;
end;
output new_data;
run;
proc iml;
use new_data;
read all var _all_;
close new_data;
%let open_d = %sysfunc (open(new_data));
%let obs_ = %sysfunc(attrn(&open_d., nobs));
%let close_d = %sysfunc(close(&open_d.));
%let obs_1 = %sysevalf(&obs_ - 1);
n1 = &ai.+&ci.;
n2 = &bi.+&di.;
m1 = &ai.+&bi.;
m2 = &ci.+&di.;
N_total = &ai. +&bi.+&ci.+&di.;
/*99% 2.576*/
eer = &ai. / n1; /*98% 2.326*/
cer = &bi. / n2;
nnt = 1/abs(eer-cer);
nnt_ci_l = 1/((eer-cer)+1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
nnt_ci_u = 1/((eer-cer)-1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
nnh = 1/abs(eer-cer);
nnh_ci_l = 1/((eer-cer)+1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
nnh_ci_u = 1/((eer-cer)-1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
/*OR*/
oddsti = &ai./&ci.;
oddsci = &bi./&di.;
ori = oddsti/oddsci;
/*Woolf's Logit estimate of OR*/
var_ln_ori = 1/&ai. + 1/&bi. + 1/&ci. + 1/&di.;
wi = 1 / var_ln_ori;
ln_orw = sum(wi#log(ori))/sum(wi);
orw = exp(ln_orw);
ln_ori = log(ori);
se_ln_or = sqrt(var_ln_ori);
ln_ori_ci_l = ln_ori - 1.96#sqrt(var_ln_ori);
ln_ori_ci_u = ln_ori + 1.96#sqrt(var_ln_ori);
ori_ci_l = exp(ln_ori_ci_l);
ori_ci_u = exp(ln_ori_ci_u);
orw_ci_l = exp(ln_orw-1.96#sqrt(1/sum(wi)));
orw_ci_u = exp(ln_orw+1.96#sqrt(1/sum(wi)));
print ori ori_ci_l ori_ci_u ln_ori ln_ori_ci_l ln_ori_ci_u orw orw_ci_l orw_ci_u;
/* orw顯著性 h0:ORw=OR0 (e.g. or0=1)*/
or0=1;
Uw = sum(wi)#((ln_orw-log(or0))##2);
Uw_pv = (1-probchi(Uw, 1));
print Uw Uw_pv;
/*ori同質性檢定 H0: or1=or2=...=or*/
Qw = sum(wi#(log(ori)##2)) - sum(wi)#(ln_orw##2);
Qw_pv = (1-probchi(Qw, &obs_1.));
print Qw Qw_pv;
i2 = (qw - &obs_1. )/ qw;
print i2;
/* Random Effect */
t_2 = max(0, (Qw-(&obs_1.)) / (sum(wi)- (sum(wi##2)) /sum(wi))) ;
wi_ = 1/(var_ln_ori+t_2);
theta_ = sum(wi_#ori) / sum(wi_);
se_theta_ = 1 / sqrt(sum(wi_));
theta__ci_l = theta_ - 1.96 # se_theta_;
theta__ci_u = theta_ + 1.96 # se_theta_;
r_pooled_or = exp(sum(wi_#log(ori))/sum(wi_));
se_r_pooled_or = 1 / sqrt(sum(wi_));
r_pooled_or_ci_l = r_pooled_or - 1.96#se_r_pooled_or;
r_pooled_or_ci_u = r_pooled_or + 1.96#se_r_pooled_or;
print t_2 wi_ theta_ se_theta_ theta__ci_l theta__ci_u ;
print r_pooled_or se_r_pooled_or r_pooled_or_ci_l r_pooled_or_ci_u;
/*theta_的有效性檢驗， h0 : theta_ =0 */
utheta_ = ((sum(wi_ # ori))##2 ) / sum(wi_);
utheta__pv = (1-probchi(utheta_, 1));
zutheta_ = theta_ / se_theta_;
zutheta__pv = (1-probnorm(abs(zutheta_)))#2;
print utheta_ utheta__pv zutheta_ zutheta__pv;
median_ln_ori = median(ln_ori);
wi_pct = wi/sum(wi);
create output var { i2
ori ori_ci_l ori_ci_u ln_ori ln_ori_ci_l ln_ori_ci_u orw orw_ci_l orw_ci_u se_ln_or
Uw Uw_pv Qw
Qw_pv wi median_ln_ori n1 n2 wi_pct wi_};
append;
close output;
GroupId=3;
create pooled var {orw orw_ci_l orw_ci_u GroupId};
append;
close output;
create ran_pooled var {r_pooled_or r_pooled_or_ci_l r_pooled_or_ci_u GroupId};
append;
close output;
quit;title;
data output2;
merge output new_data;
obs = _n_;
GroupId = 1;
run;
data output3;
set output2
pooled (rename=(orw=ori orw_ci_l=ori_ci_l orw_ci_u=ori_ci_u))
ran_pooled (rename=(r_pooled_or=ori r_pooled_or_ci_l=ori_ci_l r_pooled_or_ci_u=ori_ci_u));
if _n_ = &obs_.+1 then &references. = "Pooled_F";
if _n_ = &obs_.+2 then &references. = "Pooled_R";
run;
proc sgplot data = output2;
HIGHLOW y=obs HIGH=ori_ci_u LOW=ori_ci_l / close = ori;
REFLINE 1 / label = "1" AXIS= X;
yaxis REVERSE;
quit;
proc print data = output2(keep = &references. ori ori_ci_u ori_ci_l);
run;
/*
proc sgplot data = output2;
STEP x = obs y= wi_or / DATALABEL = wi_or;
quit;
*/
proc sgplot data = output2;
scatter y = se_ln_or x = ln_ori / datalabel = obs;
yaxis REVERSE;
refline median_ln_ori / axis = x;
quit;
%mend;
%macro meta_rd(data = , ai = , bi = , ci = , di = ,references=);
data new_data dorp_data;
set &data.;
if ai = 0 and bi = 0 then do;
output dorp_data;
delete;
end;
if ai = bi and ci = di then do;
output dorp_data;
delete;
end;
output new_data;
run;
proc iml;
use new_data;
read all var _all_;
close new_data;
%let open_d = %sysfunc (open(new_data));
%let obs_ = %sysfunc(attrn(&open_d., nobs));
%let close_d = %sysfunc(close(&open_d.));
%let obs_1 = %sysevalf(&obs_ - 1);
n1 = &ai.+&ci.;
n2 = &bi.+&di.;
m1 = &ai.+&bi.;
m2 = &ci.+&di.;
N_total = &ai. +&bi.+&ci.+&di.;
/*99% 2.576*/
eer = &ai. / n1; /*98% 2.326*/
cer = &bi. / n2;
nnt = 1/abs(eer-cer);
nnt_ci_l = 1/((eer-cer)+1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
nnt_ci_u = 1/((eer-cer)-1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
nnh = 1/abs(eer-cer);
nnh_ci_l = 1/((eer-cer)+1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
nnh_ci_u = 1/((eer-cer)-1.96*sqrt((eer#(1-eer)/n1)+(cer#(1-cer)/n2)));
/*RD*/
pti = &ai./(&ai.+&ci.);
pci = &bi./(&bi.+&di.);
RDi = pti - pci;
var_rdi = (pti#(1-pti))/n1 + (pci#(1-pci))/n2;
se_rd = sqrt(var_rdi);
wi = 1/var_rdi;
RD_pooled = sum(wi#rdi) / sum(wi);
var_rd_p = 1/sum(wi);
se_rd_p = 1/sqrt(sum(wi));
RDi_ci_l = RDi - 1.96#sqrt(var_rdi);
RDi_ci_u = RDi + 1.96#sqrt(var_rdi);
RDp_ci_l = RD_pooled - 1.96#se_rd_p;
RDp_ci_u = RD_pooled + 1.96#se_rd_p;
print rdi RDi_ci_l RDi_ci_u RD_pooled RDp_ci_l RDp_ci_u;
/* rd顯著性 h0:rd=0 if |zrd|>1.96 reject ho*/
Zrd = RD_pooled/se_rd_p;
zrd_pv = (1-probnorm(abs(Zrd)))#2;
print zrd zrd_pv;
/* rd同質性 h0:rd1=rd2...*/
Qrd = sum(wi#((RDi-RD_pooled)##2));
Qrd_pv = (1-probchi(Qrd, &obs_1.));
print Qrd Qrd_pv;
i2 = (qrd - &obs_1.) / qrd;
print i2;
/* Random Effect */
t_2 = max(0, (Qrd-(&obs_1.)) / (sum(wi)- (sum(wi##2)) /sum(wi))) ;
wi_ = 1/(var_rdi+t_2);
theta_ = sum(wi#rdi) / sum(rdi);
se_theta_ = 1 / sqrt(sum(wi_));
theta__ci_l = theta_ - 1.96 # se_theta_;
theta__ci_u = theta_ + 1.96 # se_theta_;
r_pooled_rd = exp(sum(wi_#rdi)/sum(wi_));
se_r_pooled_rd = 1 / sqrt(sum(wi_));
r_pooled_rd_ci_l = r_pooled_rd - 1.96#se_r_pooled_rd;
r_pooled_rd_ci_u = r_pooled_rd + 1.96#se_r_pooled_rd;
print t_2 wi_ theta_ se_theta_ theta__ci_l theta__ci_u ;
print r_pooled_rd se_r_pooled_rd r_pooled_rd_ci_l r_pooled_rd_ci_u;
/*theta_的有效性檢驗， h0 : theta_ =0 */
utheta_ = ((sum(wi_ # rdi))##2 ) / sum(wi_);
utheta__pv = (1-probchi(utheta_, 1));
zutheta_ = theta_ / se_theta_;
zutheta__pv = (1-probnorm(abs(zutheta_)))#2;
print utheta_ utheta__pv zutheta_ zutheta__pv;
median_rdi = median(rdi);
wi_pct = wi/sum(wi);
create output var {
i2 Qw_pv rdi RDi_ci_l RDi_ci_u RD_pooled RDp_ci_l RDp_ci_u se_rd zrd Qrd Qrd_pv
wi median_rdi n1 n2 wi_pct};
append;
close output;
GroupId=3;
create pooled var {RD_pooled RDp_ci_l RDp_ci_u GroupId};
append;
close output;
create ran_pooled var { r_pooled_rd r_pooled_rd_ci_l r_pooled_rd_ci_u GroupId};
append;
close output;
quit;title;
data output2;
merge output new_data;
obs = _n_;
GroupId = 1;
run;
data output3;
set output2
pooled (rename=(RD_pooled=rdi RDp_ci_l=rdi_ci_l RDp_ci_u=rdi_ci_u))
ran_pooled (rename=(r_pooled_rd=rdi r_pooled_rd_ci_l=rdi_ci_l r_pooled_rd_ci_u=rdi_ci_u));
if _n_ = &obs_.+1 then &references. = "Pooled_F";
if _n_ = &obs_.+2 then &references. = "Pooled_R";
run;
proc sgplot data = output2;
HIGHLOW y=obs HIGH=rdi_ci_u LOW=rdi_ci_l / close = rdi;
REFLINE 0 / label = "0" AXIS= X;
yaxis REVERSE;
quit;
proc print data = output2(keep = &references. rdi rdi_ci_u rdi_ci_l);
run;
/*
proc sgplot data = output2;
STEP x = obs y= wi / DATALABEL = wi;
quit;
*/
proc sgplot data = output2;
scatter y = se_rd x = rdi / datalabel = obs;
yaxis REVERSE;
refline median_rdi / axis = x;
quit;
%mend;
%macro my_meta_norm(data = , NT = , YT = , ST = , NC = , YC = , SC = , references=);
%let open_d = %sysfunc (open(&data.));
%let obs_ = %sysfunc(attrn(&open_d., nobs));
%let close_d = %sysfunc(close(&open_d.));
%let obs_1 = %sysevalf(&obs_ - 1);
proc iml;
use &data.;
read all var _all_;
close &data.;
/* si = sqrt(st##2+sc##2-2*0.7*st*sc) /sqrt(2*(1-0.7)); */
/* si = sqrt(((nt-1)#((st)##2) + (nc-1)#((sc)##2))/(nt+nc-2));*/
/* esi = (yt - yc)/si;*/
/* proc iml;
yt = {54 48 40};
yc = {71 72 75};
st={14 15 11};
sc = {6 8 8};
n = {36 39 37};
nt = {36 39 37};
nc = {36 39 37};
si = sqrt(st##2+sc##2-2#0.7#st#sc) /sqrt(2#(1-0.7));
esi = (yt - yc)/si;
se_esi = sqrt((1/n)+(esi##2)/(2#n)) # sqrt(2#(1-0.7));
create d1 var {yt yc st sc n si esi se_esi nt nc};
append;
close d1;
quit;*/
/* m = nt+nc -2;*/
/* m = nt -1;*/
/* cm = 1 - (3/(4#m-1));*/
/* esui = cm#esi;*/
/* var_esui = (nt+nc)/(nt#nc) + ((esui)##2)/(2#(nt+nc));*/
/* se_esui = sqrt(var_esui);*/
esui_ci_l = esui - 1.96 #se_esui;
esui_ci_u = esui + 1.96 #se_esui;
print si esi esui var_esui esui_ci_l esui_ci_u;
wi = 1/var_esui;
pooled_es = sum(wi#esui) / sum(wi);
se_pooled_es = 1 / sqrt(sum(wi));
pooled_es_ci_l = pooled_es - 1.96#se_pooled_es;
pooled_es_ci_u = pooled_es + 1.96#se_pooled_es;
print wi pooled_es se_pooled_es pooled_es_ci_l pooled_es_ci_u;
/* print wi2 pooled_es2 se_pooled_es2;*/
/*ES的顯著性檢驗, */
/* i.e. Testing H0: Effect Size = 0*/
ues = ((sum(wi#esui))##2) / sum(wi);
ues_pv = (1-probchi(ues, 1));
print ues ues_pv;
median_esui = median(esui);
/*ES的同質性檢驗, */
/* i.e. Testing H0: = = = */
qes = sum(wi# ((esui-pooled_es)##2));
qes_pv = (1-probchi(qes, &obs_1.));
print qes qes_pv;
i2 = (qes -&obs_1.)/qes;
print i2;
wi_pct = wi/sum(wi);
/* Random Effect Model*/
t_2 = max(0, (qes-(&obs_1.)) / (sum(wi)- (sum(wi##2)) /sum(wi))) ;
wi_ = 1/(var_esui+t_2);
wi__ = (wi_/sum(wi_) )*100;
theta_ = sum(wi_#esui) / sum(wi_);
se_theta_ = 1 / sqrt(sum(wi_));
theta__ci_l = theta_ - 1.96 # se_theta_;
theta__ci_u = theta_ + 1.96 # se_theta_;
r_pooled_es = sum(wi_#esui) / sum(wi_);
se_r_pooled_es = 1 / sqrt(sum(wi_));
r_pooled_es_ci_l = pooled_es - 1.96#se_r_pooled_es;
r_pooled_es_ci_u = pooled_es + 1.96#se_r_pooled_es;
print t_2 wi_ wi__ theta_ se_theta_ theta__ci_l theta__ci_u ;
print r_pooled_es se_r_pooled_es r_pooled_es_ci_l r_pooled_es_ci_u;
/*r ES的顯著性檢驗, */
/* i.e. Testing H0: Effect Size = 0*/
ues = ((sum(wi_#esui))##2) / sum(wi);
ues_pv = (1-probchi(ues, 1));
print ues ues_pv;
/*r ES的同質性檢驗, */
/* i.e. Testing H0: = = = */
qes = sum(wi_# ((esui-pooled_es)##2));
qes_pv = (1-probchi(qes, &obs_1.));
print qes qes_pv;
i2 = (qes -&obs_1.)/qes;
print i2;
wi_pct = wi/sum(wi);
/*theta_的有效性檢驗， h0 : theta_ =0 */
utheta_ = ((sum(wi_ # esui))##2 ) / sum(wi_);
utheta__pv = (1-probchi(utheta_, 1));
zutheta_ = theta_ / se_theta_;
zutheta__pv = (1-probnorm(abs(zutheta_)))#2;
print utheta_ utheta__pv zutheta_ zutheta__pv;
GroupId=3;
create pooled var {pooled_es pooled_es_ci_l pooled_es_ci_u GroupId};
append;
close pooled;
create ran_pooled var {r_pooled_es r_pooled_es_ci_l r_pooled_es_ci_u GroupId};
append;
close ran_pooled;
;
create output var { i2 si esi esui var_esui se_esui esui_ci_l esui_ci_u pooled_es se_pooled_es pooled_es_ci_l pooled_es_ci_u ues ues_pv qes qes_pv
median_esui wi_ wi__ r_pooled_es se_r_pooled_es r_pooled_es_ci_l r_pooled_es_ci_u };
append;
close output;
quit;
data output2;
merge output (rename=(esui=smd esui_CI_L=lcl esui_CI_u=ucl WI__ = weigth_pct)) &data.;
obs = _n_;
GroupId = 1;
smd_ = Round(smd, 0.01);
run;/*
data output3;
set output2
pooled (rename=(pooled_es=esui pooled_es_ci_l=esui_ci_l pooled_es_ci_u=esui_ci_u))
ran_pooled (rename=(r_pooled_es=esui r_pooled_es_ci_l=esui_ci_l r_pooled_es_ci_u=esui_ci_u));
if _n_ = &obs_.+1 then &references. = "Pooled_F";
if _n_ = &obs_.+2 then &references. = "Pooled_R";
run;
*/
proc sgplot data = output2;
HIGHLOW y=obs HIGH=esui_ci_u LOW=esui_ci_l / close = esui ;
REFLINE 0 / label = "0" AXIS= X;
yaxis REVERSE;
quit;
proc print data = output2(keep = &references. esui esui_ci_l esui_ci_u wi);
run;
proc sgplot data = output2;
scatter y = se_esui x = esui / datalabel = obs;
refline median_esui / axis = x;
yaxis REVERSE;
quit;
%mend;
%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;
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);
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;
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;
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;
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;
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;
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;
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));
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)));
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)));
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));
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;
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';
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'-----------'/
@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'--------------------------------------------------------------------------------------'//;
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;*/
title;
quit;
%Mend PubBias;
%macro ForestMacro (
Data=, /*--Data Set Name (Required)--*/
Study=, /*--Variable name for Study (Required)--*/
OddsRatio=, /*--Variable name for Odds Ratio (Required)--*/
LCL=, /*--Variable name for Lower Confidence Limit (Required)--*/
UCL=, /*--Variable name for Upper Confidence Limit (Required)--*/
Group=GroupId, /*--Variable name for Study Type--*/
Weight=, /*--Variable name for Study Weight in %--*/
StatCol1=, /*--Variable name for Stat Column 1--*/
StatCol2=, /*--Variable name for Stat Column 2--*/
StatCol3=, /*--Variable name for Stat Column 3--*/
StatCol4=, /*--Variable name for Stat Column 4--*/
DisplayCols=YES, /*--Display the columns for OR, LCL, UCL & Weight--*/
WtFactor=, /*--Multiplier factor for Study Weights--*/
/*--If not provided WtFactor is computed internally--*/
Bands=YES, /*--Draw Horizontal Alternating Bands--*/
Borders=NO, /*--Draw Borders--*/
GraphWalls=NO, /*--Draw Filled Walls behind the Graph--*/
StatWalls=NO, /*--Draw Filled Walls behind the Statistics Tables--*/
Width=6.4in, /*--Default width of the graph in pixels--*/
Height=, /*--Default height of the graph is computed based on number of observations--*/
LabelColWidth=0.2, /*--Fractional width for Label Column--*/
Label1=Favors Treatment, /*--Favorable Label--*/
Label2=Favors Control, /*--Unfavorable Label--*/
PlotTitle=Odds Ratio and 95% CL, /*--Plot Title--*/
FootNote=, /*--Graph Footnote--*/
Title2=, /*--Graph title2--*/
Title=Impact of Treatment on Mortality by Study /*--Graph Title--*/
);
%local WeightVar MarkerSize GraphColWidth StatColWidth Border DisplaySecondary GraphWallDisplay StatWallDisplay;
%local OddsLabel LowerLabel UpperLabel WeightLabel SLabel1 SLabel2 SLabel3 SLabel4;
%local GraphHeight Ratio RowHeight HeaderHeight Nobs;
/*--Data, Study, OddsRatio, LCL and UCL are required --*/
/*--Group is optional --*/
/*--Terminatethese required parameters are not supplied--*/
%if %length(&Data) eq 0 %then %do;
%put The parameter 'Data' is required - Forest Macro Terminated.;
%goto finished;
%end;
%else %if %length(&Study) eq 0 %then %do;
%put The parameter 'Study' is required - Forest Macro Terminated.;
%goto finished;
%end;
%else %if %length(&LCL) eq 0 %then %do;
%put The parameter 'LCL' is required - Forest Macro Terminated.;
%goto finished;
%end;
%else %if %length(&UCL) eq 0 %then %do;
%put The parameter 'UCL' is required - Forest Macro Terminated.;
%goto finished;
%end;
%else %if %length(&OddsRatio) eq 0 %then %do;
%put The parameter 'OddsRatio' is required - Forest Macro Terminated.;
%goto finished;
%end;
/*--Initialize GraphHeight, Height per row and Height for other graph items--*/
%let GraphHeight=&Height;
%let RowHeight=22;
%let HeaderHeight=100;
%if %length(&Footnote) ne 0 %then %do;
%let HeaderHeight=115;
%end;
/*--If the Weight column is not provided, use equal weights, and suppress display of Weight stat--*/
%if &Weight eq %then %do;
%let WeightVar = _Weight;
%let MarkerSize = 7;
%end;
%else %do;
%let WeightVar=&Weight;
%let MarkerSize = 0;
%end;
/*--Set up GTL options for borders--*/
%let DisplaySecondary = displaysecondary=none;
%let Borders=%upcase(&Borders);
%if &Borders eq YES or &Borders eq Y %then %do;
%let Border = line;
%let DisplaySecondary = displaysecondary=(line);
%end;
/*--Set up GTL options for GraphWall Display--*/
%let GraphWallDisplay = walldisplay=none;
%let GraphWalls=%upcase(&GraphWalls);
%if &GraphWalls eq YES or &GraphWalls eq Y %then %do;
%let GraphWallDisplay = walldisplay=(fill);
%end;
/*--Set up GTL options for StatWall Display--*/
%let StatWallDisplay = walldisplay=none;
%let StatWalls=%upcase(&StatWalls);
%if &StatWalls eq YES or &StatWalls eq Y %then %do;
%let StatWallDisplay = walldisplay=(fill);
%end;
/*--Create Label Columns for standard and additional columns--*/
/*--Load Stat Column Label or name into macro for label column value--*/
%let dsid=%sysfunc(open(&Data));
%if &dsid %then %do;
%let Nobs=%sysfunc(attrn(&dsid, nlobs));
%if &Nobs eq 0 %then %do;
%put The Data Set &Data has no observations - Forest Macro Terminated.;
%let rc=%sysfunc(close(&dsid));
%goto finished;
%end;
%if &Nobs gt 150 %then %do;
%put The Data Set &Data has over 100 observations - Forest Macro Terminated.;
%let rc=%sysfunc(close(&dsid));
%goto finished;
%end;
/*--Count the number of stat columns--*/
%let idx=0;
/*--Column display information for the OddsRatio column--*/
%let DisplayCols=%upcase(&DisplayCols);
%if &DisplayCols eq YES or &DisplayCols eq Y %then %do;
%let OddsLabel=%sysfunc(varlabel(&dsid, %sysfunc(varnum(&dsid,&OddsRatio))));
%if %length(&OddsLabel) eq 0 %then %let OddsLabel=&OddsRatio;
%let idx= %eval(&idx+1);
%let LowerLabel=%sysfunc(varlabel(&dsid, %sysfunc(varnum(&dsid,&LCL))));
%if %length(&LowerLabel) eq 0 %then %let LowerLabel=&LCL;
%let idx= %eval(&idx+1);
%let UpperLabel=%sysfunc(varlabel(&dsid, %sysfunc(varnum(&dsid,&UCL))));
%if %length(&UpperLabel) eq 0 %then %let UpperLabel=&UCL;
%let idx= %eval(&idx+1);
%if &Weight ne %then %do;
%let WeightLabel=%sysfunc(varlabel(&dsid, %sysfunc(varnum(&dsid,&Weight))));
%if %length(&WeightLabel) eq 0 %then %let WeightLabel=&Weight;
%let idx= %eval(&idx+1);
%end;
%end;
/*--Additional columns to be displayed--*/
%if %length(&StatCol1) ne 0 %then %do;
%let SLabel1=%sysfunc(varlabel(&dsid, %sysfunc(varnum(&dsid,&StatCol1))));
%if %length(&SLabel1) eq 0 %then %let SLabel1=&StatCol1;
%let idx= %eval(&idx+1);
%end;
%if %length(&StatCol2) ne 0 %then %do;
%let SLabel2=%sysfunc(varlabel(&dsid, %sysfunc(varnum(&dsid,&StatCol2))));
%if %length(&SLabel2) eq 0 %then %let SLabel2=&StatCol2;
%let idx= %eval(&idx+1);
%end;
%if %length(&StatCol3) ne 0 %then %do;
%let SLabel3=%sysfunc(varlabel(&dsid, %sysfunc(varnum(&dsid,&StatCol3))));
%if %length(&SLabel3) eq 0 %then %let SLabel3=&StatCol3;
%let idx= %eval(&idx+1);
%end;
%if %length(&StatCol4) ne 0 %then %do;
%let SLabel4=%sysfunc(varlabel(&dsid, %sysfunc(varnum(&dsid,&StatCol4))));
%if %length(&SLabel4) eq 0 %then %let SLabel4=&StatCol4;
%let idx= %eval(&idx+1);
%end;
%let rc=%sysfunc(close(&dsid));
/*--Set column weights based on number of stat columns--*/
%let StatColWidth=%sysevalf(&idx * 0.075);
%let GraphColWidth= %sysevalf(1.0 - &LabelColWidth - &StatColWidth);
%end;
%else %do;
%put The data set &Data does not exist - Forest Macro Terminated.;
%goto finished;
%end;
/*--Compute Weight Factor if not provided --*/
/*--Estimate height of graph if not provided--*/
data _null_;
set &Data end=last;
retain totalweight 0;
totalweight+wt;
if last then do;
%if &wtFactor eq %then %do;
if totalweight <= 0 then totalweight=1;
call symput ('wtFactor', 1 / totalweight);
%end;
/*--Estimate Ratio of Plot height by Graph Height--*/
call symput ('Ratio', (_N_* &RowHeight)/(_N_* &RowHeight + &HeaderHeight));
/*--Estimate the optimal height of the graph based on obs count--*/
%if &Height eq %then %do;
call symput ('GraphHeight', _N_ * &RowHeight + &HeaderHeight);
%end;
end;
run;
/*--Append a PX only if this internally estimated--*/
%if &Height eq %then %do;
%let GraphHeight=&GraphHeight.px;
%end;
/*--Process Data--.3*/
data _forest;
set &Data;
format _wt ;
_ObsId=_N_;
%if &Weight eq %then %do;
&WeightVar=0;
%end;
label _wt=&WeightLabel;
/*--If Group column is provided--*/
%if %length(&Group) ne 0 %then %do;
/*--Group=1 (Study) values will be drawn without a group role--*/
if &group=1 then do;
_wt=&WeightVar / 100;
_grp=10;
_1 = &OddsRatio;
_lcl1=&LCL;
_ucl1=&UCL;
/*--Compute marker width--*/
_x1=&OddsRatio / (10 ** (&WeightVar*&WtFactor/2));
_x2=&OddsRatio * (10 ** (&WeightVar*&WtFactor/2));
end;
/*--Group=2 & 3 (SubGroup and Overall) values will be drawn with groupindex=2 & 3--*/
else if &group > 1 then do;
_grp=&group;
_2 = &OddsRatio;
end;
%end;
%else %do;
_wt=&WeightVar / 100;
_grp=10;
_1 = &OddsRatio;
_lcl1=&LCL;
_ucl1=&UCL;
/*--Compute marker width--*/
_x1=&OddsRatio / (10 ** (&WeightVar*&WtFactor/2));
_x2=&OddsRatio * (10 ** (&WeightVar*&WtFactor/2));
%end;
/*--Create label columns for standard and additional statistic--*/
%if %length(&OddsRatio) ne 0 %then %do;
_OddsRatioLabel = symget('OddsLabel');
%end;
%if %length(&LCL) ne 0 %then %do;
_LowerLabel = symget('LowerLabel');
%end;
%if %length(&UCL) ne 0 %then %do;
_UpperLabel = symget('UpperLabel');
%end;
%if %length(&Weight) ne 0 %then %do;
_WeightLabel = symget('WeightLabel');
%end;
%if %length(&StatCol1) ne 0 %then %do;
_StatColLabel1 = symget('SLabel1');
_StatCol1 = &StatCol1;
%end;
%if %length(&StatCol2) ne 0 %then %do;
_StatColLabel2 = symget('SLabel2');
_StatCol2 = &StatCol2;
%end;
%if %length(&StatCol3) ne 0 %then %do;
_StatColLabel3 = symget('SLabel3');
_StatCol3 = &StatCol3;
%end;
%if %length(&StatCol4) ne 0 %then %do;
_StatColLabel4 = symget('SLabel4');
_StatCol4 = &StatCol4;
%end;
run;
/*--Reverse the order to avoid putting axis reverse--*/
proc sort data=_forest out=_forest;
by descending _ObsId;
run;
/*--Add sequence numbers to each observation--*/
data _forest;
set _forest;
studyvalue=_n_;
run;
/*--Output values and formatted strings to data set--*/
data _forestFormat;
set _forest end=last;
keep label start end fmtname type hlo;
retain fmtname '_Study' type 'n';
label=&Study;
start=studyvalue;
end=studyvalue;
output;
if last then do;
hlo='O';
label='Other';
output;
end;
run;
/*--Create Format from data set--*/
proc format library=work cntlin=_forestFormat;
run;
/*--Apply format to study values--*/
/*--Compute width of box proportional to weight in log scale--*/
data _forest;
format studyvalue _study.;
set _forest;
%let Bands=%upcase(&Bands);
%if &Bands eq YES or &Bands eq Y %then %do;
if mod(studyvalue, 2) = 0 then _StudyRef=StudyValue;
%end;
run;
/*--Compute top and bottom offsets--*/
data _null_;
pct=&Ratio/nobs;
thk=pct* 0.9 *100;
call symputx("pct", pct);
call symputx("pct2", 2*pct);
call symputx("RefThickness", thk);
call symputx("count", nobs);
set _forest nobs=nobs;
run;
/*title;*/
/*options nodate nonumber;*/
/*--Define GTL template for graph--*/
proc template;
define statgraph ForestMacro;
begingraph / designwidth=&Width designheight=&GraphHeight;
entrytitle "&Title";
entryfootnote halign=left "&FootNote";
%if %length(&title2) ne 0 %then %do;
entrytitle "&title2" / textattrs=graphLabelText;
%end;
layout lattice / columns=3 columnweights=(&LabelColWidth &GraphColWidth &StatColWidth) columngutter=0
rowdatarange=union;
/*--Column # 1 contains the Study Labels using Secondary Y axis--*/
layout overlay / walldisplay=none x2axisopts=(display=none)
yaxisopts=(linearopts=(tickvaluesequence=(start=1 end=&count increment=1))
offsetmin=&pct2 offsetmax=&pct display=none
displaysecondary=(tickvalues &border));
scatterplot y=studyvalue x=_1 / yaxis=Y xaxis=X2 markerattrs=(size=0) includemissinggroup=true;
scatterplot y=studyvalue x=_1 / yaxis=Y xaxis=X2 markerattrs=(size=0) includemissinggroup=true;
endlayout;
/*--Column # 2 contains the graph--*/
layout overlay / &GraphWallDisplay border=false
yaxisopts=(linearopts=(tickvaluesequence=(start=1 end=&count increment=1))
offsetmin=&pct2 offsetmax=&pct display=none);
/*--Draw alternating bands using referenceline--*/
%if &Bands eq YES or &Bands eq Y %then %do;
referenceline y=_StudyRef / lineattrs=(thickness=&RefThickness.PCT) datatransparency=0.9;
%end;
/*--Draw Markers for SubGroup and Overall values--*/
%if %length(&Group) ne 0 %then %do;
scatterplot y=studyvalue x=_2 / markerattrs=(symbol=diamondfilled size=10) group=_grp
includemissinggroup=true index=_grp;
%end;
/*--Draw OddsRatio and Limits for Study Values--*/
scatterplot y=studyvalue x=_1 / xerrorupper=_ucl1 xerrorlower=_lcl1
markerattrs=graphdata1(symbol=squarefilled size=&MarkerSize);
/*--Draw box representing the weight of the study--*/
vectorplot y=studyvalue x=_x2 xorigin=_x1 yorigin=studyvalue / lineattrs=GraphData1(thickness=8)
arrowheads=false;
/*--Draw Reference lines and labels--*/
referenceline x=0;
entry halign=left "&Label1" / valign=bottom;
entry halign=right "&Label2" / valign=bottom;
endlayout;
/*--Column # 2 contains the statistics data--*/
layout overlay / &StatWallDisplay border=false
x2axisopts=(display=(tickvalues &border) displaysecondary=(line))
yaxisopts=(linearopts=(tickvaluesequence=(start=1 end=&count increment=1))
offsetmin=&pct2 offsetmax=&pct
display=none &DisplaySecondary.);
/*--Draw alternating bands using referenceline--*/
%if &Bands eq YES %then %do;
referenceline y=_StudyRef / lineattrs=(thickness=&RefThickness.PCT) datatransparency=0.9;
%end;
/*--Draw standard statistics columns--*/
%if &DisplayCols eq YES or &DisplayCols eq Y %then %do;
scatterplot y=studyvalue x=_OddsRatioLabel / markercharacter=&OddsRatio xaxis=x2;
scatterplot y=studyvalue x=_LowerLabel / markercharacter=&LCL xaxis=x2;
scatterplot y=studyvalue x=_UpperLabel / markercharacter=&UCL xaxis=x2;
%if &Weight ne %then %do;
scatterplot y=studyvalue x=_WeightLabel / markercharacter=_wt xaxis=x2;
%end;
%end;
/*--Draw additional statistics columns--*/
%if %length(&StatCol1) ne 0 %then %do;
scatterplot y=studyvalue x=_StatColLabel1 / markercharacter=&StatCol1 xaxis=x2;
%end;
%if %length(&StatCol2) ne 0 %then %do;
scatterplot y=studyvalue x=_StatColLabel2 / markercharacter=&StatCol2 xaxis=x2;
%end;
%if %length(&StatCol3) ne 0 %then %do;
scatterplot y=studyvalue x=_StatColLabel3 / markercharacter=&StatCol3 xaxis=x2;
%end;
%if %length(&StatCol4) ne 0 %then %do;
scatterplot y=studyvalue x=_StatColLabel4 / markercharacter=&StatCol4 xaxis=x2;
%end;
endlayout;
endlayout;
endgraph;
end;
run;
proc sgrender data=_forest template=ForestMacro description='Forest Plot';
run;
%finished:
title;
FootNote;
%mend ForestMacro;
0 意見:
張貼留言