Standard error of reconstructed slopes (coxph/coxme)

Question

I have performed mixed effect Cox hazard regressions, and reconstructed the slopes to get group specific slopes (e.g. sex-specific responses to the explanatory variable). I aim to test whether the slopes differ from one another (e.g. do males and females respond differently to the explanatory variable?). To do this I will use Z-tests (here and here) where

$$Z= \frac{\beta_1-\beta_2}{\sqrt{{SE_{\beta_1}}^{2}+{SE_{\beta_2}}^2}}$$

However, I have performed my models in R using the coxme package which gives the following output, from which I reconstruct the sex- and group-specific slopes with the included function.

...
Fixed coefficients
                        coef exp(coef)   se(coef)      z    p
SexM             0.091305017 1.0956031 0.09085235   1.00 0.31
GroupG2         -0.036313825 0.9643376 0.08889039  -0.41 0.68
NE              -0.192009224 0.8252993 0.01317388 -14.57 0.00
SexM:GroupG2     0.009757875 1.0098056 0.12750426   0.08 0.94
SexM:NE         -0.212264676 0.8087506 0.02008058 -10.57 0.00
GroupG2:NE      -0.006933708 0.9930903 0.01814987  -0.38 0.70
SexM:GroupG2:NE  0.044999019 1.0460268 0.02756553   1.63 0.10
...


coxSlopeFunc = function(model, nfixed = 1){
    if(nfixed ==1){
        # Slope for Females + G1
        FG1 = model$coefficients[3]
 # Slope for Males + G1
 MG1 = model$coefficients[3] + model$coefficients[5]

        # Slope for Females + G2
        FG2 = model$coefficients[3] + model$coefficients[6]
        # Slope for Males + G2
        MG2 = model$coefficients[3] + model$coefficients[5] + model$coefficients[6] + model$coefficients[7]

        # Sex differences in slope
        SG1 = FG1 - MG1
        SG2 = FG2 - MG2

    matrix(c(FG1,MG1,FG2,MG2,SG1,SG2), ncol = 1, byrow = T)}
}
round(coxSlopeFunc(coxdum),3)

> round(coxSlopeFunc(coxdum),3)
       [,1]
[1,] -0.192
[2,] -0.404
[3,] -0.199
[4,] -0.366
[5,]  0.212
[6,]  0.167

However, I am unsure how to calculate the SE of the slope for each - should I just sum the standard errors of the components?

$$\frac{(-0.192009224 - (-0.192009224 + -0.212264676))}{\sqrt{0.01317388^2 + (0.01317388 + 0.02008058)^2}}$$

Answer 1

In this question and in two other recent questions (here and here) you are interested in the "slope" that represents the relation of a continuous variable, NE, to survival in a Cox proportional hazards model. You are particularly interested in whether the relation of NE to survival differs between sexes; furthermore, individuals might belong to either of 2 Groups, G1 and G2. I'll assume that linearity with respect to NEand the proportional hazards assumption are both verified. I'll ignore the random-effects term* that you included in other related questions.

Much depends on how you want to treat the group membership issue. In the example of this question,** where there is little evidence for differences in outcome with respect to Group G1 versus G2, preliminary data exploration and subject-matter knowledge might suggest that you simply remove the Group variable from your model. Tests of the interaction term SexM:NE would then provide the answer to your question about sex differences with respect to the relation of NE to survival.

If you wish to maintain the breakdown by Group, you might consider using the anova wrapper for the Cox model output in R to provide a single test of the SexM:NE coefficient, combining information from both Groups. The anova function performs a hierarchical test of coefficients in the order of entry into the model, using the equivalent of Type I sums of squares. If the data are reasonably balanced among combinations of covariates (for survival models, in terms of event numbers) then this may provide a useful test provided that you are clear about what it examines. For example, if you specified the following model:

Surv(time,status)~ NE + Sex + Group + NE:Sex + NE:Group + Sex:Group + NE:Sex:Group

then anova would first associate as much as possible with NE, then with Sex, then with Group, and then test whether the NE:Sex interaction significantly explained any residual. Note that this is a different way of evaluating the results of the model from the treatment-contrasts summaries provided by print or summary that you display above, even though it is based on the same model.

If you have a particular interest in the 3-way interaction among NE, Sex and Group then you will need to examine contrasts of the coefficients, similar to what you propose but with an important difference in implementation. The test you propose is equivalent to a Wald test. Examining whether there was a difference between males and females with respect to the value of the NE coefficient in Group G2, you would test whether $\beta_{SexM:NE}+\beta_{SexM:GroupG2:NE}$ (MG2-FG2 in your question) is different from 0. As the coefficient estimates are correlated, you need the variance of a sum of correlated variables, which in this example is: $$ \mathop{\rm var}(\beta_{SexM:NE})+\mathop{\rm var}(\beta_{SexM:GroupG2:NE}) + 2\mathop{\rm cov}(\beta_{SexM:NE},\beta_{SexM:GroupG2:NE})$$

Your proposed Z-test (based on the square root of the variance) ignores the covariance of the estimates of the coefficients, which is necessary and is provided by the corresponding off-diagonal element of the variance-covariance matrix, which you can get from the vcov function applied to your model. If you use the rms package, then there is a contrast function for that package's cph Cox models that allows tests of arbitrary contrasts, including bootstrap non-parametric tests.

One warning on this approach, however: Terry Therneau, responsible for much of the survival analysis infrastructure in R, has warned in this vignette that the apparent similarity of Cox regression models to linear regression models does not necessarily extend to tests on contrasts of coefficients. Examine those arguments carefully as you proceed.

---

*In another question you note that you are also evaluating a "random effect" that only has 2 levels. Treating a variable with so few levels as a random effect can be considered inadvisable.

**The model results in this question differ from those in other questions you have asked recently on the same general matter. It helps to specify the seed for the random-number generator, e.g., with set.seed() in R, before you generate your random data, to have reproducible "random" data for demonstration.