How to create a Cubic Spline Curve in survival analyse?

[bhr-q]

Hello All,

I wanted to create a Cubic Spline Curve with a continuous variable on the x-axis and hazard ratio for the time to event variable on y-axis with 95% CI, I was trying to use proc phreg using the effect statement, but not sure about the code, could you please help me to do that? 

Thank you,

[Ksharp]

You could use EFFECT statement of PROC PHREG to get job done.
Check Rick blogs:
https://blogs.sas.com/content/iml/2024/06/03/vize-multivar-regression-splines.html
https://blogs.sas.com/content/iml/2019/02/18/regression-restricted-cubic-splines-sas.html
https://blogs.sas.com/content/iml/2017/04/19/restricted-cubic-splines-sas.html

[bhr-q]

Thank you so much for your answer. I did try the following way:

proc phreg data=tmp;
effect spl = spline(ind_var  / knotmethod=percentilelist(5));
model time_to_hfhosp * event(0) = spl;
store out=PhregStore;
run;
proc plm restore=PhregStore;
score data=tmp   out=SplPlot / ilink ; 
run;
proc sgplot data=SplPlot;
band x=ind_var     lower=lower_ci upper=upper_ci / transparency=0.5;
series x=ind_var   y=Predicted / lineattrs=(color=blue thickness=2);
run;

However, for the last part of the code, the log shows the following message:

"ERROR: Variable lower_ci not found.
ERROR: Variable upper_ci not found." 

So, I fixed it as shown below and got the following plot. but, I’m not sure if it’s a correct way, as I couldn’t find any sample code that uses the EFFECT statement within PROC PHREG?

proc phreg data=tmp;
effect spl = spline(ind_var/ naturalcubic knotmethod=percentiles(5));
model time_to_hfhosp * event(0) = spl;
store out=PhregStore;
run;
data ScoreIn;
do ind_var = 20 to 100 by 1;   output;   end;
run;
proc plm restore=PhregStore;
score data=ScoreIn out=pred_values predicted=log_hazard stderr=se_log_hazard;
run;
data plot;
set pred_values;
hazard_ratio = exp(log_hazard);
lower_ci = exp(log_hazard - 1.96 * se_log_hazard);
upper_ci = exp(log_hazard + 1.96 * se_log_hazard);
run;
proc sgplot data=plot;
band x=ind_var lower=lower_ci upper=upper_ci / transparency=0.5;
series x=ind_var y=hazard_ratio;
refline 1 / axis=y lineattrs=(pattern=shortdash);
run;

[Ksharp]

I think you should use HazzardRatio statement to calculate it .

data bmt;
 set sashelp.bmt;
 call streaminit(123);
 ind_var=ceil(rand('normal',50,10));
run;

ods output  HazardRatios=HazardRatios;
proc phreg data=bmt;
effect spl = spline(ind_var/ degree=3);
model T * Status(0) = spl;
hazardratio 'Effect of 1-unit change in ind_var' ind_var / at(ind_var=(40 to 60 by 1));
run;

data HazardRatios2;
 set HazardRatios;
 ind_var=input(scan(Description,-1,'='),best.);
run;

proc sgplot data=HazardRatios2;
band x=ind_var lower=WaldLower upper=WaldUpper / transparency=0.5;
series x=ind_var y=HazardRatio;
refline 1 / axis=y lineattrs=(pattern=shortdash);
run;

[bhr-q]

Thank you for your answer. I think for a cubic spline curve, the Y-axis should show the predicted HR across the range of the independent variable, not the raw HR, and your way shows the raw HR for a 1% increase in the independent variable at each level. 

I was thinking of using PROC PLM, but not sure about my code.

[Ksharp]

I think it is hard to use PROC PLM to get HazzardRatio, that why the reason sas showed in LOG:
NOTE: Hazard ratios that cannot be conveniently calculated or displayed are set to missing in the ParameterEstimates table. Use the HAZARDRATIO
statement to compute the needed hazard ratios.

Anyway, you could check LCLM UCLM option of PROC PLM:
https://blogs.sas.com/content/iml/2019/02/11/proc-plm-regression-models-sas.html

And I don't understand where is different from predicted HR and raw HR. If the ind_var value does not exist in raw data,that would lead to predicted HR,right?
Or maybe other stat expert would answer your question ?
Like: @StatDave @SteveDenham @lvm .....

[bhr-q]

Thank you for your answer. I meant there is a difference between HR and predicted HR, and I was thinking the spline curve might include the predicted HR on the y-axis rather than HR. I have one independent variable in the model, named RV.

I also checked the LCLM and UCLM options of PROC PLM, and I could run the code below:

But the below plot is a bit strange

proc phreg data=tmp;
effect spl = spline(RV / knotmethod=percentilelist(5));
model time_to_hfhosp * event(0) = spl;
store out=PhregStore;  
run;
proc plm restore=PhregStore;
score   data=tmp out=NewScore   predicted LCLM UCLM / ilink; 
run;
proc sgplot data=NewScore  noautolegend;
band x=RV lower=LCLM upper=UCLM / transparency=0.5;
series x=RV y=Predicted / lineattrs=(color=blue thickness=2);
refline 1 / axis=y lineattrs=(pattern=shortdash color=gray);  
run;

[Ksharp]

"But the below plot is a bit strange"
Yes. That is the reason why I suggest to use HAZARDRATIO to calcuated it. HR is not easy to obtain.
About your other question,I have no idea about it. Maybe other sas user could help you.