1

I am using the coxme package in R to perform mixed effect cox regression analysis, and would greatly appreciate some advice on how to interpret the regression slopes. The data is a response variable (lifespan), converted to a survival object (age at death, and information on censoring). There are two sexes (M: male, and F: female) and two groups (G1 and G2), and sampling was performed in two separate years, and there is a continuous explanatory variable (NE). I use sex, group and NE as fixed effects, and year as a random effect, because the aim is to quantify the differences in response to NE ($\beta$), for the within-sex and with-group comparisons:

  • $\beta_{M,G1} - \beta_{F,G1}$
  • $\beta_{M,G2} - \beta_{F,G2}$
  • $\beta_{M,G1} - \beta_{M,G2}$
  • $\beta_{F,G1} - \beta_{F,G2}$

So here is some dummy data and the model (in R):

n = 2000
set.seed(1239)
dumDF = data.frame(
    "Sex" = sample(c("M","F"), size = n, replace = T), 
    "Group" = sample(c("G1","G2"), replace = T),
    "Year" = sample(c("1994","1995"), replace = T), 
    "NE" = rnorm(n,0,5), 
    "Status" = sample(c(1,2,2,2), replace = T), 
    "Life" = rnorm(n, 50, 5))

dumDF$Life = ifelse(dumDF$Sex=="M", dumDF$Life -10 + dumDF$NE, dumDF$Life+dumDF$NE)

dumDF$SurvObj = with(dumDF, Surv(Life, Status == 2))
coxdum <- coxme(SurvObj ~ Sex*Group*NE + (1|Year), data = dumDF)
coxdum

> coxdum
Cox mixed-effects model fit by maximum likelihood
  Data: dumDF
  events, n = 1000, 2000
  Iterations= 1 7 
                    NULL Integrated   Fitted
Log-likelihood -6597.047   -6104.89 -6104.89

                   Chisq df p    AIC    BIC
Integrated loglik 984.32  8 0 968.32 929.05
 Penalized loglik 984.32  7 0 970.32 935.96

Model:  SurvObj ~ Sex * Group * NE + (1 | Year) 
Fixed coefficients
                        coef exp(coef)   se(coef)      z      p
SexM             1.971525293 7.1816222 0.10196386  19.34 0.0000
GroupG2         -0.035924373 0.9647132 0.08891258  -0.40 0.6900
NE              -0.176770390 0.8379722 0.01356340 -13.03 0.0000
SexM:GroupG2     0.004847530 1.0048593 0.12708703   0.04 0.9700
SexM:NE         -0.053164686 0.9482238 0.01909755  -2.78 0.0054
GroupG2:NE      -0.001227879 0.9987729 0.01866811  -0.07 0.9500
SexM:GroupG2:NE  0.043730917 1.0447012 0.02699633   1.62 0.1100

Random effects
 Group Variable  Std Dev Variance
 Year  Intercept 2e-02   4e-04   

I then reconstruct the slopes from the estimated coefficients:

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

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

    matrix(c(FG1,MG1,MG2,FG2), ncol = 1)
    }
}
coxSlopeFunc(coxdum)

> coxSlopeFunc(coxdum)
             [,1]
[1,] -0.176770390
[2,] -0.229935076
[3,] -0.187432037
[4,] -0.177998268

As you can see, the slopes are all negative ([1,] through [4,]). It is understanding the meaning of these slopes that I am stuck on.

Would it be correct to say that a negative slope means there is a decreasing hazard (and therefore reduced mortality) with increasing NE? In other words, if the slope is negative, the explanatory variable reduces the rate of mortality?

Furthermore, if the value of NE increases by a unit of 1, is the hazard ($H$) expected to change by $H \times (1 + \beta)$?


Update:

As I read this I interpret the response to be the log of the individual hazard function, $log \space h_i(t)$, as given on page 3. Thus, a negative slope would indicate that the hazard (mortality) reduces with increasing NE.

rg255
  • 752
  • 2
  • 8
  • 27

0 Answers0