Plotting PROC PHREG after Multiple Imputation

[fikel]

Hello, I had previously asked about plotting Proc PHREG for restricted splines analysis and was able to use this code:

proc phreg data = have;
title "Colorectal-cancer incidence attempt";
	
effect resid_total_polys = spline(resid_total_poly / basis=tpf(noint) NATURALCUBIC details knotmethod=Percentilelist(5 35 65 95));

class menopause_with_males (ref="0") race (ref="1") hhincome (ref="1")
		enrollment_source (ref="C") fh_colorectalcancer (ref="0")
		BMI_category_expanded (ref="0") alcohol_category (ref="0") smokestatus_packyear (ref="0") hei10_category (ref="0");
 
model eof_age_CRC*Inc_CRC(0) = resid_total_polys 
		menopause_with_males race hhincome enrollment_source 
		BMI_category_expanded totalactivitymethr smokestatus_packyear Comorbidity_Index 
		ffq_kcal hei10_category alcohol_category fh_colorectalcancer
		/ entry=enrollment_agemonths rl=wald;
		
store work.coxr;
run;



%macro est(ref=212, start=0, end=3000, by=1);

%Do i = 1 %To %eval(%SysFunc( Ceil( %SysEvalF( ( &End - &Start ) / &By ) ) ) +1) ;

   %Let value=%SysEvalF( ( &Start - &By ) + ( &By * &I ) ) ;

    estimate "&value." resid_total_polys [-1, &ref] [1, &value] / exp cl;

    %end;

%mend est;

ods exclude all;

ods dataset Estimates=Estimates;

proc plm restore=work.coxr;

    %est(ref=212, start=0, end=3000, by=1);

run;

ods exclude none;

data estimates;

    set estimates;

    resid_total_poly=label*1;

run;

proc sgplot data=estimates NOAUTOLEGEND Noborder;

    Series y=ExpEstimate x=resid_total_poly / LINEATTRS=(color=black thickness=3px);

    Series y=LowerExp x=resid_total_poly / LINEATTRS=(pattern=ShortDash color=Black THICKNESS=1);

    Series y=UpperExp x=resid_total_poly / LINEATTRS=(pattern=ShortDash color=Black THICKNESS=1);
    
    band x=resid_total_poly upper=UpperExp lower=LowerExp / fillattrs=(color=graydd) transparency=.50;

    REFLINE 1 / axis=y;

    REFLINE 3000 / Axis=X LINEATTRS=(pattern=ThinDot color=Black THICKNESS=1);;

    yaxis Values=(0.5 0.8 1 1.2 1.5) Label="Hazard ratio" Type=LOG LABELATTRS=(weight=BOLD);

    xaxis min=0 VALUES=(0 to 3000 by 100) LABELATTRS=(weight=BOLD);

run;

To produce this graph:

Now, however, I am using multiple imputation to handle missing values, but I am not sure how to adapt the code to be able to handle this. This is the code I have:

proc phreg data = mult_imp.mi_mvn namelen=50;
title "All-cancer incidence attempt";

effect resid_total_polys = spline(resid_total_poly / basis=tpf(noint) NATURALCUBIC details knotmethod=Percentilelist(5 35 65 95));

model eof_age_CRC*Inc_CRC(0) = resid_total_polys
		menopause_males_dum2 menopause_males_dum3 racedum2 racedum3 incomedum2 incomedum3 incomedum4 
		enrolldum2 BMI_dum2 BMI_dum3 totalactivitymethr packyears_dum2 packyears_dum3 packyears_dum4 packyears_dum5 Comorbidity_Index 
		ffq_kcal hei10_category alcohol_dum2 alcohol_dum3 fh_colorectalcancer 
		/ entry=enrollment_agemonths rl=wald;
by _imputation_;

store work.coxr;

run;
quit;

%macro est(ref=212, start=0, end=3000, by=1);
%Do i = 1 %To %eval(%SysFunc( Ceil( %SysEvalF( ( &End - &Start ) / &By ) ) ) +1) ;
   %Let value=%SysEvalF( ( &Start - &By ) + ( &By * &I ) ) ;
    estimate "&value." resid_total_polys [-1, &ref] [1, &value] / exp cl;
    %end;
%mend est;

ods exclude all;
ods dataset Estimates=Estimates;

proc plm restore=work.coxr;
    %est(ref=212, start=0, end=3000, by=1);
run;

ods exclude none;
data estimates;
    set estimates;
    resid_total_poly=label*1;
run;

proc sgplot data=estimates NOAUTOLEGEND Noborder;
    Series y=ExpEstimate x=resid_total_poly / LINEATTRS=(color=black thickness=3px);
    Series y=LowerExp x=resid_total_poly / LINEATTRS=(pattern=ShortDash color=Black THICKNESS=1);
    Series y=UpperExp x=resid_total_poly / LINEATTRS=(pattern=ShortDash color=Black THICKNESS=1); 
    band x=resid_total_poly upper=UpperExp lower=LowerExp / fillattrs=(color=graydd) transparency=.50;
    REFLINE 1 / axis=y;
    REFLINE 3000 / Axis=X LINEATTRS=(pattern=ThinDot color=Black THICKNESS=1);;
    yaxis Values=(0.5 0.8 1 1.2 1.5) Label="Hazard ratio" Type=LOG LABELATTRS=(weight=BOLD);
    xaxis min=0 VALUES=(0 to 3000 by 100) LABELATTRS=(weight=BOLD);

run;

Where the set "mult_imp.mi_mvn" is the dataset obtained after running Proc MI with 5 imputations. I am stuck trying to figure out how to get a single graph for this model that combines the 5 imputed sets. However, this produces this plot: 

I was thinking that I somehow need to use proc mianalyze to combine the estimates for the models from each imputation but was not sure how to accomplish this.