topic Re: How do I plot the spline effect in PHREG? in SAS Health and Life Sciences

How do I plot the spline effect in PHREG?

wenxinxu — Wed, 13 Feb 2013 16:13:06 GMT

Hi,

I am new to SAS and am having trouble graphing the SPLINE constructed effect in PHREG.

My question is, what is the best way to draw a spline using the output from the SPLINE effect in PHREG? This is a survival analysis so it seems to me that using TRANSREG would not work as well.

/* Here is the PHREG statement I am using: */

PROC PHREG DATA=WORK.TMP0TempTableInput
;
CLASS smoking(ref = 'N');
CLASS AgeCategory(ref = '1');
EFFECT spl = spline(BMI / degree=3 knotmethod=rangefractions(.05 .25 .75 .95) Details);
MODEL MonthsToDeathEnd * CensoredDeath (1) = spl Female Smoking AgeCategory HxCAD HxCHF HxCVA HxPVD HxDiabetes HxCOPD Charlson weightedMPR weightedtownsend /
TIES=BRESLOW
RISKLIMITS ALPHA=0.05
SELECTION=NONE
;

RUN;TITLE;

/* Here is an example of the type of output I am trying to generate: */

Thanks!

Re: How do I plot the spline effect in PHREG?

DanielleSchoenaker — Thu, 30 May 2013 01:02:19 GMT

Hi,

Just wondering whether you know by now how to generate the plot you showed? I am aiming for a similar plot after a similar PHREG statement.

Many thanks in advance!

Re: How do I plot the spline effect in PHREG?

JacobSimonsen — Sun, 04 May 2014 09:43:41 GMT

I sometime also deals with this problem. I have not found a way to plot the spline-curve created by the effect-estimates. I therefore only use the spline-feature in the effect-statement when I want to adjust for some covariate which own effect is not of interest.

If I want to plot the spline, then I do all the coding wihout use of the effectstatement. That is (1) create the coefficients which will be used as regressors in phreg (or some other regression procedure), (2) from the estimates a plot dataset can be created and at last (3) plot the curves with proc gplot. It can be it can be done more easily, so I encourage others to come with better solutions. With some small modifications, it is also possible to use the code below to plot a time-dependent effect in a cox-regression.

There are some matrix-calculations, and it is thefore neccesary to have run the matrix-package first. The matrixalgebra file which are included is the one that is attached this post:

Here is an example:

*k is is the number of knots;

%let k=9;

%include "" ; /*insert path to the matrix package */

option cmplib=work.func; /*this should be the same location as the one used in the matrix-package*/

*for this example I also simulate data;

data simulation;

array spl{%eval(&k.-1)};

array knots{&k.} _temporary_ (-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4) ; /*this is my knot-points*/

do j=1 to 10000;

* a is my covariate, which is here some number between zero and one.;

a=ranuni(-1)-0.5;

*t is my survival time;

t=-log(ranuni(-1))/(exp(-(a)**2));

*then create the spline-coefficients;

ref=0; *0 is here my reference value;

spl[1]=a-ref;

do i=1 to &k.-2;

spl[i+1]=max(a-knots,0)**3

-((knots[&k]-knots)/(knots[&k]-knots[&k-1]))*max(a-knots[&k.-1],0)**3

+((knots[&k-1]-knots)/(knots[&k]-knots[&k-1]))*max(a-knots[&k.],0)**3

-(max(ref-knots,0)**3

-((knots[&k]-knots)/(knots[&k]-knots[&k-1]))*max(ref-knots[&k.-1],0)**3

+((knots[&k-1]-knots)/(knots[&k]-knots[&k-1]))*max(ref-knots[&k.],0)**3);

end;

output;

end;

run;

*because I also wants pointwise confidence-limits I put the covariance-matrix into a dataset "outest";

proc phreg data=simulation covout outest=outest;

model t=spl1-spl%eval(&k.-1);

run;

data plot;

array knots{&k.} _temporary_ (-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4) ;

array spl{%eval(&k.-1)};

array variance{%eval(&k.-1),%eval(&k.-1)} _temporary_;

array beta{%eval(&k.-1),1} _temporary_;

array design{1,%eval(&k.-1)} _temporary_;

array designt{%eval(&k.-1),1} _temporary_;

array logest{1,1} _temporary_;

array var{1,1} _temporary_;

*read in the estimates;

do until (eof);

set outest end=eof;

do i=1 to &k-1.;

if _type_='PARMS' then beta[i,1]=spl;

else variance[i,input(substr(_name_,4,length(_name_)-3),best.)]=spl;

end;

call show(variance);

*create the curve and confidencelimits;

ref=0;

do a=-0.5 to 0.5 by 0.001;

design[1,1]=a-ref;

do i=1 to &k.-2;

design[1,i+1]=max(a-knots,0)**3

-((knots[&k]-knots)/(knots[&k]-knots[&k-1]))*max(a-knots[&k.-1],0)**3

+((knots[&k-1]-knots)/(knots[&k]-knots[&k-1]))*max(a-knots[&k.],0)**3

-(max(ref-knots,0)**3

-((knots[&k]-knots)/(knots[&k]-knots[&k-1]))*max(ref-knots[&k.-1],0)**3

+((knots[&k-1]-knots)/(knots[&k]-knots[&k-1]))*max(ref-knots[&k.],0)**3);

end;

call trans(design,designt);

call multiplicer(design,beta,logest);

call multiplicer3(design,variance,designt,var);

rr=exp(logest[1,1]);

lower=exp(logest[1,1]-1.96*sqrt(var[1,1]));

upper=exp(logest[1,1]+1.96*sqrt(var[1,1]));

output;

end;

keep a rr lower upper;

run;

symbol1 i=join v=none l=1 c=black;

symbol2 i=join v=none l=2 c=black;

symbol3 i=join v=none l=2 c=black;

proc gplot data=plot;

plot (rr lower upper)*a/overlay;

run;

Re: How do I plot the spline effect in PHREG?

TessaRH — Wed, 01 Apr 2015 12:46:10 GMT

Hi,

I want to do the same approximately.

Current code:

ods graphics on;

proc phreg data=splinedata ;

effect lala=spline(log_lab/knotmethod=list(1 3 5));

model time*event(0)=lala;

run;

Anyone knows how to plot this? I want to plot on the y-axis HR, and a continuous variable on the x-axis.

I'm using SAS 9.3.

Anyone familiar with the %rcs_reg macro by Desquibet (Dose-response analyses using restricted cubic spline functions in public health research - Desquilbet - 2010 - Statistics in Medicine - Wiley Online Library)?

Thanks!

Tessa

Re: How do I plot the spline effect in PHREG?

LinusS_ — Wed, 03 Jun 2015 15:13:45 GMT

Hi,

I think something like this should do it:

proc phreg data=temp;

effect bmiS = spline(bmi / basis=tpf(noint) NATURALCUBIC details knotmethod=list(18.0175 20.5289 22.4059 27.4693) );

model time*event(0)=bmiS / rl=wald ties=EFRON ;

store sasuser.coxr;

run;

%macro est(ref=20, start=15, end=35, by=0.1);

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

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

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

%end;

%mend est;

ods html select none;

ods rtf select none;

ods dataset Estimates=Estimates;

proc plm restore=sasuser.coxr;

%est(ref=20, start=16, end=35, by=0.01);

run;

ods html select all;

ods rtf select all;

data estimates;

set estimates;

BMI=label*1;

run;

proc sgplot data=estimates NOAUTOLEGEND Noborder;

Series y=ExpEstimate x=BMI / ;

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

Series y=UpperExp x=BMI / LINEATTRS=(pattern=ShortDash color=Black THICKNESS=1);

REFLINE 1 / axis=y;

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

yaxis Values=(0.8 1 1.5 2 3 4 6 8 10) Label="Hazard ratio" Type=LOG LABELATTRS=(weight=BOLD);

xaxis min=16 VALUES=(16 to 35 by 1) LABELATTRS=(weight=BOLD);

run;

Re: How do I plot the spline effect in PHREG?

bb_heihei — Wed, 26 Aug 2015 03:53:14 GMT

Hi LinusS.,

Thank you! This works well for PHREG. I am also looking for a way to draw a spline based on the survey data. I think the above code would not work well for the PROC SURVEYPHREG, because there is not the effect statement in SURVEYPHREG. Is there a best way to draw a spline after SURVEYPHREG?

Thanks a lot!

Re: How do I plot the spline effect in PHREG?

deepnorth3 — Wed, 08 Nov 2017 01:07:47 GMT

Hi, can you please tell me what is matrix package? Is it a R package?

Re: How do I plot the spline effect in PHREG?

SriniR — Thu, 20 Feb 2020 21:28:00 GMT

@JacobSimonsen wrote:
I sometime also deals with this problem. I have not found a way to plot the spline-curve created by the effect-estimates. I therefore only use the spline-feature in the effect-statement when I want to adjust for some covariate which own effect is not of interest.

If I want to plot the spline, then I do all the coding wihout use of the effectstatement. That is (1) create the coefficients which will be used as regressors in phreg (or some other regression procedure), (2) from the estimates a plot dataset can be created and at last (3) plot the curves with proc gplot. It can be it can be done more easily, so I encourage others to come with better solutions. With some small modifications, it is also possible to use the code below to plot a time-dependent effect in a cox-regression.

There are some matrix-calculations, and it is thefore neccesary to have run the matrix-package first. The matrixalgebra file which are included is the one that is attached this post:

Here is an example:

*k is is the number of knots;
%let k=9;

%include "" ; /*insert path to the matrix package */
option cmplib=work.func; /*this should be the same location as the one used in the matrix-package*/

*for this example I also simulate data;
data simulation;
array spl{%eval(&k.-1)};
array knots{&k.} _temporary_ (-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4) ; /*this is my knot-points*/
do j=1 to 10000;
   * a is my covariate, which is here some number between zero and one.;
    a=ranuni(-1)-0.5;
*t is my survival time;
t=-log(ranuni(-1))/(exp(-(a)**2));

*then create the spline-coefficients;
ref=0; *0 is here my reference value;
spl[1]=a-ref;
do i=1 to &k.-2;
spl[i+1]=max(a-knots,0)**3
     -((knots[&k]-knots)/(knots[&k]-knots[&k-1]))*max(a-knots[&k.-1],0)**3
    +((knots[&k-1]-knots)/(knots[&k]-knots[&k-1]))*max(a-knots[&k.],0)**3
-(max(ref-knots,0)**3
   -((knots[&k]-knots)/(knots[&k]-knots[&k-1]))*max(ref-knots[&k.-1],0)**3
   +((knots[&k-1]-knots)/(knots[&k]-knots[&k-1]))*max(ref-knots[&k.],0)**3);
end;
    output;
end;
run;

*because I also wants pointwise confidence-limits I put the covariance-matrix into a dataset "outest";
proc phreg data=simulation covout outest=outest;
model t=spl1-spl%eval(&k.-1);
run;

data plot;
array knots{&k.} _temporary_ (-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4) ;
array spl{%eval(&k.-1)};
array variance{%eval(&k.-1),%eval(&k.-1)} _temporary_;
array beta{%eval(&k.-1),1} _temporary_;
array design{1,%eval(&k.-1)} _temporary_;
array designt{%eval(&k.-1),1} _temporary_;
array logest{1,1} _temporary_;
array var{1,1} _temporary_;

*read in the estimates;
do until (eof);
    set outest end=eof;
    do i=1 to &k-1.;
      if _type_='PARMS' then beta[i,1]=spl;
      else variance[i,input(substr(_name_,4,length(_name_)-3),best.)]=spl;
    end;
end;
call show(variance);

*create the curve and confidencelimits;
ref=0;
do a=-0.5 to 0.5 by 0.001;
design[1,1]=a-ref;
do i=1 to &k.-2;
design[1,i+1]=max(a-knots,0)**3
-((knots[&k]-knots)/(knots[&k]-knots[&k-1]))*max(a-knots[&k.-1],0)**3
+((knots[&k-1]-knots)/(knots[&k]-knots[&k-1]))*max(a-knots[&k.],0)**3
-(max(ref-knots,0)**3
-((knots[&k]-knots)/(knots[&k]-knots[&k-1]))*max(ref-knots[&k.-1],0)**3
+((knots[&k-1]-knots)/(knots[&k]-knots[&k-1]))*max(ref-knots[&k.],0)**3);
end;
    call trans(design,designt);
    call multiplicer(design,beta,logest);
    call multiplicer3(design,variance,designt,var);
    rr=exp(logest[1,1]);
    lower=exp(logest[1,1]-1.96*sqrt(var[1,1]));
    upper=exp(logest[1,1]+1.96*sqrt(var[1,1]));
    output;
end;
keep a rr lower upper;
run;
symbol1 i=join v=none l=1 c=black;
symbol2 i=join v=none l=2 c=black;
symbol3 i=join v=none l=2 c=black;
proc gplot data=plot;
plot (rr lower upper)*a/overlay;
run;

Hello, It has been a while since post appeared. I was recently doing to search for an approach to plot HR against a continuous variable and I came across this thread. You are referring to a matrix package attached to your post ("There are some matrix-calculations, and it is therefore necessary to have run the matrix-package first. The matrixalgebra file which are included is the one that is attached this post"). There does not seem to be any attachment. Could you re-post your message with the attachment? Or, if the matrix package is embedded in the post, could you point to that? Would greatly appreciate your help.

Regards,

Srini Rajagopalan

Re: How do I plot the spline effect in PHREG?

sseb — Wed, 23 Dec 2020 16:37:18 GMT

Thanks @LinusS_ for the code! I'm wondering if you know how to modify it to get the estimates by sex? I tried the following but it didn't work, it looks especially wonky for females.

%macro est(ref=20, start=15, end=35, by=0.1);
   %Do i = 1 %To %eval(%SysFunc( Ceil( %SysEvalF( ( &End - &Start ) / &By ) ) ) +1) ;
      %Let value=%SysEvalF( ( &Start - &By ) + ( &By * &I ) ) ;
         estimate "&value. male" sex 0 bmis [-1, &ref] [1, &value] / exp cl;
         estimate "&value. female" sex 1 bmis [-1, &ref] [1, &value] / exp cl;
   %end;
%mend est;

Re: How do I plot the spline effect in PHREG?

LinusS_ — Wed, 23 Dec 2020 20:50:11 GMT

@sseb What you're trying to estimate doesn't really make sense unless you have an interaction between sex and BMI (without an interaction the HR is the same for both males and females). If you include an interaction I think it would be something like:

estimate "&value. &year. sex level 1" bmis [-1, &ref.] [1, &value.] sex*bmis [-1, &ref. 1] [1, &value. 1]/ exp cl cov;
estimate "&value. &year. sex level 2" bmis [-1, &ref.] [1, &value.] sex*bmis [-1, &ref. 2] [1, &value. 2]/ exp cl cov

where sex is a class variable. This should be relatively easy to check since the interaction only matters for the non-reference level, so for the reference level of sex the two lines

estimate "&value. &year. sex level 2 (ref)" bmis [-1, &ref.] [1, &value.] sex*bmis [-1, &ref. 2] [1, &value. 2]/ exp cl cov
estimate "&value. &year." bmis [-1, &ref.] [1, &value.] / exp cl cov

should produce identical results, so you can just change which level is the reference and use the simpler code.

Re: How do I plot the spline effect in PHREG?

sseb — Thu, 24 Dec 2020 23:56:57 GMT

Original post:

@LinusS_ Thanks for the quick reply! Yes, I forgot to mention in my previous comment that I had added an interaction term between bmi and sex.

I tried the code you commented, but I'm getting an error saying 'the level 21 is not valid for CLASS variable sex. The level specifications for this variable must range from 1 to 2'. (Note that 21 is the reference I set for BMI). I assume this has something to do with the order of the estimate statement and SAS thinking the 21 is for sex and not the 1/2 at the end, but I don't know how to change the estimate statement to reflect this. Any ideas?

I tried adding a comma between &ref/&value and 1/2 and then tried switching the interaction term from sex*bmi to bmi*sex, but no luck. (By the way, if you have any resources on coding with the estimate statement with the square brackets, please let me know!)

Thanks again for all your help!

Edit:

Nevermind! Just found some SAS resources on non-positional syntax (p.455) and figured out that it should be

bmis*sex[-1, 1 &ref.][1, 1 &value.]

and not

bmis*sex[-1, &ref. 1][1, &value. 1]

(and so on for 2) because of the order of my variables. It works now. Thanks!

Re: How do I plot the spline effect in PHREG?

maomao0803 — Wed, 22 Sep 2021 19:49:36 GMT

Hi there,

I am a student using SAS for a project. I used the macro code you listed here to try to draw a cubic spline curve for cox model with multiple covariates adjusted. For unknown reasons, it didn't work and showed the following error. The file "estimates" was correctly generated and can be opened. Could you help with troubleshooting, please? Thank you very much!

"94 ods html select none;

ERROR: The HTML destination is not active; no select/exclude lists are available.

96 ods rtf select none;

ERROR: The RTF destination is not active; no select/exclude lists are available.

.......

108 ods html select all;

ERROR: The HTML destination is not active; no select/exclude lists are available.

109

110 ods rtf select all;

ERROR: The RTF destination is not active; no select/exclude lists are available.".

******here are the codes I used ***********

ods graphics on;

PROC PHREG DATA = BB.omega3_cvd ; /* or other plots */
effect spl = spline(diet_avg / basis=tpf(noint) NATURALCUBIC details knotmethod=rangefractions(1, 1.4, 2.0) );
model PERMTH_INT*MORTSTAT(0) = spl RIDAGEYR RIAGENDR RIDRETH1 DMDEDUC2 DMDMARTL PIR dietenergyavg active
smokestatus drinkstatus SDDSRVYR BMXBMI / rl=wald ties=EFRON ;
store BB.coxr;
run;

%macro est(ref=1.4, start=0, end=20, by=0.1);

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

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

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

%end;

%mend est;

ods html select none;

ods rtf select none;

ods dataset Estimates=Estimates;

proc plm restore=BB.coxr;

%est(ref=1.4, start=0, end=20, by=0.1);

run;

ods html select all;

ods rtf select all;

data estimates;

set estimates;

diet_avg=label*1;

run;

proc sgplot data=estimates NOAUTOLEGEND Noborder;

Series y=ExpEstimate x=diet_avg / ;

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

Series y=UpperExp x=diet_avg / LINEATTRS=(pattern=ShortDash color=Black THICKNESS=1);

REFLINE 1 / axis=y;

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

yaxis Values=(0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 3 4) Label="Hazard ratio" Type=LOG LABELATTRS=(weight=BOLD);

xaxis min=0 VALUES=(0 to 20 by 1) LABELATTRS=(weight=BOLD);

run;

Re: How do I plot the spline effect in PHREG?

sophiec — Sat, 29 Apr 2023 01:44:26 GMT

I have the same question! Did you ever figure out how to plot splines with PROC SURVEYPHREG?

Thanks!

Sophie

Re: How do I plot the spline effect in PHREG?

HenrikLJ — Wed, 09 Aug 2023 08:25:24 GMT

This is a great program. Thank you!