I am working on a dataset to test the association between empirical antibiotics (variable emp, the antibiotics are cefuroxime or ceftriaxone compared with a reference antibiotic) and 30-day mortality (variable mort30). The data comes from patients admitted in 6 hospitals (variable site2) with a specific type of infection. Therefore, I would like to adjust for this clustering of patients on hospital level.

First I did this using the coxme() function for mixed models. However, based on visual inspection of the Schoenfeld residuals there were violations of the proportional hazards assumption and I tried adding a time transformation (tt) to the model. Unfortunately, the coxme() does not offer the possibility for time transformations.

Therfore, I tried other options to adjust for the clustering, including coxph() combined with frailty() and cluster. Surprisingly, the standard errors I get using the cluster() option are much smaller than using the coxme() or frailty().

**Does anyone know what is the explanation for this and which option would provide the most reliable estimates? **

1) Using coxme:

> uni.mort <- coxme(Surv(FUdur30, mort30num) ~ emp + (1 | site2), data = total.pop)
> summary(uni.mort)
Cox mixed-effects model fit by maximum likelihood
  Data: total.pop
  events, n = 58, 253
  Iterations= 24 147 
                    NULL Integrated    Fitted
Log-likelihood -313.8427  -307.6543 -305.8967

                  Chisq   df         p  AIC  BIC
Integrated loglik 12.38 3.00 0.0061976 6.38 0.20
 Penalized loglik 15.89 3.56 0.0021127 8.77 1.43

Model:  Surv(FUdur30, mort30num) ~ emp + (1 | site2) 
Fixed coefficients
                    coef exp(coef)  se(coef)    z    p
empCefuroxime  0.5879058  1.800214 0.6070631 0.97 0.33
empCeftriaxone 1.3422317  3.827576 0.5231278 2.57 0.01

Random effects
 Group Variable  Std Dev   Variance 
 site2 Intercept 0.2194737 0.0481687
> confint(uni.mort)
                    2.5 %   97.5 %
empCefuroxime  -0.6019160 1.777728
empCeftriaxone  0.3169202 2.367543

2) Using frailty()

uni.mort <- coxph(Surv(FUdur30, mort30num) ~ emp + frailty(site2), data = total.pop)
> summary(uni.mort)
Call:
coxph(formula = Surv(FUdur30, mort30num) ~ emp + frailty(site2), 
    data = total.pop)

  n= 253, number of events= 58 

               coef   se(coef) se2    Chisq DF  p     
empCefuroxime  0.6302 0.6023   0.6010 1.09  1.0 0.3000
empCeftriaxone 1.3559 0.5221   0.5219 6.75  1.0 0.0094
frailty(site2)                        0.40  0.3 0.2900

               exp(coef) exp(-coef) lower .95 upper .95
empCefuroxime      1.878     0.5325    0.5768     6.114
empCeftriaxone     3.880     0.2577    1.3947    10.796

Iterations: 7 outer, 27 Newton-Raphson
     Variance of random effect= 0.006858179   I-likelihood = -307.8 
Degrees of freedom for terms= 2.0 0.3 
Concordance= 0.655  (se = 0.035 )
Likelihood ratio test= 12.87  on 2.29 df,   p=0.002

3) Using cluster()

uni.mort <- coxph(Surv(FUdur30, mort30num) ~ emp, cluster = site2, data = total.pop)
> summary(uni.mort)
Call:
coxph(formula = Surv(FUdur30, mort30num) ~ emp, data = total.pop, 
    cluster = site2)

  n= 253, number of events= 58 

                 coef exp(coef) se(coef) robust se     z Pr(>|z|)    
empCefuroxime  0.6405    1.8975   0.6009    0.3041 2.106 0.035209 *  
empCeftriaxone 1.3594    3.8937   0.5218    0.3545 3.834 0.000126 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

               exp(coef) exp(-coef) lower .95 upper .95
empCefuroxime      1.897     0.5270     1.045     3.444
empCeftriaxone     3.894     0.2568     1.944     7.801

Concordance= 0.608  (se = 0.027 )
Likelihood ratio test= 12.08  on 2 df,   p=0.002
Wald test            = 15.38  on 2 df,   p=5e-04
Score (logrank) test = 10.69  on 2 df,   p=0.005,   Robust = 5.99  p=0.05

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).
>