Differences between Wald Statistics and P-values obtained with Dummy Coding vs. Direct Coding in Cox Models

Question

Using the larynx dataset (source: Survival Analysis Techniques for Censored and Truncated Data) to illustrate, supposing you want to estimates the hazard rates of event for stages 2 and stages 3 relative to stage 1, the analysis can proceed as follows:

> library(survival)
> data("larynx", package="KMsurv")
> 
> 
> #Recoding and reffiting model 
> larynx$Stg2vs1<- with(larynx,ifelse(stage == 2,1,
+                                     ifelse(stage==1,0,NA)))
> 
> larynx$Stg3vs1<- with(larynx,ifelse(stage == 3,1,
+                                     ifelse(stage==1,0,NA)))
> 
> mod1<- coxph(Surv(time,delta)~ Stg2vs1, data=larynx)
> summary(mod1)
Call:
coxph(formula = Surv(time, delta) ~ Stg2vs1, data = larynx)

  n= 50, number of events= 22 
   (40 observations deleted due to missingness)

       coef exp(coef) se(coef)     z Pr(>|z|)
Stg2vs1 0.07603   1.07899  0.45892 0.166    0.868

        exp(coef) exp(-coef) lower .95 upper .95
Stg2vs1     1.079     0.9268    0.4389     2.652

Concordance= 0.516  (se = 0.055 )
Rsquare= 0.001   (max possible= 0.948 )
Likelihood ratio test= 0.03  on 1 df,   p=0.9
Wald test            = 0.03  on 1 df,   p=0.9
Score (logrank) test = 0.03  on 1 df,   p=0.9

And for stage 3:

> mod2<- coxph(Surv(time,delta)~ Stg3vs1, data=larynx)
> summary(mod2)
Call:
coxph(formula = Surv(time, delta) ~ Stg3vs1, data = larynx)

  n= 60, number of events= 32 
   (30 observations deleted due to missingness)

          coef exp(coef) se(coef)    z Pr(>|z|)  
Stg3vs1 0.6115    1.8431   0.3556 1.72   0.0855 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 
1

    exp(coef) exp(-coef) lower .95 upper .95
Stg3vs1     1.843     0.5426    0.9181       3.7

Concordance= 0.606  (se = 0.048 )
Rsquare= 0.048   (max possible= 0.978 )
Likelihood ratio test= 2.95  on 1 df,   p=0.09
Wald test            = 2.96  on 1 df,   p=0.09
Score (logrank) test = 3.05  on 1 df,   p=0.08

However to using dummy coding gives the following:

> #Using Dummies
> 
> larynx$stage <- as.factor(larynx$stage)
> mod2 <- coxph(Surv(time,delta)~ stage, data=larynx)
> summary(mod2)
Call:
coxph(formula = Surv(time, delta) ~ stage, data = larynx)

  n= 90, number of events= 50 

          coef exp(coef) se(coef)     z Pr(>|z|)    
stage2 0.06481   1.06696  0.45843 0.141   0.8876    
stage3 0.61481   1.84930  0.35519 1.731   0.0835 .  
stage4 1.73490   5.66838  0.41939 4.137 3.52e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 
1

       exp(coef) exp(-coef) lower .95 upper .95
stage2     1.067     0.9372    0.4344      2.62
stage3     1.849     0.5407    0.9219      3.71
stage4     5.668     0.1764    2.4916     12.90

Concordance= 0.668  (se = 0.043 )
Rsquare= 0.167   (max possible= 0.987 )
Likelihood ratio test= 16.49  on 3 df,   p=9e-04
Wald test            = 19.24  on 3 df,   p=2e-04
Score (logrank) test = 22.88  on 3 df,   p=4e-05

Why are the estimates and p-values different?

Answer 1

They are different as your models for the instantaneous hazard are

\begin{aligned} \lambda_i(t) &= \lambda_{0}(t)\exp\left( \beta_2 1_{\text{is stage 2}}(i)\right) \\ \hat\lambda_i(t) &= \hat\lambda_{0}(t)\exp\left( \hat\beta_3 1_{\text{is stage 3}}(i)\right) \\ \tilde\lambda_i(t) &= \tilde\lambda_{0}(t)\exp\left( \tilde\beta_21_{\text{is stage 2}}(i) + \tilde\beta_31_{\text{is stage 3}}(i) + \tilde\beta_41_{\text{is stage 4}}(i)\right) \\ \end{aligned}

where the $1$ functions are an indicator which is one if individual $i$ has the state given group the subscript. Thus, you would only get the same slopes in the latter case if all the baseline hazards where equal for all four stages.