/*定义纳入试验的个数*/
%let trialnum=4;
/*创建数据集*/
data origin;
input id trial study$ time treat n mean sd @@;
cards;
1 1 Gupta 0 0 50 5.2 4.4
2 1 Gupta 1 0 50 23.1 1.2
3 1 Gupta 2 0 50 25.1 1.1
4 1 Gupta 0 1 50 4.5 4.7
5 1 Gupta 1 1 50 20.7 1.3
6 1 Gupta 2 1 50 23.7 1.6
7 2 Tan 0 0 30 8.4 0.5
8 2 Tan 1 0 30 26.4 1.8
9 2 Tan 2 0 30 21.8 2.1
10 2 Tan 0 1 30 8.3 0.4
11 2 Tan 1 1 30 20.8 2.3
12 2 Tan 2 1 30 18.4 2.8
13 3 Kuntz 0 0 100 4.9 3.8
14 3 Kuntz 1 0 100 25.1 6.9
15 3 Kuntz 2 0 100 27.9 9.9
16 3 Kuntz 0 1 100 5.9 3.9
17 3 Kuntz 1 1 100 25.1 9.4
18 3 Kuntz 2 1 100 27.7 12.2
19 4 Montorsi 0 0 52 8.2 3.2
20 4 Montorsi 1 0 52 23.1 8.6
21 4 Montorsi 2 0 52 25.1 7.2
22 4 Montorsi 0 1 48 7.8 3.6
23 4 Montorsi 1 1 48 26.5 15.5
24 4 Montorsi 2 1 48 24.7 10
;
run;
/*Proc print data=origin;run;代码用于查看创建的数据集是否正确*/
/*对每一个试验进行回归系数的估计*/
%macro reg(isize);
%do i=1 %to &isize.;
data origin&i.;
set origin;
if trial = &i.;*提取试验;
w=1/(sd/sqrt(n));*计算权重;
timetrt=time*treat;*产生时间和治疗交互作用的因子;
run;
proc reg data= origin&i.;
model mean=treat time timetrt;
weight w;
ods output ParameterEstimates=PEst&i.;
run;
%end;
%mend; /*宏程序结束*/
%reg(&trialnum);
data PEst;
set PEst1-PEst4;/*合并参数*/
run;
data PEst;
set PEst;
if Variable = "Intercept" then do; var=1; end;
else if Variable ="treat" then do; var=2; end;
else if Variable ="time" then do; var=3; end;
else do;var=4; end;
keep Variable var Estimate StdErr;
run;
proc print data=PEst noobs; run;
data start;/*设置初始值*/
input est;
cards;
0.01
run;
/*对每一个回归系数进行meta分析*/
%macro meta(isize);
%do i=1 %to &isize.;
data PEst&i.;
set PEst;
if var = &i.;
run;
data PEst&i.;
set PEst&i.;
study=_N_;
run;
proc print data=PEst&i. noobs; run;
data covvars&i.;
set PEst&i.;
est=StdErr**2;
keep est;
run;
data covvars&i.;
set start covvars&i.;
run;
proc print data=covvars&i. noobs; run;
ODS TRACE ON;
proc mixed data=PEst&i. method = ml cl;
class study;
model Estimate=/ s cl;
random int /subject = study;
repeated/group=study;
parms/parmsdata= covvars&i. eqcons=2 to 5;
ods output SolutionF=solF&i.;
run;
ODS TRACE OFF;
%end;
%mend; /*宏程序结束*/
%meta(4);/*此处的4代表回归系数的个数*/
/*合并回归系数的结果，并打印显示*/
data solF;
set solF1-solF4;
run;
proc print data=solF noobs; run;
%macro forestdata(isize);
%do i=1 %to &isize.;
data PEst&i.;
set PEst&i.;
keep Estimate StdErr study;
run;
data solF&i.;
set solF&i.;
keep Estimate StdErr Probt;
run;
data forest&i.(rename=(study=studyNUM));
set PEst&i. solF&i.;
run;
proc print data=forest&i. noobs; run;
%end;
%mend; /*宏程序结束*/
%forestdata(4);/*此处的4代表回归系数的个数*/
%macro Anno_data(isize);
%do i=1 %to &isize.;
data forest&i.;
set forest&i.;
if studyNUM=. then do;wi=0;end;
else do;wi=1/(StdErr)**2;end;
run;
proc sql;
create table forest&i. as
select StudyNUM,Estimate,StdErr,wi,Probt, sum(wi) as sumwi
from forest&i.;
quit;
data forest&i.(rename=(Estimate=OddsRatio));
set forest&i.;
format weight percent5.;
Estimate=round(Estimate,.01);
LowerCL=round(Estimate-1.96*StdErr,.01);
UpperCL=round(Estimate+1.96*StdErr,.01);
Weight=round(wi/sumwi,.0001);
if Weight=0 then do;Weight=.; end;
run;
data forest&i.;
set forest&i.;
length Study $ 16;
length studyref $ 16;
format Q1 Q3 4.2;
if studyNUM=1 then do;Study="Gupta";end;
else if studyNUM=2 then do;Study="Tan";end;
else if studyNUM=3 then do;Study="Kuntz";end;
else if studyNUM=4 then do;Study="Montorsi";end;
else do;Study="";end;
if Study="" then do; grp=2;studyref="Overall";study2="Overall";end;
else do;grp=1;end;
ObsId=_N_;
if mod (obsid, 2) = 0 then studyref=Study;
lcl2=lowercl;
ucl2=uppercl;
Q1=min(round(OddsRatio-weight,.01),round(OddsRatio+weight,.01));
Q3=max(round(OddsRatio-weight,.01),round(OddsRatio+weight,.01));
ES='ES'; LCL='LCL'; UCL='UCL'; WT='Weight';Pvalue='Pvalue';
drop wi studyNUM StdErr sumwi;
run;
/**/
proc print data=forest&i.;run;
%end;
%mend; /*宏程序结束*/
%Anno_data(4);/*此处的4代表回归系数的个数*/
%macro Meta_forest(isize);
%do i=1 %to &isize.;
data forest;
set forest&i.;
run;
data _null_;
pct=0.75/nobs;
height=3;
dpi=100;
call symputx("pct", pct);
call symputx("pct2", 2*pct);
call symputx("dpi", dpi);
call symputx("height", height);
call symputx("heightin", height || "in");
call symputx("thickness", floor(height*dpi*pct));
set forest nobs=nobs;
run;
/*--Create custom style--*/
/*--创建自定义样式--*/
proc template;
define style styles.ForestColor93;
parent = Styles.htmlBlue;
style GraphFonts from GraphFonts /
'GraphDataFont' = ("<sans-serif>, <MTsans-serif>",7pt)
'GraphValueFont' = ("<sans-serif>, <MTsans-serif>",7pt);
style GraphData1 from GraphData1 /
contrastcolor = GraphColors('gcdata2')
color = GraphColors('gdata2');
style GraphData2 from GraphData2 /
contrastcolor = GraphColors('gcdata1')
color = GraphColors('gdata1');
end;
run;
/*--Set Style, DPI, image parameters and titles--*/
/*--设置样式，DPI，图像参数和标题--*/
ods html off;/*ods html close可以只生成图片到文件夹;*/
ods graphics / reset width=6in height=&heightin imagename="ForestPlotColor";
ods listing sge=off style=Styles.ForestColor93 image_dpi=&dpi;
title h=12pt "Impact of Covariates on Survival by Study";
title2 h=8pt 'The estimated regression coefficients and 95% CL';
/*--Create the plot--*/
proc sgplot data=forest noautolegend nocycleattrs;
/*--Draw alternate reference line用于绘制研究与研究之间间隔显色--*/
refline studyref / lineattrs=(thickness=&thickness) transparency=0.85;/*0.85设置颜色的深浅程度，越小颜色越深，1为透明*/
/*--Display overall values (Study2) using scatter plot使用散点图显示总体值（Study2）--*/
scatter y=study2 x=oddsratio / markerattrs=(symbol=diamondfilled size=10);
/*--Display individual values (Study) using highLow plot使用highLow图显示单个值（Study）--*/
highlow y=study low=lcl2 high=ucl2 / type=line;
highlow y=study low=q1 high=q3 / type=bar;
/*--Display statistics columns on X2 axis在X2轴上显示统计列--*/
scatter y=study x=ES / markerchar=oddsratio x2axis;
scatter y=study x=lcl / markerchar=lowercl x2axis;
scatter y=study x=ucl / markerchar=uppercl x2axis;
scatter y=study x=wt / markerchar=weight x2axis;
/*--Display statistics columns on X2 axis在X2轴上显示统计列--*/
scatter y=study2 x=ES / markerchar=oddsratio x2axis;
scatter y=study2 x=lcl / markerchar=lowercl x2axis;
scatter y=study2 x=ucl / markerchar=uppercl x2axis;
scatter y=study2 x=wt / markerchar=weight x2axis;
scatter y=study2 x=Pvalue / markerchar=Probt x2axis;
/*--Draw other details in the graph--*/
refline 0 100 / axis=x;
refline -10 10 / axis=x lineattrs=(pattern=shortdash) transparency=0.5;/*设置参考线，transparency透明度*/
inset ' Unfavorable factor' / position=bottomleft;
inset 'Favorable factor' / position=bottom;
/*--Set X, X2 axis properties with fixed offsets--*/
xaxis type=log offsetmin=0 offsetmax=0.35 min=-30 max=20 minor display=(nolabel) ;
x2axis offsetmin=0.7 display=(noticks nolabel);
/*--Set Y axis properties (including reverse) using offsets computed earlier--*/
yaxis display=(noticks nolabel) offsetmin=&pct offsetmax=&pct2 reverse;
run;
%end;
%mend; /*宏程序结束*/
%Meta_forest(1);/*此处的4代表回归系数的个数*/
%Meta_forest(2);/*此处的4代表回归系数的个数*/
扫码加好友，拉您进群
全部版块
我的主页
收藏
