0
$\begingroup$

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!

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Tomas Bencomo Apr 30 at 22:40
  • $\begingroup$ @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. $\endgroup$ – EdM May 1 at 2:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.