longitudinal data table in gtsummary

library(data.table)
library(gtsummary)    
mm<-structure(list(time = c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 
    2, 3, 1, 2, 3, 1, 2, 3), meas1 = c(59, 0, 0, 109, 6.8, 0, 98, 
    11, 2.7, 121, 9.9, 0, 136, 63, 22.9, 130, 0, 0, 86, 23, 0), ID = c("001", 
    "001", "001", "002", "002", "002", "003", "003", "003", "004", 
    "004", "004", "005", "005", "005", "006", "006", "006", "007", 
    "007", "007"), gr = c("A", "A", "A", "B", "B", "B", "A", "A", 
    "A", "C", "C", "C", "C", "C", "C", "B", "B", "B", "B", "B", "B"
    )), row.names = c(NA, -21L), class = c("data.table", 
    "data.frame"))

Now, I´d like to generate a table with cells containing the median and IQR for the different levels of group (gr) at each of the different time points.

g <- function(x) {
  k <- as.list(quantile(x, (1 : 3) / 4, na.rm=TRUE))
  names(k) <- .q(Q1, Q2, Q3)
  k
}
mm2 <- mm[, g(meas1), by=.(gr, time)]
mm2[,iqr:=paste0(Q2,"(",Q1,"-",Q3,")")]
dcast(mm2,time~gr,value.var = "iqr")

Is there an easy way to do this in gtsummary? Preferably without inflating the N in the column header and extra bonus if the range could be with a newline centered on the median?

===================================================================

Edit after including the answer from @Daniel D. Sjoberg

tba <- unique(mm$time) %>% 
  sort() %>%
  purrr::map(
    ~ mm %>%
      dplyr::filter(time %in% .x) %>%
      tbl_summary(
        by = gr, 
        include = meas1,
        type = meas1 ~ "continuous2",
        statistic = meas1 ~ c("{median} ({p25}, {p75})", "{min}, {max}"),
        label = list(meas1 = paste("Measure 1 at Timepoint", .x))
      ) 
  ) %>%
  tbl_stack() 

fac <- as.list(as.data.frame(combn(levels(as.factor(mm$gr)),2)))
tbb <- lapply(fac,function(x){
  hh <- mm[gr%in%x]
  unique(hh$time) %>% 
    sort() %>% 
    purrr::map(
      ~hh %>% dplyr::filter(time %in% .x) %>% 
        tbl_summary( by = gr,
                     include = meas1,
                     type = meas1 ~ "continuous2",
                     statistic = meas1 ~ c("{median}","({p25}, {p75})"),
                     label = list(meas1 = paste("Measure 1 at Timepoint", .x))) %>%
        add_difference()) %>% 
    tbl_stack() %>% modify_column_hide(all_stat_cols())
})

tb_b <- tbb %>% tbl_merge(tab_spanner = unlist(lapply(fac,function(x) paste0("**",x[[1]]," vs. ",x[[2]],"**"))))  

tbl_merge(list(tba,tb_b)) 

tba results in desired output:

tb_b also results in desired output:

But- when merging these, the spanner from tb_b is lost identifying the contrasts. Any way of that I can retain the original spanners from tb_b when doing tbl_merge?

Also - the add_q helper is likely wrong when doing these stacked tables. Any way of fixing this?