Combining mstate results using gtsummary

Question

I am trying to nicely collate the output from mstate using gtsummary. The following code does the job, but is quite unwieldly. I'm sure there is an elegant way to do this -- but this is my current best attempt.

Is there some way to nicely iterate through the list of tbl_regression objects to change the variable and labels so that tbl_merge recognizes that these the same variables? Is there some way to change the name of the variables before passing the model to tbl_regression? Is there a way to efficiently keep track of which transition the parameters originated from?

Thank you so much!

library(mstate)
#> Loading required package: survival
library(dplyr)
#> Warning: package 'dplyr' was built under R version 4.0.5
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(gtsummary)
#> Warning: package 'gtsummary' was built under R version 4.0.5
#> #Uighur
library(stringr)

data("ebmt4")
ebmt <- ebmt4

## Data prep
# Create transition matrix
tmat <- transMat(x = list(c(2, 3), c(3), c()),
                 names = c("Tx", "Rec", "Death"))
# tmat

# Create msdata object (wide -> long)
msebmt <- msprep(data = ebmt, trans = tmat, 
                 time = c(NA, "rec", "srv"), 
                 status = c(NA, "rec.s","srv.s"), 
                 keep = c("match", "proph", "year"))

covs <- c("match", "proph", "year")
msebmt.expand <- expand.covs(msebmt, covs, longnames = FALSE)

# convert time to years
msebmt.expand[, c("Tstart", "Tstop", "time")] <- msebmt.expand[, c("Tstart", "Tstop", "time")]/365.25

## mstate model
cfull <- coxph(Surv(Tstart, Tstop, status) ~ 
match.2 + match.3 + proph.1 + proph.2 + proph.3 + year1.1 + 
year1.3 + strata(trans), data = msebmt.expand, method = "breslow")
cfull
#> Call:
#> coxph(formula = Surv(Tstart, Tstop, status) ~ match.2 + match.3 + 
#>     proph.1 + proph.2 + proph.3 + year1.1 + year1.3 + strata(trans), 
#>     data = msebmt.expand, method = "breslow")
#> 
#>             coef exp(coef) se(coef)      z        p
#> match.2 -0.01239   0.98769  0.10627 -0.117  0.90722
#> match.3  0.31782   1.37413  0.11716  2.713  0.00667
#> proph.1 -0.38050   0.68352  0.07169 -5.308 1.11e-07
#> proph.2  0.24675   1.27986  0.09823  2.512  0.01201
#> proph.3  0.16939   1.18458  0.12474  1.358  0.17449
#> year1.1  0.05683   1.05847  0.05840  0.973  0.33052
#> year1.3 -0.13797   0.87112  0.10825 -1.275  0.20246
#> 
#> Likelihood ratio test=48.79  on 7 df, p=2.498e-08
#> n= 5776, number of events= 2056

## gtsummary attempt to prettify

# make a large tbl_regression
cfull.tbl <- cfull %>% tbl_regression(exponentiate = TRUE)
# split the large tbl_regression into individual rows
names.cfull.tbl <- cfull.tbl$table_body$variable
cfull.tbl.split <- tbl_split(cfull.tbl, c(all_of(names.cfull.tbl)))

# change names so that they no longer indicate which transition state they originated from and can thus be combined (this code thanks to Arun on R4DS slack)
cfull.tbl.split.nonums <- purrr::map(
  cfull.tbl.split, 
  function(x) {
    x$table_body$label <- str_sub(x$table_body$label, start = 1L, end = -3L)
    x$table_body$variable <- str_sub(x$table_body$variable, start = 1L, end = -3L)
    return(x)
  }
)

# combine transition states
cfull.tbl.trans1 <- 
  tbl_stack(
    tbls = list(
      cfull.tbl.split.nonums[[3]],
      cfull.tbl.split.nonums[[6]])
  )

cfull.tbl.trans2 <- 
  tbl_stack(
    tbls = list(
      cfull.tbl.split.nonums[[1]],
      cfull.tbl.split.nonums[[4]])
  )

cfull.tbl.trans3 <- 
  tbl_stack(
    tbls = list(
      cfull.tbl.split.nonums[[2]],
      cfull.tbl.split.nonums[[5]],
      cfull.tbl.split.nonums[[7]])
  )
# combine into one object
  tbl_merge(
    tbls = list(
      cfull.tbl.trans1,
      cfull.tbl.trans2,
      cfull.tbl.trans3),
    tab_spanner = c("tx -> recov",
                    "tx -> death",
                    "recov -> death")
  )