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!
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!
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;
Hi, can you please tell me what is matrix package? Is it a R package?
@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
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 - Desquilbe...)?
Thanks!
Tessa
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;
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!
Xi
I have the same question! Did you ever figure out how to plot splines with PROC SURVEYPHREG?
Thanks!
Sophie
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;
@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.
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!
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;
******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;
This is a great program. Thank you!
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 16. Read more here about why you should contribute and what is in it for you!
Need courses to help you with SAS Life Sciences Analytics Framework, SAS Health Cohort Builder, or other topics? Check out the Health and Life Sciences learning path for all of the offerings.