I'm new in fitting multivariable Cox models w/ restricted cubic splines.
This is my code:
d <- datadist(df)
options(datadist = "d")
mod_full <- cph(Surv(time, status) ~ rcs(age, 3) + sex + weight, data = df, x = TRUE, y = TRUE, surv = TRUE)
mod_full
summary(mod_full)
anova(mod_full)
However, I'm not quite sure how to interpret & report the model. In a model w/o RCS, I would just report HR + 95% CI + p. Would be very thankful for some insights.
Let's again start with the lung data from the {survival} package and fit a similar model. We will expand age in a restricted cubic spline and also model the effect of calories consumed at meals. The model is
library(rms)
library(survival)
dd <- datadist(lung)
options('datadist'=dd)
fit <- cph(Surv(time, status) ~ rcs(age, 3) + meal.cal, data=lung, x=T, y=T)
Let $X$ be the design matrix containing spline basis functions and let $W$ be the calories consumed at meals. The model for the hazard is
$$ \lambda(t; X, W) = \lambda_0(t) \exp(X\beta + W\gamma) $$
Here, $\lambda_0(t)$ is the baseline hazard, $\beta$ are the regression coefficients for the spline and $\gamma$ is the regression coefficient for calories.
Since we're using a spline basis with 3 knots for age, we should have 2 basis functions in $X$ meaning we have 3 total regression coefficients: 2 for age, 1 for calories. We can verify this by print fit to the console
> fit
Frequencies of Missing Values Due to Each Variable
Surv(time, status) age meal.cal
0 0 47
Cox Proportional Hazards Model
cph(formula = Surv(time, status) ~ rcs(age, 3) + meal.cal, data = lung,
x = T, y = T)
Model Tests Discrimination
Indexes
Obs 181 LR chi2 4.00 R2 0.022
Events 134 d.f. 3 R2(3,181)0.006
Center 0.5728 Pr(> chi2) 0.2611 R2(3,134)0.007
Score chi2 4.16 Dxy 0.108
Pr(> chi2) 0.2449
Coef S.E. Wald Z Pr(>|Z|)
age 0.0077 0.0218 0.35 0.7236
age' 0.0149 0.0257 0.58 0.5624
meal.cal 0.0000 0.0002 -0.07 0.9406
The coefficients presented here are the log hazard ratios and their standard errors. This is typically what you put in a table 2; it will allow other people to fit a similar model on their data. One could evaluate a hypothesis test for calories consumed from meals from this (in this case, its clear there is a negligible effect of calories on survival, if there is an effect at all). However, the spline is more complicated. The spline is the set of coefficients all together, so interpreting a hypothesis test for any one coefficient is meaningless.
Instead, we can use anova to perform a joint test for the spline coefficients.
> anova(fit)
Wald Statistics Response: Surv(time, status)
Factor Chi-Square d.f. P
age 3.80 2 0.1498
Nonlinear 0.34 1 0.5624
meal.cal 0.01 1 0.9406
TOTAL 4.12 3 0.2490
Its been a while since I've used rms so I hope someone can correct the following if I am wrong. The statistics reported here are likelihood ratio chi-square statistics Wald statistics estimated from the estimated covariance matrix. You can see that the chi square statistic and p value for calories is very similar to the wald statistic and p value printed out in the last step. There is a good reason for this I will not get into.
Age has 2 chi-square statistics. The first tests for any effect of age (i.e. $H_0: \beta = 0$ vs $H_A: \beta \neq 0$). The second tests of there is a non-linear effect of age by testing just the coefficients which would produce a non-linear relative hazard.
What do you report? Report both. You should report confidence intervals and point estimates for all coefficients in the model, but also report relevant tests of hypothesis from the anova call.
It would also be sensible to plot these effects while holding other variables constant, something we covered in this answer.
