%macro fine_stratification (in_data= , exposure= , PS_provided= no , ps_var= , ps_class_var_list= , ps_cont_var_list= , PSS_method=exposure,
						    n_of_strata= , out_data= PS_FS, id_var= , effect_estimate= ,
							outcome= , survival_time= , out_excel= , ps_dist_plots= no, ps_dist_figure= );

*************************************************************************************************************

** step 0- Create the crude Table 1 **;

*************************************************************************************************************;

%table1 (in_for_table1= &in_data, treatment_var= &exposure, 
categorical_var_list= &ps_class_var_list, continuous_var_list= &ps_cont_var_list, out_table1= Crude_T1);

*************************************************************************************************************

** step 1- build and calculate the PS- if not already done **;

*************************************************************************************************************;
title 'PS model'; 

%let ps_provided1 = %upcase(&PS_provided); %if "&PS_provided1" = "NO" %then %do;

ods output Association=save_c; 

proc logistic data=&in_data desc ;
class &ps_class_var_list;
model &exposure= &ps_cont_var_list &ps_class_var_list; 
output out=sample_ps p=&ps_var;
run;

ods output clear; 
data save_c; set save_c; where label2='c'; keep label2 nvalue2; run;

%end;

%let ps_provided1 = %upcase(&PS_provided); %if "&PS_provided1" = "YES" %then %do;

data sample_ps; set &in_Data; run;

%end;

title; 

*************************************************************************************************************

** step 2- trim the non-overlapping regions of the PS **;

*************************************************************************************************************;

	*figure out the region of overlap;
	proc sql noprint;
	create table bounds_pscore as
	select max(&ps_var) as max_p, min(&ps_var) as min_p, &exposure
	from sample_ps
	group by &exposure;
	quit;

	*Save upper and lower bound of overlapping PS regions;

	proc sql noprint;
	select max(min_p) into: overlap_lower_bound from bounds_pscore; quit;
	proc sql noprint;
	select min(max_p) into: overlap_upper_bound from bounds_pscore; quit;

	*delete non-overlapping observations;

	data sample_ps_trimmed; set sample_ps;
	where &overlap_lower_bound<=&ps_var<=&overlap_upper_bound;
	run;

*************************************************************************************************************

** step 3- Conduct propensity score stratification 

*************************************************************************************************************;

%let PSS_method1 = %upcase(&PSS_method); %if "&PSS_method1" = "EXPOSURE" %then %do;

		** separate exposed only **;

		data only_exp; set sample_ps_trimmed; where &exposure=1; run;

		** Create strata based on the PS distribution of the exposed **;

		proc rank data=only_exp groups=&n_of_strata out=fs_exp;
		ranks strata_sas;
		var &ps_var;
		run;

		data fs_exp; set fs_exp; strata=1+strata_sas; * number of strata indicator, since proc rank assigns strata beginning from 0, we add 1 here **;
		drop strata_sas;
		run;

		*Get bounds of the PS for strata created;

				proc means data=fs_exp noprint;  
				 class strata; 
				 var &ps_var; 
				 output out=bounds_for_fs n=n 
				 min(&ps_var)=PS_min max(&ps_var)=PS_max; 
				run; 

				data bounds_for_fs; set bounds_for_fs;
				where strata ne .; *delete overall mean; 
				run;

		** Now, using the bounds of PS from the exposed distribution, assign strata in the unexposed patients;

		%let n=&n_of_strata; 

				%macro strata_in_unexp;

				** lets save the bounds for each separate strata into SAS datasets, which will be merged to the unexposed only 
				   cohort for the ease of strata creation in the next step **;

				%do strata= 1 %to &n;

					proc sql ; create table PSvalue&strata as
					select (ps_min) as ps_low&strata
					from bounds_for_fs
				    where strata=&strata; 
					quit;

				%end;

				data only_unexp; set sample_ps_trimmed; where &exposure=0; run;

				data only_unexp1; 
		              merge only_unexp psvalue1-psvalue&n;
				run;

				** the above merge creates a dataset with extra columns with each column containing 
				value of the lower bound of the stratum boundary calculated from the ps distribution of the exposed (eg 5 strata= 5 extra columns). 

				This information is used below to determine stratum membership of each unexposed subject. however, before
				moving on, we fill all the rows with the values of this variable by using missing=mean because the merge only fills one row since
				the datasets psvalue1 through n only has one row; 

				proc stdize out=only_unexp2 reponly missing=mean; run; 

				** create strata **;

				data unexp_strata; set only_unexp2;
				array abc(*) ps_low: ;
				do i=3 to dim(abc) while (strata=.);
				if 0<=&ps_var<ps_low2 then strata=1; ** the first strata between 0 and the minimum PS for the second stratum (stratum 1) **;
				if abc(i-1)<=&ps_var<abc(i) then strata=(i-1); ** rest of the strata based on the PS between minimum PS for two neighbouring 
				                                             strata. Note that since the last strata will be defined based on PS values
				                                             greater than the min PS in the penultimate strata, this step will code that
				                                             value as . In the next step, we give this strata the value equal to the 
															 highest strata**;
				end;
				run;

				data fs_unexp; set unexp_strata; if strata=. then strata=&n; 
				drop i ps_low:;
				run;

				/*
				proc means; var &ps_var; class strata; run;
				*/
				%mend;

				 %strata_in_unexp;

		** now concatenate the exposed and unexposed patients **;

		proc sort data=fs_exp; by &id_var; run;
		proc sort data=fs_unexp; by &id_var; run;

		data stratified;
		set fs_exp fs_unexp;
		by &id_var;
		run;

%end;

%let PSS_method1 = %upcase(&PSS_method); %if "&PSS_method1" = "COHORT" %then %do;

		proc rank data=sample_ps_trimmed groups=&n_of_strata out=fs;
		ranks strata_sas;
		var &ps_var;
		run;

		data stratified; set fs; 
		strata=1+strata_sas; 
		* number of strata indicator, since proc rank assigns strata beginning from 0, we add 1 here **;
		drop strata_sas;
		run;
			*Get bounds of the PS for strata created;

				proc means data=stratified noprint;  
				 class strata; 
				 var &ps_var; 
				 output out=bounds_for_fs n=n 
				 min(&ps_var)=PS_min max(&ps_var)=PS_max; 
				run; 

				data bounds_for_fs; set bounds_for_fs;
				where strata ne .;
				run;

%end;

*************************************************************************************************************

** step 4 Finally, assign weights to unexposed in each strata based on the distribution of the exposed

*************************************************************************************************************;

** assign weights to the unexposed according to the distribution of the exposed in each strata **;

		proc sort data=stratified out=matched_ps1; by strata; run;

		*total exposed; 

		data total_exposed; set matched_ps1; where &exposure=1; run;

		data total_exposed1; set total_exposed; by strata; if first.strata then total_Exp=0; total_Exp+1;
		if last.strata then output;
		keep strata total_exp;
		run;

		*total unexposed; 

		data total_unexposed; set matched_ps1; where &exposure=0; run;

		data total_unexposed1; set total_unexposed; by strata; 
		if first.strata then total_unExp=0; total_unExp+1;
		if last.strata then output;
		keep strata total_unexp;
		run;

		*merge the two files;

		data  strata_n; 
		merge  total_exposed1  total_unexposed1 ;
		by strata;
		if total_exp=. then total_exp=0; if total_unexp=. then total_unexp=0;
		run;

		*Drop the strata with either 0 exposed or 0 unexposed;

		data  strata_n; set  strata_n; if total_exp=0 or total_unexp=0 then delete; run;

		proc sql noprint; select sum(total_exp) into: sum_exp from  strata_n; quit; 
		
		proc sql noprint; select sum(total_unexp) into: sum_unexp from  strata_n; quit; 

		**Create weights for unexposed in each strata based on the distribution of exposed patients**;

		data  strata_n; set  strata_n; 
		weight=(total_exp/&sum_exp)/(total_unexp/&sum_unexp);
		weighted_n_unexp=total_unexp*((total_exp/&sum_exp)/(total_unexp/&sum_unexp));
		run;

		* QC step ;

		** Check the total- should be equal to the total number of unexposed we started with originally **;

		proc sql noprint; select sum (weighted_n_unexp) into: sum_unexp_weighted  from  strata_n; quit; 

		title 'Counts in each stratum- unweighted and weighted'; 

		proc print data=strata_n; run;

		** check numbers **;

		data numbers; 
		length n $50;
		n='After dropping uninformative strata';
		total_exposed= &sum_exp; total_unexposed=&sum_unexp; total_weighted_unexp=&sum_unexp_weighted;
		run;

		title 'Check if n of unexposed weighted = unweighted'; 

		proc print data=numbers noobs; run; 
		
		title; 

		*/

*** merge these weights with original cohort***;
		
proc sql;
create table with_weights as
select stratified.*, strata_n.weight
from stratified
left join strata_n
on stratified.strata=strata_n.strata
where weight ne .; *delete the uninformative strata;
quit;

		*** these weights only apply to unexposed, every exposed person get the weight of 1 ***;

		data &out_data; set with_weights;
		psweight=1; if &exposure=0 then psweight=weight; 
		drop weight;
		run;


*************************************************************************************************************

** step 5- Create the propensity score weighted Table 1 & calculate post-weighting c statistic
		   as a measure of improvement in balance between the two groups in the covariate vector**;

*************************************************************************************************************;

%table1 (in_for_table1= &out_data, treatment_var= &exposure, 
categorical_var_list= &ps_class_var_list, continuous_var_list= &ps_cont_var_list, weight=psweight, out_table1= Weighted_T1);

*title 'Post-weighting PS model'; 

ods graphics off;           
ods exclude all;             
ods output Association=save_c1; 

proc logistic data=&out_data desc ;
class &ps_class_var_list;
model &exposure= &ps_cont_var_list &ps_class_var_list; 
weight psweight;
run;

ods output clear; 

ods exclude none; 

data save_c1; set save_c1; where label2='c'; keep label2 nvalue2; run;

%if "&PS_provided1" = "NO" %then %do;
data balance; length measure $30; set save_c (in=r) save_c1 (in=s); 
if r=1 then measure='PS model c-statistic'; else measure='Post weighting c-statistic';
run;
%end; 

%if "&PS_provided1" = "YES" %then %do;
data balance; length measure $30; set save_c1; 
measure='Post weighting c-statistic';
run;
%end; 


*************************************************************************************************************

** step 6- Estimate the effect of treatment using weighted log reg model or weighted PH model

*************************************************************************************************************;

** generalized linear model with identity link for risk differences; 

%let effect_estimate1 = %upcase(&effect_estimate); %if "&effect_estimate1" = "RD" %then %do;

* crude and weighted event counts; 

proc sql; create table crude_counts as 
select &exposure, sum (&outcome) as n_crude_events, count (&id_var) as n_crude_total, (calculated n_crude_events/calculated n_crude_total) as crude_risk
from &in_data
group by &exposure;
quit;

proc sql; create table weighted_counts as 
select &exposure, round(sum (&outcome*psweight)) as n_weighted_events, sum(psweight) as n_weighted_total, (calculated n_weighted_events/calculated n_weighted_total) as weighted_risk
from &out_data
group by &exposure;
quit;

data event_counts; merge crude_counts weighted_counts; by &exposure; run;

*crude analysis; 

title 'Unadjusted model'; 

ods output ParameterEstimates=crude_beta;
proc genmod data=&in_data desc;
class &exposure/param=ref ref=first;
model &outcome= &exposure/dist=bin link=identity;
run;
ods output clear; 

data result1; set crude_beta; where parameter not in ('Intercept','Scale'); 
RD=(Estimate); LCL=(Estimate-1.96*StdErr); UCL=(Estimate+1.96*StdErr);
keep parameter estimate stderr RD LCL UCL;
run;


*ps-adjusted analysis; 

title 'PS-weighted model'; 

ods output ParameterEstimates=weighted_beta ;
proc genmod data=&out_data desc;
class &exposure/param=ref ref=first;
model &outcome= &exposure/dist=bin link=identity;
weight psweight;
run;
ods output clear; 

data result2; set weighted_beta; where parameter not in ('Intercept','Scale'); 
RD=(Estimate); LCL=(Estimate-1.96*StdErr); UCL=(Estimate+1.96*StdErr);
keep parameter estimate stderr RD LCL UCL;
run;

data results; length method $20;  set result1 (in=r) result2; 
if r=1 then Method= 'Crude'; else Method='PS-weighted'; 
run;

title; 

%end;

** generalized linear model with log link for risk ratios; 

%let effect_estimate1 = %upcase(&effect_estimate); %if "&effect_estimate1" = "RR" %then %do;

* crude and weighted event counts; 

proc sql; create table crude_counts as 
select &exposure, sum (&outcome) as n_crude_events, count (&id_var) as n_crude_total, (calculated n_crude_events/calculated n_crude_total) as crude_risk
from &in_data
group by &exposure;
quit;

proc sql; create table weighted_counts as 
select &exposure, round(sum (&outcome*psweight)) as n_weighted_events, sum(psweight) as n_weighted_total, (calculated n_weighted_events/calculated n_weighted_total) as weighted_risk
from &out_data
group by &exposure;
quit;

data event_counts; merge crude_counts weighted_counts; by &exposure; run;

** crude analysis; 

title 'Unadjusted model'; 

ods output ParameterEstimates=crude_beta;
proc genmod data=&in_data desc;
class &exposure/param=ref ref=first;
model &outcome= &exposure/dist=bin link=log;
run;
ods output clear; 

data result1; set crude_beta; where parameter not in ('Intercept','Scale'); 
RR=exp(Estimate); LCL=exp(Estimate-1.96*StdErr); UCL=exp(Estimate+1.96*StdErr);
keep parameter estimate stderr RR LCL UCL;
run;

** ps adjusted; 

title 'PS-weighted model'; 

ods output ParameterEstimates=weighted_beta ;
proc genmod data=&out_data desc;
class &exposure/param=ref ref=first;
model &outcome= &exposure/dist=bin link=log;
weight psweight;
run;
ods output clear; 

data result2; set weighted_beta; where parameter not in ('Intercept','Scale'); 
RR=exp(Estimate); LCL=exp(Estimate-1.96*StdErr); UCL=exp(Estimate+1.96*StdErr);
keep parameter estimate stderr RR LCL UCL;
run;

data results; length method $20;  set result1 (in=r) result2; 
if r=1 then Method= 'Crude'; else Method='PS-weighted'; 
run;

title;

%end;

** generalized linear model with logit link for odds ratios; 

%let effect_estimate1 = %upcase(&effect_estimate); %if "&effect_estimate1" = "OR" %then %do;

* crude and weighted event counts; 

proc sql; create table crude_counts as 
select &exposure, sum (&outcome) as n_crude_events, count (&id_var) as n_crude_total, (calculated n_crude_events/calculated n_crude_total) as crude_risk
from &in_data
group by &exposure;
quit;

proc sql; create table weighted_counts as 
select &exposure, round(sum (&outcome*psweight)) as n_weighted_events, sum(psweight) as n_weighted_total, (calculated n_weighted_events/calculated n_weighted_total) as weighted_risk
from &out_data
group by &exposure;
quit;

data event_counts; merge crude_counts weighted_counts; by &exposure; run;

**Crude analysis;

title 'Unadjusted model'; 

ods output ParameterEstimates=crude_beta;
proc genmod data=&in_data desc;
class &exposure/param=ref ref=first;
model &outcome= &exposure/dist=bin link=logit;
run;
ods output clear; 

data result1; set crude_beta; where parameter not in ('Intercept','Scale'); 
OR=exp(Estimate); LCL=exp(Estimate-1.96*StdErr); UCL=exp(Estimate+1.96*StdErr);
keep parameter estimate stderr OR LCL UCL;
run;

** ps adjusted; 

title 'PS-weighted model'; 

ods output ParameterEstimates=weighted_beta ;
proc genmod data=&out_data desc;
class &exposure/param=ref ref=first;
model &outcome= &exposure/dist=bin link=logit;
weight psweight;
run;
ods output clear; 

data result2; set weighted_beta; where parameter not in ('Intercept','Scale'); 
OR=exp(Estimate); LCL=exp(Estimate-1.96*StdErr); UCL=exp(Estimate+1.96*StdErr);
keep parameter estimate stderr OR LCL UCL;
run;

data results; length method $20;  set result1 (in=r) result2; 
if r=1 then Method= 'Crude'; else Method='PS-weighted'; 
run;

title; 

%end;

** cox proportional hazard regression model with time to event data for hazard ratios; 

%let effect_estimate1 = %upcase(&effect_estimate); %if "&effect_estimate1" = "HR" %then %do;

* crude and weighted event counts; 

proc sql; create table crude_counts as 
select &exposure, sum (&outcome) as n_crude_events, sum(&survival_time) as cumulative_ptime_crude, (calculated n_crude_events/calculated cumulative_ptime_crude) as crude_incidence_rate
from &in_data
group by &exposure;
quit;

proc sql; create table weighted_counts as 
select &exposure, round(sum (&outcome*psweight)) as n_weighted_events, round(sum(&survival_time*psweight)) as cumulative_ptime_weighted, (calculated n_weighted_events/calculated cumulative_ptime_weighted) 
as weighted_incidence_rate
from &out_data
group by &exposure;
quit;

data event_counts; merge crude_counts weighted_counts; by &exposure; run;

** crude analysis**;

title 'Unadjusted model'; 

ods output ParameterEstimates=crude_beta; 
  proc phreg data=&in_data;
  class &exposure/param=ref ref=first;
  model &survival_time*&outcome(0) = &exposure/rl;
  title 'Crude';
  run;
ods output clear; 

** adjusted analysis **;

title 'PS-weighted model'; 

ods output ParameterEstimates=weighted_beta; 

  proc phreg data=&out_data ;
  class &exposure/param=ref ref=first;
     model &survival_time*&outcome(0) = &exposure/rl;
     weight psweight;
     run;
ods output clear; 

data results; length method $20;  set crude_beta (in=r) weighted_beta; 
if r=1 then Method= 'Crude'; else Method='PS-weighted'; 
drop chisq probchisq classval0 df;
run;

title; 

%end;

*************************************************************************************************************

** step 7- Save all the outputs in an excel file

*************************************************************************************************************;

ods tagsets.excelxp file="&out_excel..xml" style=journal
options(sheet_interval='yes');

title 'Crude Table 1';

proc print data=Crude_T1 noobs;
var Table Category Treated_freq Treated_percent Ref_freq Ref_percent Treatment_column Ref_column crude_diff std_diff;
run;

title 'Weighted Table 1';

proc print data=Weighted_T1 noobs;
var Table Category Treated_freq Treated_percent Ref_freq Ref_percent Treatment_column Ref_column crude_diff std_diff;
run;

title 'Balance';

proc print data=balance noobs;
run;

title 'Event counts';

proc print data=event_counts noobs;
run;

title 'Effect estimates';

proc print data=results noobs;
run;

ods tagsets.excelxp close;

title; 

*************************************************************************************************************

** step 8- Draw the PS distribution plots

*************************************************************************************************************;

%let ps_dist_plots1 = %upcase(&ps_dist_plots); %if "&ps_dist_plots1" = "YES" %then %do;

%macro ps_overlap_plots (in_for_plot= , exp_var= &exposure, prop_score= &ps_var , bin_width= 0.1, max_ps=1, ps_boundaries_dataset= dummy, 
                         figure_name= , figure_title= );

** First step is outputting exposed and uexposed groups separately **;

data cohort_ex cohort_unex; 
set &in_for_plot; 
if &exp_var=1 then output  cohort_ex;
if &exp_var=0 then output  cohort_unex;
keep &exp_var &prop_score psweight;
run;

** second is create separate datasets with % of patients in desired histogram width **;
 
%let n_bins=(1/&bin_width); ** internal macro variable for number of bins based on the width provided;

data s1; set cohort_ex; 
his_bin=1;
	do i=1 to (&n_bins-1);
	if (i*&bin_width)< &prop_score <=((i*&bin_width)+&bin_width) then his_bin=(i+1);
	end;
run; 

proc freq noprint; tables his_bin/out=freq_exposed;
run;

data s0; set cohort_unex; 
his_bin=1;
	do i=1 to (&n_bins-1);
	if (i*&bin_width)< &prop_score <=((i*&bin_width)+&bin_width) then his_bin=(i+1);
	end;
run; 

proc freq noprint; tables his_bin/out=freq_unexposed; weight psweight;
run;

** merge the two output tables and make necessary changes **;

data bins;
merge freq_exposed (rename= (percent= p_exp)) freq_unexposed (rename= (percent= p_unexp));
by his_bin;
drop count;
run;

data zero; his_bin=0; * a zero category, needed to create smooth curves ;
p_exp=0; p_unexp=0;
run;

data bins1; set zero bins ;
run;

data bins1; set bins1;
PS=his_bin*&bin_width; ** create the variable indicating the propensity score since the variable his_bin
					      is a categorical variable showing the number of bin, which is not meaningful **;
label p_exp= 'Treated';
label p_unexp='Referent';
run;

proc stdize reponly missing=0;
run;

** dummy dataset to create fake values for reference lines to make sure the lines are plotted even when fewer values of PS strata than 100 are provided **;

data dummy; do i=1 to 100; dvar=i+2; output; end; run; 
data dummy1; set dummy; ps_min=99; run;
proc sql noprint; select dvar into: r1-:r100 from dummy; quit; 

** Create macro variables from the PS bounds dataset to plot reference lines on the plot **; 

proc sql noprint; select ps_min into: r1-:r100 from &ps_boundaries_dataset; quit; 

        goptions reset=all rotate=landscape ftitle=swiss ftext=swiss
        gunit=pct htitle=3.5 htext =2.5 ;

    symbol1 v=none i=spline w=3 color=gray00;   
   	symbol2 v=none i=spline line=3 color=gray00 w=3;  
   	axis1 label=('Propensity scores') w=3 ;
	axis2 label = ('% patients' ) w = 3 ;
	legend1 label=none position=(top right inside);
	proc gplot data=bins1; where ps<=1;
		plot p_exp*PS=1 p_unexp*PS=2/ legend=legend1 overlay frame vaxis=axis2 haxis=axis1 
    href=&r1 &r2 &r3 &r4 &r5 &r6 &r7 &r8 &r9 &r10 
         &r11 &r12 &r13 &r14 &r15 &r16 &r17 &r18 &r19 &r20 
		 &r21 &r22 &r23 &r24 &r25 &r26 &r27 &r28 &r29 &r30 
		 &r31 &r32 &r33 &r34 &r35 &r36 &r37 &r38 &r39 &r40 
	     &r41 &r42 &r43 &r44 &r45 &r46 &r47 &r48 &r49 &r50 
		 &r51 &r52 &r53 &r54 &r55 &r56 &r57 &r58 &r59 &r60 
		 &r61 &r62 &r63 &r64 &r65 &r66 &r67 &r68 &r69 &r70 
		 &r71 &r72 &r73 &r74 &r75 &r76 &r77 &r78 &r79 &r80 
		 &r81 &r82 &r83 &r84 &r85 &r86 &r87 &r88 &r89 &r90 
		 &r91 &r92 &r93 &r94 &r95 &r96 &r97 &r98 &r99 &r100
    name=&figure_name;
    title &figure_title;
	RUN;
	quit;

%mend;

data sample_ps1; set sample_ps; psweight=1; run; *since this is unweighted plot, assign weight of 1 to everyone; 

%ps_overlap_plots (in_for_plot= sample_ps1, exp_var= &exposure, prop_score= &ps_var , bin_width= 0.1, max_ps=1, ps_boundaries_dataset= dummy1, 
                         figure_name= 'Figure1', figure_title= 'Crude PS distribution'); 

* the following macro will created a weighted PS distribution by accounting for PS weights while creating histogram bins; 

%ps_overlap_plots (in_for_plot= &out_data, exp_var= &exposure, prop_score= &ps_var , bin_width= 0.1, max_ps=1, ps_boundaries_dataset= Bounds_for_fs, 
                         figure_name= 'Figure2', figure_title= 'Weighted PS distribution'); 

** combine the two plots on the same page; 


options orientation=landscape nodate nonumber;

ods listing close;
ods pdf file="&ps_dist_figure..pdf" notoc ;  

proc greplay nofs tc=work.templt igout=work.gseg; **     1) fs-full screen
                                                    2) TC=template-catalog
                                                    specifies the template catalog to use with the GREPLAY procedure,
                                                    and identifies the template catalog where the template entry is
                                                    stored to replay your graphics,
                                                    3) IGOUT=<libref.>input-catalog
                                                    specifies the input catalog that stores the
                                                    graphics output that you want to use with the GREPLAY procedure.
                                                    If you omit the libref, SAS/GRAPH looks for the catalog in the
                                                    temporary WORK library.;

          tdef graf2by1

          1 / ulx = 0    uly = 100
              urx = 100    ury = 100
              llx = 0   lly = 50
              lrx = 100   lry = 50

          2 / ulx = 0    uly = 47
              urx = 100    ury = 47
              llx = 0   lly = 0
              lrx = 100   lry = 0;

    template graf2by1; ** tell sas to use the template you just created;

    treplay 1:FIGURE1 2:FIGURE2;

    *Replays graphics entries into template panels. TREPLAY copies one or more entries from the graphics input
    catalog into a template-entry in the graphics output catalog, using positioning information provided by the
    template.;

run;


ods pdf close; 
ods listing;

proc greplay igout = work.gseg nofs;
delete _all_;
run;

quit;

%end;

title; 

%mend;