Mixed effect Cox models using coxme package

References

Load packages

library(coxme)

Data

data(package = "coxme")
data(eortc)

Null model

## Model without random effects
efit0 <- coxph(Surv(y, uncens) ~ 1, eortc)
summary(efit0)

## Call:  coxph(formula = Surv(y, uncens) ~ 1, data = eortc)
## 
## Null model
##   log likelihood= -10639 
##   n= 2323

Analysis without random effects

## Model without random effects
efit1 <- coxph(Surv(y, uncens) ~ trt, eortc)
summary(efit1)

## Call:
## coxph(formula = Surv(y, uncens) ~ trt, data = eortc)
## 
##   n= 2323, number of events= 1463 
## 
##       coef exp(coef) se(coef)    z Pr(>|z|)    
## trt 0.6183    1.8557   0.0634 9.75   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##     exp(coef) exp(-coef) lower .95 upper .95
## trt      1.86      0.539      1.64       2.1
## 
## Concordance= 0.563  (se = 0.006 )
## Rsquare= 0.044   (max possible= 1 )
## Likelihood ratio test= 106  on 1 df,   p=0
## Wald test            = 95  on 1 df,   p=0
## Score (logrank) test = 98  on 1 df,   p=0

## Log likelihood
efit1$loglik

## [1] -10639 -10586

## Chi-squared value comparing this model to the null model
-2 * (efit0$loglik - logLik(efit1))

## 'log Lik.' 105.7 (df=1)

Analysis with random effects

## Model with random effects
efit2 <- coxme(Surv(y, uncens) ~ trt + (1|center), eortc)
efit2

## Cox mixed-effects model fit by maximum likelihood
##   Data: eortc
##   events, n = 1463, 2323
##   Iterations= 9 49 
##                  NULL Integrated Fitted
## Log-likelihood -10639     -10521 -10479
## 
##                   Chisq    df p   AIC   BIC
## Integrated loglik 236.1  2.00 0 232.1 221.5
##  Penalized loglik 319.7 28.69 0 262.4 110.7
## 
## Model:  Surv(y, uncens) ~ trt + (1 | center) 
## Fixed coefficients
##       coef exp(coef) se(coef)     z p
## trt 0.7086     2.031  0.06424 11.03 0
## 
## Random effects
##  Group  Variable  Std Dev Variance
##  center Intercept 0.3292  0.1084

## Log likelihood
eft2logLik <- efit2$loglik + c(0, 0, efit2$penalty)
eft2logLik

##       NULL Integrated  Penalized 
##     -10639     -10521     -10479

## Chi-squared value comparing the penalized model to the null model
-2 * (eft2logLik["NULL"] - eft2logLik["Penalized"])

##  NULL 
## 319.7

Testing for random effects

Using integrated log likelihood

## -2 logLik difference
(logLikDiffNeg2 <- -2 * (logLik(efit1) - eft2logLik["Integrated"]))

## 'log Lik.' 130.5 (df=1)

## Degree of freedom
## efit1: df = 1
## efit2: df = 2
efit2$df[1]

## [1] 2

## Degree of freedom difference
(dfDiff <- efit2$df[1] - 1)

## [1] 1

## P value for the random effects: significant random effects across centers
pchisq(q = as.numeric(logLikDiffNeg2), df = dfDiff, lower.tail = FALSE)

## [1] 3.261e-30

Using penalized log likelihood

## -2 logLik difference
(logLikDiffNeg2 <- -2 * (logLik(efit1) - eft2logLik["Penalized"]))

## 'log Lik.' 214.1 (df=1)

## Degree of freedom
## efit1: df = 1
## efit2: penalized degrees of freedom (1 fixed effect and df for random effects (<= 37 centers))
efit2$df[2]

## [1] 28.69

## Degree of freedom difference
(dfDiff <- efit2$df[2] - 1)

## [1] 27.69

## P value for the random effects: significant random effects across centers
pchisq(q = as.numeric(logLikDiffNeg2), df = dfDiff, lower.tail = FALSE)

## [1] 1.036e-30

Simplification using a function

## Function to test random effects
TestRanef <- function(coxphModel,coxmeModel) {

    if (!class(coxphModel) == "coxph" | !class(coxmeModel) == "coxme") {
        stop("Wrong models")
    }

    ## Degrees of freedom
    coxphDf <- sum(!is.na(coef(coxphModel)))
    coxmeDf <- coxmeModel$df
    names(coxmeDf) <- c("Integrated","Penalized")
    ## DF differnces
    dfDiff <- coxmeDf - coxphDf


    ## Log likelihodds
    coxphLogLik <- coxphModel$loglik[2]
    coxmeLogLik <- coxmeModel$loglik + c(0, 0, coxmeModel$penalty)
    coxmeLogLik <- coxmeLogLik[2:3]
    ## -2 logLik difference
    logLikNeg2Diff <- c(-2 * (coxphLogLik - coxmeLogLik["Integrated"]),
                        -2 * (coxphLogLik - coxmeLogLik["Penalized"]))

    ## p-values
    pVals <- pchisq(q = logLikNeg2Diff, df = dfDiff, lower.tail = FALSE)

    ## Combine
    outDf <- data.frame(dfDiff, logLikNeg2Diff, pVals)
    colnames(outDf) <- c("df diff", "-2logLik diff", "p-values")
    outDf
}

TestRanef(efit1, efit2)

##            df diff -2logLik diff  p-values
## Integrated    1.00         130.5 3.261e-30
## Penalized    27.69         214.1 1.036e-30