Cox regression hazard ratio (group 1 vs. group 4) changes calculation changes after combining group 2 and 3?

Question

I am working on a project which involves fitting a Cox proportional hazard model to a time to event situation in SAS. We have the variables:

1. y - event or censor (0/1)
2. timeto - time to event y
3. x - categorical variable that takes the value (1,2,3,4) which represents group 1,2,3,4

The SAS code I am using is:

proc phreg data=data1;
class x;
model timeto*y(0) = x / ties=efron rl;
hazardratio x / diff=all;
run;

I would then get the desired hazard ratios (HR) outputs for

• group 1 vs. 2
• group 1 vs. 3
• group 1 vs. 4 (this is the HR of interest, HR=0.5151, let's say)
• group 2 vs. 3
• group 2 vs. 4
• group 3 vs. 4

Then when I was asked to combine group 2 and 3 for simplicity purpose. So I did a simple data step:

data data1;
set data1;
if x=2 then x=3; *simply put all group 2 observations to group 3;
run;

Then when I run the proc phreg code again, I get the HR outputs

• group 1 vs. 3
• group 1 vs. 4 (this becomes different, HR=0.4988)
• group 3 vs. 4

From what I understand the HR for group 1 vs group 4 should not change based on the formula in https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.3/statug/statug_phreg_details24.htm.

I am not an expert, and I really wonder why the HR changes after regrouping. Might someone be willing to explain to me?

Answer 1

Following Therneau and Grambsch, Equation 3.4 of Modeling Survival Data--Extending the Cox Model (Springer, 2000), if there are no tied event times the score equation solved to estimate the parameter-value vector $\theta$ is:

$$ U(\theta) = \sum_{i=1}^n \int_0^{\infty} \left[X_i(s) - \bar x(\theta,s)\right] dN_i(s) = \sum_{i=1}^n U_i(\theta)$$

where $dN_i(s)$ is 1 if case $i$ has an event at time $s$ and 0 otherwise, $X_i$ represents the covariate values for case $i$ and $\bar x$ is a risk-weighted mean of $X$ over observations at risk:

$$\bar x(\theta,s) = \frac{\sum Y_i(s) r_i(s)X_i(s)}{\sum Y_i(s) r_i(s)}. $$

Here, $Y_i(s)$ is the at-risk indicator for cases at time $s$ and the risk-weighting is $r_i(\theta,s) =\exp[\theta' X_i(s)]$.

When you group 2 dummy variables into 1, you thus are changing the score equations solved to obtain the Cox model coefficients. With the nonlinearity introduced by the exponentiation needed to calculate the risk-weighted average, combining your original groups 2 and 3 together could change the corresponding risk-weighted covariate averages and thus the estimates of the other coefficients.