# Interpretting Cox Regression ANOVA

I'm having difficulty interpreting the results from anova() in the rms package. My confusion arises from what information the summary() function provides compared to anova().

My model results:

ddist <- datadist(df)
options(datadist="ddist")
S <- Surv(df$$OS.time, df$$OS)
fit <- cph(S ~ rcs(age, 4) + gender + rcs(TP53.rna, 4) + rcs(RET.rna, 4) + rcs(TP53.rna, 4) * rcs(RET.rna, 4))
fit
Cox Proportional Hazards Model

cph(formula = Surv(df$$OS.time, df$$OS) ~ rcs(age, 4) + gender +
rcs(TP53.rna, 4) + rcs(RET.rna, 4) + rcs(TP53.rna, 4) * rcs(RET.rna,
4), data = df)

Model Tests       Discrimination
Indexes
Obs       511    LR chi2     49.95    R2       0.094
Events    216    d.f.           19    Dxy      0.263
Center 7.7451    Pr(> chi2) 0.0001    g        0.534
Score chi2  59.35    gr       1.707
Pr(> chi2) 0.0000

Coef     S.E.    Wald Z Pr(>|Z|)
age                      0.0002  0.0217  0.01  0.9930
age'                    -0.0267  0.0610 -0.44  0.6610
age''                    0.3301  0.2746  1.20  0.2292
gender=MALE             -0.0664  0.1618 -0.41  0.6816
TP53.rna                 0.7518  2.2612  0.33  0.7395
TP53.rna'               -0.5187  4.0546 -0.13  0.8982
TP53.rna''               0.9551 25.7561  0.04  0.9704
RET.rna                  0.5396  4.5774  0.12  0.9062
RET.rna'                10.4676 11.4171  0.92  0.3592
RET.rna''              -54.9637 44.0895 -1.25  0.2125
TP53.rna * RET.rna      -0.0239  0.5413 -0.04  0.9648
TP53.rna' * RET.rna     -0.1546  0.9908 -0.16  0.8760
TP53.rna'' * RET.rna     0.9757  6.3974  0.15  0.8788
TP53.rna * RET.rna'     -1.2545  1.3585 -0.92  0.3558
TP53.rna' * RET.rna'     2.1108  2.5886  0.82  0.4148
TP53.rna'' * RET.rna'   -8.9685 17.0594 -0.53  0.5991
TP53.rna * RET.rna''     6.2989  5.2387  1.20  0.2292
TP53.rna' * RET.rna''   -8.7830  9.8230 -0.89  0.3713
TP53.rna'' * RET.rna''  33.4958 63.1735  0.53  0.5960
summary(fit)
Effects              Response : Surv(df$$OS.time, df$$OS)

Factor               Low     High   Diff.   Effect    S.E.    Lower 0.95 Upper 0.95
age                  53.4490 69.095 15.6450  0.017118 0.19751 -0.37000   0.40424
Hazard Ratio        53.4490 69.095 15.6450  1.017300      NA  0.69073   1.49820
TP53.rna              9.1969 11.032  1.8348 -0.343310 0.25449 -0.84209   0.15548
Hazard Ratio         9.1969 11.032  1.8348  0.709420      NA  0.43081   1.16820
RET.rna               4.4128  7.128  2.7152  0.027138 0.27591 -0.51363   0.56791
Hazard Ratio         4.4128  7.128  2.7152  1.027500      NA  0.59832   1.76460
gender - FEMALE:MALE  2.0000  1.000      NA  0.066374 0.16180 -0.25075   0.38350
Hazard Ratio         2.0000  1.000      NA  1.068600      NA  0.77822   1.46740

Adjusted to: TP53.rna=10.3185 RET.rna=5.7292
anova(fit)
Wald Statistics          Response: Surv(df$$OS.time, df$$OS)

Factor                                            Chi-Square d.f. P
age                                               23.29       3   <.0001
Nonlinear                                        12.18       2   0.0023
gender                                             0.17       1   0.6816
TP53.rna  (Factor+Higher Order Factors)           16.77      12   0.1583
All Interactions                                 11.75       9   0.2276
Nonlinear (Factor+Higher Order Factors)           7.72       8   0.4613
RET.rna  (Factor+Higher Order Factors)            16.98      12   0.1503
All Interactions                                 11.75       9   0.2276
Nonlinear (Factor+Higher Order Factors)          10.02       8   0.2636
TP53.rna * RET.rna  (Factor+Higher Order Factors) 11.75       9   0.2276
Nonlinear                                        11.54       8   0.1730
Nonlinear Interaction : f(A,B) vs. AB            11.54       8   0.1730
f(A,B) vs. Af(B) + Bg(A)                          3.01       4   0.5568
Nonlinear Interaction in TP53.rna vs. Af(B)       7.29       6   0.2949
Nonlinear Interaction in RET.rna vs. Bg(A)        6.97       6   0.3233
TOTAL NONLINEAR                                   29.02      14   0.0104
TOTAL NONLINEAR + INTERACTION                     29.02      15   0.0160
TOTAL                                             55.26      19   <.0001

According to my understanding of Cox Regression, summary() tells us that we cannot conclude whether any of the covariates are associated with survival due to all the confidence intervals including 1. Furthermore, we need to use summary to interpret results and not the raw model output because the cubic splines render p-values for individual coefficients uninformative.

Assuming that interpretation is correct, why do the anova results show age as highly significant when the confidence interval includes. Are summary() and anova() testing different hypotheses, and if so, what's the difference? The rms documentation states The anova function automatically tests most meaningful hypotheses in a design, but I'm having a hard time understanding which hypotheses are being tested, how to read the table, and how these differ from summary().

Any help on interpreting the anova table and model results is greatly appreciated!

## 1 Answer

The statistical issue here is the distinction between a global significance test and a test of differences between two specific scenarios.

Take-home message: you have to be careful not to over-generalize the results seen from tests of particular scenarios when continuous predictors are modeled nonlinearly.

The anova() test in the rms package is a global test of whether the indicated variables contribute overall to the model. The test reported in the first line of each predictor is effectively a test of whether there is a significant difference between models that contain or do not contain all terms involving the predictor; it's done via a Wald test on the coefficients for all such terms. Subsequent lines for each predictor similarly compare models including it but with or without interactions or nonlinear terms involving it. As explained here, this is different from the anova() test in the survival package, which by default adds predictors sequentially to the model in the order included in the formula and uses likelihood ratio tests.

The combination of coefficients including age was significant by this test, and so were the specifically nonlinear coefficients associated with age, so both age and its nonlinear terms qualify as "significant" in the rms anova() test.

The particular summary() function you used from rms compares two specific scenarios for each of the predictors at fixed values of the other predictors. From the documentation:

By default, inter-quartile range effects (odds ratios, hazards ratios, etc.) are printed for continuous factors, and all comparisons with the reference level are made for categorical factors.

So the result of summary() for age represents the difference between 53.4490 and 69.095 years of age. You are correct in that there is no significant difference in hazard between those ages.

The nonlinear modeling of age, however, allows for both observations. The coefficients for the nonlinear fit of age suggest that the hazard associated with age is not monotone. I have seen situations in which the hazard is relatively high for younger patients (presumably due to their having some more severe form of disease), drops for patients of intermediate ages, then rises again at advanced age. If that's the case in your data, then it might just be that the particular values of age chosen for the comparison in summary() happened not to differ in hazard.

You can examine this in your data. As you are using rms, you can use plot(Predict()) to see how the hazard changes with values of continuous predictors that were modeled with nonlinear terms.

• Thanks for the really helpful explanation! For clarification, how does anova() compare models with and without the variable to assess importance? I've read that methods like stepwise feature selection based on p-values are biased, so is it using some kind of AIC comparison? – Tomas Bencomo Apr 30 at 22:40
• @TomasBencomo it is a Wald test on the vector of regression coefficients involving the predictor. In this context for a single predictor without nonlinear terms or interactions it’s the Z-test as in the original printout of your fit for each of the individual terms, using the variance of each coefficient estimate to test whether the individual coefficient value is significantly different from 0. For a vector of coefficients it’s a test that all coefficients are different from 0, using both the individual variances and the covariances among them. – EdM May 1 at 2:18