ERROR: Unable to transcode data to/from UCS-2 encoding.

可試試此MACRO(連結)
或使用KCOUNT() 找中文 檢查長度夠不夠

以下舊文
















append失敗時，可能為以下情況:

1.非NULL欄位，觀測值為NULL。
2.文字欄位長度不夠，產生亂碼。
    (數值欄位也有可能，SAS大多為8，但如append到SQL DB就須注意)
3.亂碼。

不過可以先去檢查程式碼內是否有使用  TEST_VAR =  "內有中文" ，之類的引號內有中文，
要將編碼改成 UTF-8，下圖方式即可修改編碼: 將程式碼另存選擇編碼覆蓋即可!

可使用KCOUNT() 找中文

如是資料原本身就給出亂碼的話，
可試試此MACRO(連結)

以下MACRO為一筆一筆的APPEND，看哪幾筆資料會發生ERROR，再去檢查。

 

%macro Debug_append(base=,
   data= );
%let open_d = %sysfunc (open(&data.)); 
%let  obs = %sysfunc(attrn(&open_d., nobs));
%let close_d = %sysfunc(close(&open_d.));
 
%do _n = 1  %to &obs.;
 data _data_n&_n.;
  set &data.;
  if _n_ = &_n.;
 run;

 %put ==================================================================;
 %put 第 &_n. 筆觀測值異常;
 %put 1. 非NULL欄位   為  NULL;
 %put 2.文字欄位長度不夠 產生亂碼 ;
 %put 3.亂碼;
 %put ==================================================================;

 %let open_d = %sysfunc (open(_data_n&_n.)); 
 %let  obs2 = %sysfunc(attrn(&open_d., nobs));
 %let close_d = %sysfunc(close(&open_d.));
 %if &obs2. > 0 %then %do;
  proc append base=&base.  data=_data_n&_n.   force;
  run;
 %end;
%end;

%mend;

%Debug_append(base=db_core.FSC_ENTITY_WATCH_LIST_DIM_Jimmy,
   data=dailprep._0304_ETLS_NEWRECORDS22222 )

SOURCE_CHECK_2

SAS AML 的 SOURCE_CHECK 已經有資料筆數、NULL、UNIQUE與LOOKUP的檢核了

SOURCE_CHECK_2 檢核觀測值的MACRO。

範例:
data DB_A;
input Pty $ score  number;
datalines;
A 100 1
B  90 2
C  80 4
C .  .
D  50 5
F  72 6
G . .
;
run;

data DB_B;
input Pty $ score number;
datalines;
A 100 1
B 90 2
C 82 4
D 51 6
E  59 7
G . .
;
run;

%CK_VALUE(DB1_DATA=work.DB_A,
DB2_DATA=work.DB_B,
BY=Pty,
ALL_VAR=score | number,
VAR_NUM=2)
結果:
NG:兩邊DB，VAR的觀測值不一樣。
NOT FOUND 兩邊MISSING。
OK:一至。

/*====================================================
DATA1     DB_A 資料集1    work.也要輸入
DATA2    DB_B 資料集2
BY  KEY
ALL_VAR 全部變數名稱 | 分隔
VAR_NUM 幾個變數
JIMMY CHEN
======================================================*/
%macro CK_VALUE(DB1_DATA=,
   DB2_DATA=,
   BY =,
   ALL_VAR = ,
   VAR_NUM=);
%let DATA1 = %scan(&DB1_DATA., 2,%str(.));
%let DATA2 = %scan(&DB2_DATA., 2,%str(.));

proc sort data = &DB1_DATA. out = DB1_REN nodupkey; by &BY.;run;
proc sort data = &DB2_DATA. out = DB2_REN nodupkey; by &BY.;run;

%do d = 1 %to 2;
 %do VAR_N = 1 %to &VAR_NUM.;
  %let VAR&VAR_N. = %scan(&ALL_VAR., &VAR_N.,%str(|));
/*  %put  &&VAR&VAR_N.. ;*/
/*  %put  &&VAR&VAR_N.._&&DATA&D..;*/
  data DB&d._REN;
   length &&VAR&VAR_N.._ck $30.;
   set DB&d._REN(rename =(&&VAR&VAR_N.. =&&VAR&VAR_N.._DB&d.));
  run;
 %end;
%end;                                  

%do VAR_N = 1 %to &VAR_NUM.;
 %let VAR&VAR_N. = %scan(&ALL_VAR., &VAR_N.,%str(|));
 data CK_&&VAR&VAR_N.. (keep =&BY. &&VAR&VAR_N.._ck);
  merge DB1_REN
     DB2_REN;
  by &BY.;
  if  &&VAR&VAR_N.._DB1 =  &&VAR&VAR_N.._DB2 then &&VAR&VAR_N.._ck = "OK";
   else  &&VAR&VAR_N.._ck = "NG";
  if missing(&&VAR&VAR_N.._DB1) then &&VAR&VAR_N.._ck = "NOT FOUND &DB1_DATA.";
  if missing(&&VAR&VAR_N.._DB2) then &&VAR&VAR_N.._ck = "NOT FOUND &DB2_DATA.";
  if missing(&&VAR&VAR_N.._DB1)  and  missing(&&VAR&VAR_N.._DB2)  then &&VAR&VAR_N.._ck = "NOT FOUND"; 

 run;
%end;


data CK_VALUE_ALL;
  merge DB1_REN
      DB2_REN;
  by &BY.;
run;

%do VAR_N = 1 %to &VAR_NUM.;
 %let VAR&VAR_N. = %scan(&ALL_VAR., &VAR_N.,%str(|));
 data  CK_VALUE_ALL;
  merge  CK_VALUE_ALL
     CK_&&VAR&VAR_N..;
    by &BY.;
 run;

 proc freq data = CK_VALUE_ALL;
  title "檢核觀測值: &&VAR&VAR_N..";
  table &&VAR&VAR_N.._ck;
 run;title;

proc print data = Ck_value_all;
 where &&VAR&VAR_N.._ck = "NG";
run;
%end;

%mend;

%CK_VALUE(DB1_DATA=amlland.Account_trade_extract,
  DB2_DATA=keepok.account_trade_extract,
  BY=TRANSACTION_REFERENCE_NUMBER,
  ALL_VAR=PRINCIPAL_AMOUNT_IN_BASE_CCY  | ACCOUNT_NUMBER, 
  VAR_NUM=2)

     
              
 

SAS Macro Meta-Analysis 統合分析

%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;