rstargazercox-regressionxtabletexreg

How to export coxme regression output to Latex?


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.


Solution

  • 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()
    ```
    

    enter image description here