Standardisation des sorties de modeles

[DanChaltiel]

TODO

• Review des sorties déjà utilisées (Macros SAS ?)
• Utiliser broom::tidy() ou broom.helpers::tidy_plus_plus()
• Output plutôt flextable pour être consistant ? gtsummary propose des choses sympa mais je ne suis pas 100% convaincu (+ c'est galère à intégrer à Word avec officer)
• Ajouter des tests (usethis::use_test("models"), utiliser un dataset comme survival::lung ou autre)
• Check des hypothèses sous-jacentes ? Genre HHP et linéarité pour Cox ?

Modèles à implémenter (dans l'ordre de priorité)

• survival::coxph()
• glm() modèles logistique
• modèles mixtes -> plutôt mmrm que les vieux lme4 et nlme ?

Ressources

https://bsky.app/profile/yuzar-data-science.bsky.social/post/3lbqxjiaq622q

[DanChaltiel]

Exemple de fonction utiles :

library(ggsurvfit)
library(cli)
par = list()

fit_os  = survfit2(Surv(time_os, event_os) ~ ARM, data=db_prim)
fit_pfs = survfit2(Surv(time_pfs, event_pfs) ~ ARM, data=db_prim)

cox_os = coxph(Surv(time_os, event_os) ~ ARM + 
                 strata(strat_figo) + strata(strat_node) +
                 strata(strat_age)  + strata(strat_histo),
               data = db_prim)

report_coxph = function(fit, label){
  a = broom::tidy(fit, exponentiate=TRUE, conf.int=TRUE) %>% 
    filter(term=="armTTT") %>% 
    mutate(across(-term, ~round(.x, 2)))
  # a$estimate
  # fit$xlevels
  glue("The Hazard Ratio of treatment for {label} was {a$estimate} [95%CI {a$conf.low}; {a$conf.high}] (adjusted p-value: {a$p.value})")
}

report_survfit = function(fit, time, label, ...){
  diff = survdiff(formula(fit), data=get(fit$call$data))
  quants = 12*quantile(fit, prob=c(0.25, 0.50, 0.75), conf.int=FALSE) %>% round(1)
  #TODO ajouter médiane ?
  pars =
    tidy_survfit(fit, time=time, ...) %>% 
    select(strata, estimate, conf.low, conf.high) %>% 
    mutate(across(-strata, ~scales::percent(.x, 1)),
           label = glue("{estimate} [95%CI {conf.low}; {conf.high}] in the {strata} arm")) %>% 
    pull(label)
  
  format_inline("{label} was {pars} (Log-rank p-value = {round(diff$pvalue, 2)}).")
}

par$os = report_survfit(fit_os, time=2, label="2-year OS")
par$pfs = report_survfit(fit_pfs, time=2, label="2-year PFS")

par$cox_os =  report_coxph(cox_os, label="death")
par$cox_pfs = report_coxph(cox_pfs, label="RECIST progression or death")