0
$\begingroup$

Based on the results table of the existing CoxPH model, is it possible to calculate the standard error when using the relevel function?

Like this...

when 2 vs. 1(ref),

exp(coef) 0.96712/0.97700 = 0.9898874

se(coef) ??/?? = ??


> coxph(Surv(time, status) ~ Trt, data = data)
Call:
coxph(formula = Surv(time, status) ~ Trt, data = data)

         coef exp(coef) se(coef)      z     p
Trt1 -0.02327   0.97700  0.02451 -0.949 0.342
Trt2 -0.03343   0.96712  0.02856 -1.171 0.242

Likelihood ratio test=1.48  on 2 df, p=0.4777
n= 10000, number of events= 10000 
> data$Trt <- relevel(data$Trt, ref = '1')
> coxph(Surv(time, status) ~ Trt, data = data)
Call:
coxph(formula = Surv(time, status) ~ Trt, data = data)

         coef exp(coef) se(coef)      z     p
Trt0  0.02327   1.02354  0.02451  0.949 0.342
Trt2  "Unknown" "Unknown"  "Unknown"  "Unknown"  "Unknown"

Likelihood ratio test=1.48  on 2 df, p=0.4777
n= 10000, number of events= 10000 
$\endgroup$

1 Answer 1

0
$\begingroup$

I don't know what's going on with relevel() here (and the answer is probably a software-specific matter off-topic on this site). There's a more generally useful way to get the standard errors of combinations of coefficient estimates in a Cox or other regression model, which doesn't require refitting.

The formula for the variance of a weighted sum of variables is the key. I find it least error-prone to work in the original coefficient scale and leave the transformation to hazard ratios until the end. For your difference between Trt2 and Trt1 with corresponding coefficients $\beta_2$ and $\beta_1$, the variance estimate is:

$$\text{Var}(\beta_2-\beta_1) = \text{Var} (\beta_2) + \text{Var} (\beta_1) -2\text{Cov}(\beta_1,\beta_2),$$

where $\text{Cov}$ is the covariance between the estimates. The vcov() function applied to the model output returns the variance-covariance matrix for the coefficient estimates. The individual variances will be the squares of the se(coef) values, but you'll need to use vcov() to get the covariance.

Once you get the variance for the difference, the square root is the standard error. The point estimate of $(\beta_2-\beta_1)$, plus/minus 1.96 times the standard error, gives the standard 95% confidence interval in the coefficient scale for a Cox model. Exponentiation puts those into the hazard-ratio scale.

Working with the single original model, instead of trying to refit with releveled multi-level factors, is a general strategy that extends better to more complex models than this. Post-modeling tools, for example those in the R car and emmeans packages, are designed to help with this.

$\endgroup$
1
  • $\begingroup$ Thanks for answering! I will try your suggestion. $\endgroup$
    – Chan
    Apr 9, 2023 at 2:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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