How to draw cubic splines after Cox regression in Stata 11.1?

[xiaoyanhuang]

Hi everyone,

In a cohort study, I would like to draw cubic splines, including HR and 95% CI, after Cox regression adjusted for age and sex. Yet, I have not found out the solution.

Suppose that there are variables as follows: observetime, censor, variablex (the independent variable we are interested in, continuous), age, sex. Knots would be set as 5.

Could you tell me how to perform it in Stata 11.1, please? I appreciate it very much!

Best regards,

Xiaoyan

 

[bukharin]

I'm going to assume that if censor is 1 they were censored, and if it's 0 they died. If it's the other way around you'll need to modify the code. First you need to set up your survival analysis with:

stset observetime, fail(censor==0)

Now you need to model the effect of variablex in a flexible way. If I were you I'd use fractional polynomials instead of cubic splines. It's much easier - all you need to do is:

fracpoly: stcox variablex age sex
fracplot

(what you're seeing on the y axis is the log hazard ratio)

If you really want to use cubic splines, one option would be to use the recently published -xblc- command. To install -xblc- use the following commands:

net sj 11-3 st0215_1
net install st0215_1

Here's an example of -xblc- using the cancer dataset that comes with Stata:

sysuse cancer, clear

* create spline variables for age
mkspline spl_age=age, cubic nknots(5)

* model effect of spline variables adjusting for treatment
stcox spl_age* i.drug

* calculate effect of age; we'll set the reference at 47 (youngest) but could also use 56 (mean)
quietly levelsof age, local(ages)
quietly xblc spl_age*, covname(spl_age1) at(`ages') reference(47) eform gen(ptage hr lb ub)

* graph the result
twoway (line hr ptage, sort) (line lb ub ptage, sort lc(black black) lp(- -)), ///
	legend(off) yscale(log) ylab(0.25 .5 1 2 4 8 16 32) ///
	xtitle(Age) ytitle(Hazard ratio) yline(1) ///
	title(Effect of age) subtitle(Adjusted for treatment)

 

[xiaoyanhuang]

Thank you very much, Bukharin.

You are right. The fractional polynomial is very easy to draw according to your suggestion. But is it different with a cubic spline?

Unfortunately, I still do not understand how to draw cubic splines. But thanks all the same for posting the example.

 

[epi_hsph]

Hi there,

I am using the xblc command after logistic (and Cox) regression with cubic splines. The xblc command has been working really well, until I added several variables to my adjusted models. The xblc command will not work with more than 7 variables! It works with any combination of 7, but returns an error message when I add an 8th variable.

Has anyone else had this problem?
Does anyone know of a solution?
Any insight would be much appreciated.

Thank you!

I'm going to assume that if censor is 1 they were censored, and if it's 0 they died. If it's the other way around you'll need to modify the code. First you need to set up your survival analysis with:

stset observetime, fail(censor==0)

Now you need to model the effect of variablex in a flexible way. If I were you I'd use fractional polynomials instead of cubic splines. It's much easier - all you need to do is:

fracpoly: stcox variablex age sex
fracplot

(what you're seeing on the y axis is the log hazard ratio)

If you really want to use cubic splines, one option would be to use the recently published -xblc- command. To install -xblc- use the following commands:

net sj 11-3 st0215_1
net install st0215_1

Here's an example of -xblc- using the cancer dataset that comes with Stata:

sysuse cancer, clear

* create spline variables for age
mkspline spl_age=age, cubic nknots(5)

* model effect of spline variables adjusting for treatment
stcox spl_age* i.drug

* calculate effect of age; we'll set the reference at 47 (youngest) but could also use 56 (mean)
quietly levelsof age, local(ages)
quietly xblc spl_age*, covname(spl_age1) at(`ages') reference(47) eform gen(ptage hr lb ub)

* graph the result
twoway (line hr ptage, sort) (line lb ub ptage, sort lc(black black) lp(- -)), ///
	legend(off) yscale(log) ylab(0.25 .5 1 2 4 8 16 32) ///
	xtitle(Age) ytitle(Hazard ratio) yline(1) ///
	title(Effect of age) subtitle(Adjusted for treatment)

 

[zh_zhang1984]

Dear

However, after execute fracplot, I find it is partial predictor+residual of sth. on Y-axis. what does that mean? how can I correlate it with hazard ratio?

I'm going to assume that if censor is 1 they were censored, and if it's 0 they died. If it's the other way around you'll need to modify the code. First you need to set up your survival analysis with:

stset observetime, fail(censor==0)

Now you need to model the effect of variablex in a flexible way. If I were you I'd use fractional polynomials instead of cubic splines. It's much easier - all you need to do is:

fracpoly: stcox variablex age sex
fracplot

(what you're seeing on the y axis is the log hazard ratio)

If you really want to use cubic splines, one option would be to use the recently published -xblc- command. To install -xblc- use the following commands:

net sj 11-3 st0215_1
net install st0215_1

Here's an example of -xblc- using the cancer dataset that comes with Stata:

sysuse cancer, clear

* create spline variables for age
mkspline spl_age=age, cubic nknots(5)

* model effect of spline variables adjusting for treatment
stcox spl_age* i.drug

* calculate effect of age; we'll set the reference at 47 (youngest) but could also use 56 (mean)
quietly levelsof age, local(ages)
quietly xblc spl_age*, covname(spl_age1) at(`ages') reference(47) eform gen(ptage hr lb ub)

* graph the result
twoway (line hr ptage, sort) (line lb ub ptage, sort lc(black black) lp(- -)), ///
	legend(off) yscale(log) ylab(0.25 .5 1 2 4 8 16 32) ///
	xtitle(Age) ytitle(Hazard ratio) yline(1) ///
	title(Effect of age) subtitle(Adjusted for treatment)

 

[ashley111]

What if your age variable has too many values to use? Are there any other alternative options?

Thanks!

 

[maartenbuis]

You can use postrcspline available from SSC. You will need to specify the link(log) option after stcox in order to get the odds ratios instead of the log(odds ratios).

 

[xiexf]

Hello
I analysed some clinical data and I used the command of "fractional polynomial" but I am not very sure about the result of y axis "partial predictor+residual of sth", is it the log HR of mortality? If so, what is the command if I specify the link(log) option after stcox in order to get the odds ratios instead of the log(odds ratios)?
thanks a lot

 

[bukharin]

Yes it's log(HR). My usual trick is not to change the y axis to HR, but to relabel it. For example, relabel the log HR of 0 as a HR of 1 (since log(1)=0).

webuse cancer
fracpoly: stcox age
fracplot, ylab(`=log(0.25)' "0.25" `=log(0.5)' "0.5" 0 "1" `=log(2)' "2" `=log(5)' "4") ytitle(Hazard ratio)

 

[xiexf]

Thanks very much for your reply.
I have some other questions.
First, as I adjusted several variables in the cox regression, when I used multiply fractional polynomials model to analysis, the shape of the graph changed, almost a straight line of HR appeared. I don’t know how to explain it, there must be some restriction when choose the fractional or multivariable fractional polynomials model. Would you please tell me?
Second, when I used the command of restrict cubic spline in the post, I got another different graph with very wide CIs range, that’s may because of the trick of yaxis HR instead of logHR. But I read a paper about the two methods “ fractional polynomials VS restric cubic spline” Actually, there is some difference between them, I want to consult how to choose an appropriate and right method in analysis. Or is there any rules in choosing these two methods?

Thanks a lot again.

 

[bukharin]

It's not too surprising that the shape would change when you add other variables, it will depend on the relationship between your other variables and the one of interest. Yes FP and RC splines are different techniques so will occasionally give you slightly different results. Did you try centering and scaling the variable before the RC spline transformation? Sometimes this makes a difference. Your subject matter expertise should also help to decide best fit (ie what kind of relationship would you expect between your variable and mortality?)

 

[xiexf]

Thank you, I think I didn’t stata my problems clearly, the first question I told yesterday is the models about “FP” and “MFP”. Please see the command and plot below, I adjust the same 9 variables in the two models, the graphs are different, so I am not sure which command FP or MFP is right for this case.
1. FP:
stset studytime, fail(died==1)
fracpoly: stcox sbp age var5 var6 var7 var8 var9 var10 var11
fracplot, addplot (histogram drug, frequency yaxis(2))
figure 1 View attachment 4594

2. MFP:
stset studytime, fail(died==1)
MFP: stcox sbp age var5 var6 var7 var8 var9 var10 var11
fracplot, addplot (histogram drug, frequency yaxis(2))
figure 2 View attachment 4595
I did not center and scale the variable before the RC spline transformation, I excuted the following command for RC, and the plot changed greatly.

stset studytime, fail(died==1)
mkspline spl_drug=drug, cubic nknots(5)
stcox spl_drug* age var5 var6 var7 var8 var9 var10 var11
quietly levelsof drug, local(drugs)
quietly xblc spl_drug*, covname(spl_drug1) at(`drugs') reference(1.922) eform gen(ptdrug hr lb ub)
twoway (line hr ptdrug, sort) (line lb ub ptdrug, sort lc(black black) lp(- -)), ///
legend(off) yscale(log) ylab(-0.25 .5 1 2 4 8 16 32) ///
xtitle(spb) ytitle(HR) yline(1) ///
figure 3 View attachment 4596

 

[abcstat]

Thanks for this discussion. I'm still having trouble understanding the y-axis. What do negative values mean? I see there are some negative values in the plots presented here, and my plot is entirely negative.