# Calculate Robust CI's for Stratified-Cox Regression Interaction

I have a multivariable stratified cox regression model with correction for correlation with Robust SE. There is interaction by the stratified variable, and so I would like to report adjusted Hazard Ratios with CIs for this interaction. Below is an example:

library(survival)
data("kidney")

age_cuts <- quantile(kidney$$age, probs = c(.33, .66)) kidney$$agegroup <- cut(kidney$$age, breaks = c(-Inf, age_cuts, Inf), labels = c('t1', 't2', 't3')) kidney$$age2 <- ifelse(kidney$$agegroup=='t2', 1, 0) kidney$$age3 <- ifelse(kidney$$agegroup=='t3', 1, 0) kidney$$disease2 <- ifelse(kidney\$disease %in% c('PKD', 'Other'), 1, 0)

mod <- coxph(Surv(time, status) ~ strata(agegroup) + disease2 +
disease2:age2 + disease2:age3 + frail + cluster(id),
data = kidney)
summary(mod)
coef exp(coef) se(coef) robust se      z Pr(>|z|)
disease2      -0.2028    0.8165   0.5625    0.3891 -0.521    0.602
frail          1.4463    4.2473   0.2720    0.2500  5.786 7.22e-09
disease2:age2 -0.6402    0.5272   0.7904    0.5612 -1.141    0.254
disease2:age3 -1.3371    0.2626   0.8380    0.8968 -1.491    0.136

disease2
frail         ***
disease2:age2
disease2:age3
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

exp(coef) exp(-coef) lower .95 upper .95
disease2         0.8165     1.2248   0.38085     1.750
frail            4.2473     0.2354   2.60215     6.933
disease2:age2    0.5272     1.8968   0.17550     1.584
disease2:age3    0.2626     3.8081   0.04528     1.523


The exponentiated coefficient for 'disease2' (with CIs from Robust SE) would be the Hazard Ratio for disease within the first age category: HR 0.82 (0.38-1.75).

For age categories 2 and 3, the HRs would be exponentiated coefficients disease2 + disease2:age2 and disease2 + disease2:age3, respectively. But how could I calculate the robust SE for this linear combination in order to report CIs? Is there a way to use the model vcov?

I'm aware of the emmeans package, but it's unclear how it handles interaction with a stratified term or robust SEs. So it might be ideal to calculate manually.

(EDIT)...would it be?

terms <- names(mod$$means) vcov <- mod$$var
var <- diag(vcov)
names(var) <- terms
dimnames(vcov) <- list(terms, terms)
sqrt(var["disease2"] + var["disease2:age2"] + 2*vcov["disease2", "disease2:age2"])
sqrt(var["disease2"] + var["disease2:age3"] + 2*vcov["disease2", "disease2:age3"])


For starters, the model code shown does not produce the summary that is shown. I obtained your summary after removing the term cluster(id) and putting it in a separate argument:

mod <- coxph(Surv(time, status) ~ strata(agegroup) + disease2 +
disease2:age2 + disease2:age3 + frail, cluster = id,
data = kidney)


That point resolved, let's examine the SEs obtained via the vcov() function:

> sqrt(diag(vcov(mod)))
disease2         frail disease2:age2 disease2:age3
0.3890715     0.2499800     0.5611917     0.8967877


and note that these are the same as what is labeled "robust se" in the summary.

Accordingly, we can obtain the needed estimates via:

> L <- data.frame(age1=c(1,0,0,0), age2 = c(1,0,1,0), age3 = c(1,0,0,1))
> bhat <- coef(mod)
> sapply(L, function(x) sum(x * bhat))
age1       age2       age3
-0.2027748 -0.8429491 -1.5399127


and their (robust) SEs via:

> V <- vcov(mod)
> sapply(L, function(x) sqrt(t(x) %*% V %*% x))
age1      age2      age3
0.3890715 0.3940111 0.8462397


The emmeans support for coxph objects does not seem to work right when there is a strata term in the model. However, you can manually construct the needed emmGrid object using the model estimates, plus the L matrix created above:

> library("emmeans")
> emm <- emmobj(bhat = coef(mod), V = vcov(mod),
+     levels = c("a1","a2","a3"), linfct = t(as.matrix(L)),
+     df = Inf, tran = "log")

> summary(emm)
level estimate    SE  df asymp.LCL asymp.UCL
a1      -0.203 0.389 Inf    -0.965    0.5598
a2      -0.843 0.394 Inf    -1.615   -0.0707
a3      -1.540 0.846 Inf    -3.199    0.1187

Results are given on the log (not the response) scale.
Confidence level used: 0.95

> ### exponentiated results
> summary(emm, type = "response")
level response    SE  df asymp.LCL asymp.UCL
a1       0.816 0.318 Inf    0.3809     1.750
a2       0.430 0.170 Inf    0.1989     0.932
a3       0.214 0.181 Inf    0.0408     1.126

Confidence level used: 0.95
Intervals are back-transformed from the log scale


Note that the CIs for age 1 match those for disease2 in the outputs in the summary.

## Update

I have repaired the support for coxph models in emmeans so that it can handle a strata() term (previously those types of models were not supported). The updated code will be in versions after 1.4.5, and are available from the github site (see the DESCRIPTION).

Here is how it works out for this example. First, create the agegroup and disease2 variables as in the question. (disease2 is a binary variable, and it acts like a factor in emmeans().) Fit this model:

mod2 <- coxph(Surv(time, status) ~ strata(agegroup) * disease2 +
frail + cluster(id), data = kidney)


You can verify that the summary of mod2 is virtually identical to what is shown in the OP. The emmeans() function predicts log hazard for all combinations of agegroup and disease2:

> emm <- emmeans(mod2, ~ disease2 | agegroup)
> emm
agegroup = t1:
disease2 emmean    SE  df asymp.LCL asymp.UCL
0  1.713 0.296 Inf    1.1325      2.29
1  1.510 0.448 Inf    0.6327      2.39

agegroup = t2:
disease2 emmean    SE  df asymp.LCL asymp.UCL
0  1.713 0.296 Inf    1.1325      2.29
1  0.870 0.475 Inf   -0.0618      1.80

agegroup = t3:
disease2 emmean    SE  df asymp.LCL asymp.UCL
0  1.713 0.296 Inf    1.1325      2.29
1  0.173 0.729 Inf   -1.2555      1.60

Results are given on the log (not the response) scale.
Confidence level used: 0.95


To get log hazard ratios, use comparisons of those:

> prs <- pairs(emm, reverse = TRUE)
> confint(prs, by = NULL)
contrast agegroup estimate    SE  df asymp.LCL asymp.UCL
1 - 0    t1         -0.203 0.389 Inf    -0.965    0.5598
1 - 0    t2         -0.843 0.394 Inf    -1.615   -0.0707
1 - 0    t3         -1.540 0.846 Inf    -3.199    0.1187

Results are given on the log (not the response) scale.
Confidence level used: 0.95


(Summarizing with by = NULL just makes for a more compact display.)

To get the hazard ratios themselves, back-transform:

> confint(prs, by = NULL, type = "response")
contrast agegroup ratio    SE  df asymp.LCL asymp.UCL
1 / 0    t1       0.816 0.318 Inf    0.3809     1.750
1 / 0    t2       0.430 0.170 Inf    0.1989     0.932
1 / 0    t3       0.214 0.181 Inf    0.0408     1.126

Confidence level used: 0.95
Intervals are back-transformed from the log scale