2
$\begingroup$

I have been trying to understand a bit better polynomial contrasts using this resource/example : https://library.virginia.edu/data/articles/understanding-ordered-factors-in-a-linear-model I know this is for linear model but I was trying to follow the same steps/logic using Cox regressions for time-to-event outcomes. My understanding from the example on the webpage above (please correct me if I am wrong) is that for linear models, using an ordered factor is equivalent to using a polynomial of degree # of categories -1 and that both approaches lead to the exact same p-values despite different coefficients/SE (that I assume have to do with different design matrices when using poly() vs ordered()).

I have been trying to apply this to an example of Cox regression, and in that case, the two approaches do not seem to be equivalent. Here is a reproducible example:

library(survival)
data(lung)
lung <- lung[complete.cases(lung),]

lung$surv_obj <- with(lung, Surv(time, status == 2))
lung$ph.ecog_o <- ordered(lung$ph.ecog)

cox1 <- coxph(surv_obj~ph.ecog_o,data=lung)
cox2 <- coxph(surv_obj~poly(unclass(ph.ecog_o), 3),data=lung)

I am not able to insert images, but the results are for cox1:

Call:
coxph(formula = surv_obj ~ ph.ecog_o, data = lung)

  n= 167, number of events= 120 

               coef exp(coef) se(coef)     z Pr(>|z|)  
ph.ecog_o.L 1.57386   4.82525  0.69609 2.261   0.0238 *
ph.ecog_o.Q 0.48669   1.62692  0.52681 0.924   0.3556  
ph.ecog_o.C 0.08127   1.08467  0.27246 0.298   0.7655  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

            exp(coef) exp(-coef) lower .95 upper .95
ph.ecog_o.L     4.825     0.2072    1.2331    18.881
ph.ecog_o.Q     1.627     0.6147    0.5794     4.569
ph.ecog_o.C     1.085     0.9219    0.6359     1.850

Concordance= 0.615  (se = 0.028 )
Likelihood ratio test= 13.65  on 3 df,   p=0.003
Wald test            = 15.48  on 3 df,   p=0.001
Score (logrank) test = 16.99  on 3 df,   p=7e-04

While for cox2:

Call:
coxph(formula = surv_obj ~ poly(unclass(ph.ecog_o), 3), data = lung)

  n= 167, number of events= 120 

                                coef exp(coef) se(coef)     z Pr(>|z|)    
poly(unclass(ph.ecog_o), 3)1  4.2891   72.9001   1.2037 3.563 0.000366 ***
poly(unclass(ph.ecog_o), 3)2  1.2843    3.6122   1.1649 1.103 0.270224    
poly(unclass(ph.ecog_o), 3)3  0.3106    1.3643   1.0413 0.298 0.765477    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                             exp(coef) exp(-coef) lower .95 upper .95
poly(unclass(ph.ecog_o), 3)1    72.900    0.01372    6.8886    771.48
poly(unclass(ph.ecog_o), 3)2     3.612    0.27684    0.3683     35.42
poly(unclass(ph.ecog_o), 3)3     1.364    0.73299    0.1772     10.50

Concordance= 0.619  (se = 0.028 )
Likelihood ratio test= 13.65  on 3 df,   p=0.003
Wald test            = 15.48  on 3 df,   p=0.001
Score (logrank) test = 16.99  on 3 df,   p=7e-04

What leads to the discordances in p-values in that context that is not the case in LM? Let's say one would want to test linear or quadratic trends in ordered categorical variables for which we do not have the numerical raw data, which of the 2 models would be more accurate and why?

Thanks in advance

Improve this question
$\endgroup$
1
  • $\begingroup$ Thank you for raising that. I have corrected the call (I tested many stuff and had copied the wrong call but the summary is the one with a cubic polynomial). I am still unsure I understand why is it that in the linear example (from the website) the p-values associated with each parameter are the same with both poly and ordered but not here for Cox? $\endgroup$ Commented Feb 20 at 20:04

1 Answer 1

Reset to default
1
$\begingroup$

What leads to the discordances in p-values in that context that is not the case in LM?

The third coefficient (cubic) has the same p-value for both models; the difference is in what gets assigned to the first two terms (linear and quadratic).

The difference probably comes from the lack of an intercept in a Cox survival regression model (although I haven't worked through all the details). In a Cox model, the baseline survival function (which isn't even estimated directly) is the logical equivalent of the intercept in a linear regression model. That means that the model matrix differs between a coxph() and an lm() model. To see that, try the following commands:

head(model.matrix(cox1))
head(model.matrix(cox2))

If you look more closely into the model matrices, you will find that there are only two unique values for ph.ecog_o.Q in model cox1, but 4 for the corresponding poly(unclass(ph.ecog_o), 3)2 in model cox2 (one for each of the 4 distinct values of ph.ecog). In the models from the cited web page, both models have only 2 unique values for the quadratic terms (e.g., in poly(unclass(powerF),3) on that page).

What's critical is that both models overall are identical (plus or minus some rounding error in the concordance estimates).

Let's say one would want to test linear or quadratic trends in ordered categorical variables for which we do not have the numerical raw data, which of the 2 models would be more accurate and why?

First, it depends on what's underlying the ordered factor. Quoting from the cited web page:

Fox and Weisberg (2019) tell us contr.poly() is appropriate when the levels of your ordered factor are equally spaced and balanced... If we have an ordered factor that does not have equal spacing and/or is not balanced, we’re better off performing polynomial regression via the poly() function.

If you have more than this handful of levels of an ordered predictor, you would be even better off with a flexible spline fit. Polynomial fits force a single polynomial to cover the entire range of data. A spline lets the data help tell you the function form of the association between a predictor and outcome. See Section 2.4 of Frank Harrell's Regression Modeling Strategies.

Second, I'd be very skeptical of any Cox model that returns an apparent hazard ratio (HR) of 72.9, as you get with cox2. That would make me trust cox1 more (although again I haven't figured out the details of what's going on).

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you very much for your detailed thoughts. I was also wary of the HR of 72. $\endgroup$ Commented Feb 21 at 6:26
  • $\begingroup$ I just stumbled across this: poly() returns normalized columns, i.e. colSums(model.matrix(cox2)^2)= (1, 1, 1). Therefore the sd in the column shrinks with growing n and so the coefficient grows... $\endgroup$ Commented Nov 10 at 11:04

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.