How to export coxme regression output to Latex?

Question

I have the results of a regression from a coxme model in R, and would like to export it as a Latex table. How do I do this?

Code:

> library(coxme)
>
> m1 <- coxme(Surv(exit, status) ~ 1 + (1 | region/id) + air_quality, df)
> summary(m1)

Cox mixed-effects model fit by maximum likelihood
  Data: df
  events, n = 163, 10738
  Iterations= 29 128 
                    NULL Integrated    Fitted
Log-likelihood -1390.143  -1287.026 -1005.886

                   Chisq     df p    AIC     BIC
Integrated loglik 206.23   3.00 0 200.23  190.95
 Penalized loglik 768.51 180.64 0 407.24 -151.61

Model:  Surv(exit, status) ~ 1 + (1 | region/id) + air_quality 
Fixed coefficients
                   coef exp(coef)  se(coef)    z    p
air_quality  -0.0535707 0.9478389 0.5251661 -0.1 0.92

Random effects
 Group           Variable    Std Dev  Variance
 region/id (Intercept) 2.850487 8.125278
 region          (Intercept) 1.205821 1.454004

I am aware that stargazer doesn't support coxme objects, and I'd ideally like the functionality it offers. I only wish to report the fixed effects coefficients from the model, but want to add my own goodness-of-fit statistics under the table which is possible using the add.lines parameter in stargazer. Any help is appreciated.

Answer 1

Here is a way using kable and ehahelper::tidy.coxme:

test.qmd

---
format: pdf
---

```{r}
library(coxme)
library(broom)
library(kableExtra)

m1 <- coxme(Surv(time, status) ~ ph.ecog + age + (1|inst), lung)

  
tidy.coxme <- function(x, exponentiate = FALSE, conf.int = 0.95, ...){
    beta <- x$coefficients
    nvar <- length(beta)
    nfrail <- nrow(x$var) - nvar
    nn <- c("estimate", "exp()", "std.error", "statistic", "p.value")
    se <- sqrt(diag(as.matrix(x$var))[nfrail + 1:nvar])
    z <- qnorm((1 + conf.int)/2, 0, 1)
    ret <- data.frame(
      "term"      = names(beta),
      "estimate"  = beta,
      "std.error" = se,
      "statistic" = beta/se,
      "p.value"   = 1 - pchisq((beta/se)^2, 1),
      "conf.low"  =  beta - z * se,
      "conf.high" =  beta + z * se
    )
    if (exponentiate) {
      ret$estimate <- exp(ret$estimate)
      ret$conf.low <- exp(ret$conf.low)
      ret$conf.high <- exp(ret$conf.high)
    }
    rownames(ret) <- c(1:nrow(ret))
    ret
  }

tidy.coxme(m1) |> 
  kbl(format = "latex", digits = 2, booktabs = TRUE) |>
  kable_paper()
```