ddsjoberg / gtsummary

Presentation-Ready Data Summary and Analytic Result Tables
http://www.danieldsjoberg.com/gtsummary
Other
1.04k stars 114 forks source link

Feature request: emmeans::emmeans() like contrasts, emmeans() output into {gtsummary} #1383

Closed yuryzablotski closed 1 year ago

yuryzablotski commented 1 year ago

Dear Daniel,

It would be amazing if it would be possible to include capabilities of “emmeans”, particularly, calculating “contrasts” into {gtsummary}?

What I meant by the contrasts is a simple difference between pairwise levels of categorical variables. Because if use tbl_regression(), we see the difference (in ORs) between the reference level and other levels, but not between levels II and III in the Grade predictor, or not the difference between Stage levels pairwisely (I never understood, why it’s missing in the base r function, like summary(model) ):

m <- glm(response ~ stage + age, trial, family = binomial)

tbl_regression(m, exponentiate = T, pvalue_fun = ~style_pvalue(.x, digits = 3))

emmeans::emmeans(m, pairwise ~ stage, type = "response", adjust = "none") %>% summary(infer = T)

$emmeans stage prob SE df asymp.LCL asymp.UCL null z.ratio p.value T1 0.361 0.0686 Inf 0.240 0.503 0.5 -1.923 0.0545 T2 0.247 0.0604 Inf 0.148 0.383 0.5 -3.429 0.0006 T3 0.347 0.0768 Inf 0.215 0.508 0.5 -1.864 0.0623 T4 0.309 0.0716 Inf 0.188 0.464 0.5 -2.395 0.0166

Confidence level used: 0.95 Intervals are back-transformed from the logit scale Tests are performed on the logit scale

$contrasts contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value T1 / T2 1.718 0.756 Inf 0.725 4.07 1 1.231 0.2185 T1 / T3 1.061 0.478 Inf 0.439 2.57 1 0.132 0.8948 T1 / T4 1.260 0.565 Inf 0.523 3.03 1 0.516 0.6062 T2 / T3 0.618 0.288 Inf 0.247 1.54 1 -1.033 0.3017 T2 / T4 0.733 0.343 Inf 0.293 1.83 1 -0.664 0.5066 T3 / T4 1.187 0.567 Inf 0.466 3.03 1 0.359 0.7192

Confidence level used: 0.95 Intervals are back-transformed from the log odds ratio scale Tests are performed on the log odds ratio scale

Besides, the first part of emmeans (see below) is also useful, and it would be really grate to be able to beautifully gt-summaries the results of both, $emmeans and $contrasts

The final tbl_regression(m) should look like this:

Screenshot 2022-11-09 at 18 04 08
ddsjoberg commented 1 year ago

Can you please update your post with an example of what you want the final table to look like?

larmarange commented 1 year ago

I need to explore it. An approach could be a custom tidier handled adequately by broom.helpers

Note: currently traveling for work

yuryzablotski commented 1 year ago

Well, I did not find a better and simpler way to get to contrasts as with emmeans. And since I realized I use emmens every time I model, I wander, why it (post-hoc pairwise comparisons) isn't possible with classic model-results functions, like summary()? Anyway, it would be soooo amazing, if we would be able to use emmeans-arguments, like "type = "response" ", or "adjust = "none" " in {gtsummary}

larmarange commented 1 year ago

I have drafted quickly a possible solution through a new function tidy_add_pairwise_contrasts() in broom.helpers.

tidy_plus_plus() gains accordingly additional arguments

library(gtsummary)
library(broom.helpers)
#> 
#> Attachement du package : 'broom.helpers'
#> Les objets suivants sont masqués depuis 'package:gtsummary':
#> 
#>     all_continuous, all_contrasts

mod <- glm(response ~ age + trt + stage + grade, data = trial, family = binomial)

mod %>%
  tidy_plus_plus(
    add_pairwise_contrasts = TRUE, 
    exponentiate = TRUE
  ) %>%
  knitr::kable()
term variable var_label var_class var_type var_nlevels contrasts contrasts_type reference_row label n_obs n_event estimate std.error statistic p.value conf.low conf.high
age age Age numeric continuous NA NA NA NA Age 183 58 1.0195893 0.0115076 1.6858397 0.0918267 0.9971774 1.043416
trtDrug A trt Chemotherapy Treatment character dichotomous 2 contr.treatment treatment TRUE Drug A 89 27 1.0000000 NA NA NA NA NA
trtDrug B trt Chemotherapy Treatment character dichotomous 2 contr.treatment treatment FALSE Drug B 94 31 1.1500613 0.3232223 0.4325668 0.6653295 0.6106565 2.176274
Drug B / Drug A trt Chemotherapy Treatment character dichotomous 2 pairwise pairwise NA Drug B / Drug A NA NA 1.1500613 0.3717255 0.4325668 0.6653295 0.6103707 2.166947
stageT1 stage T Stage factor categorical 4 contr.treatment treatment TRUE T1 50 18 1.0000000 NA NA NA NA NA
stageT2 stage T Stage factor categorical 4 contr.treatment treatment FALSE T2 51 13 0.5628967 0.4437259 -1.2950770 0.1952937 0.2318960 1.333426
stageT3 stage T Stage factor categorical 4 contr.treatment treatment FALSE T3 39 14 0.9049542 0.4573824 -0.2183534 0.8271538 0.3651204 2.213548
stageT4 stage T Stage factor categorical 4 contr.treatment treatment FALSE T4 43 13 0.7633870 0.4537622 -0.5950036 0.5518411 0.3093673 1.850377
T2 / T1 stage T Stage factor categorical 4 pairwise pairwise NA T2 / T1 NA NA 0.5628967 0.2497718 -1.2950770 0.5660480 0.1800348 1.759952
T3 / T1 stage T Stage factor categorical 4 pairwise pairwise NA T3 / T1 NA NA 0.9049542 0.4139101 -0.2183534 0.9963269 0.2794586 2.930459
T3 / T2 stage T Stage factor categorical 4 pairwise pairwise NA T3 / T2 NA NA 1.6076737 0.7542474 1.0120081 0.7424042 0.4816755 5.365884
T4 / T1 stage T Stage factor categorical 4 pairwise pairwise NA T4 / T1 NA NA 0.7633870 0.3463962 -0.5950036 0.9336213 0.2379440 2.449146
T4 / T2 stage T Stage factor categorical 4 pairwise pairwise NA T4 / T2 NA NA 1.3561761 0.6346940 0.6509987 0.9152375 0.4075285 4.513093
T4 / T3 stage T Stage factor categorical 4 pairwise pairwise NA T4 / T3 NA NA 0.8435643 0.4045959 -0.3546908 0.9847078 0.2460316 2.892314
gradeI grade Grade factor categorical 3 contr.treatment treatment TRUE I 65 21 1.0000000 NA NA NA NA NA
gradeII grade Grade factor categorical 3 contr.treatment treatment FALSE II 58 17 0.8395776 0.4027265 -0.4341815 0.6641566 0.3779389 1.845925
gradeIII grade Grade factor categorical 3 contr.treatment treatment FALSE III 60 20 1.0384054 0.3892503 0.0968175 0.9228713 0.4824457 2.232749
II / I grade Grade factor categorical 3 pairwise pairwise NA II / I NA NA 0.8395776 0.3381201 -0.4341815 0.9013229 0.3266955 2.157637
III / I grade Grade factor categorical 3 pairwise pairwise NA III / I NA NA 1.0384054 0.4041996 0.0968175 0.9948455 0.4170288 2.585638
III / II grade Grade factor categorical 3 pairwise pairwise NA III / II NA NA 1.2368188 0.5043448 0.5212242 0.8609693 0.4756077 3.216350

Created on 2022-11-12 with reprex v2.0.2

larmarange commented 1 year ago

You can play with it here: https://github.com/larmarange/broom.helpers/pull/192

My feeling is that it could be displayed as usual by tbl_regression(). However, it would require additional arguments to be added to tbl_regression().

A more generic approach could be to have a way to pass additional arguments to tidy_plus_plus(). It could be simply by using ... approach. Otherwise, by providing a list with something like tidy_plus_plus_args

yuryzablotski commented 1 year ago

I can't get it to work. Unfortunately. It does not produce contrasts. Do I need a developer version of gtsummary or broom.helpers? Because with usual versions it does not produce contrasts? Something is missing. Do you have any suggeestion what it might be?

larmarange commented 1 year ago

So far, it is only implemented in a PR. So to install that version of broom.helpers, you need to use remotes::install_github("larmarange/broom.helpers#192")

Then, you will be able to test the function broom.helpers::tidy_plus_plus() with is used internally by tbl_regression(). So you could see how pairwise contrasts could be computed.

However at that stage, it would not be possible to test with tbl_regression(), as a change in gtsummary would also be required, to pass new arguments to tidy_plus_plus()

larmarange commented 1 year ago

As a proof of concept, I have prepared a PR in gtsummary allowing to pass additional arguments to tidy_plus_plus() cf. https://github.com/ddsjoberg/gtsummary/pull/1387

If you install also this PR with remotes::install_github("ddsjoberg/gtsummary#1387"), both changes are working together.

library(gtsummary)
#> #BlackLivesMatter
mod <- glm(response ~ age + trt + stage + grade, data = trial, family = binomial)
mod |> 
  tbl_regression(exponentiate = TRUE, add_pairwise_contrasts = TRUE) |> 
  as_kable()
Characteristic OR 95% CI p-value
Age 1.02 1.00, 1.04 0.092
Chemotherapy Treatment
Drug B / Drug A 1.15 0.61, 2.17 0.7
T Stage
T2 / T1 0.56 0.18, 1.76 0.6
T3 / T1 0.90 0.28, 2.93 >0.9
T3 / T2 1.61 0.48, 5.37 0.7
T4 / T1 0.76 0.24, 2.45 >0.9
T4 / T2 1.36 0.41, 4.51 >0.9
T4 / T3 0.84 0.25, 2.89 >0.9
Grade
II / I 0.84 0.33, 2.16 >0.9
III / I 1.04 0.42, 2.59 >0.9
III / II 1.24 0.48, 3.22 0.9

Created on 2022-11-19 with reprex v2.0.2

yuryzablotski commented 1 year ago

Now it has worked! That's amazing! 🤩 Contrasts are working with gaussian family of glm also! I love all of gtsummary functions, but the possibility to get contrasts puts it on another level! Love it! Will use it daily! Thanks Joseph and Daniel!

Just for your information:

Wishful thinking of my is - it would be really useful, if we could use the functionality of {emmeans} in {gtsummary} somehow, particularly, to determine which contrasts we'll use after an interaction:

emmeans::emmeans(mod, pairwise ~ stage | grade, type = "response", adjust = "none") %>% summary(infer = T)

I'll use it only without interactions so far, that is already cool enough!

larmarange commented 1 year ago

For interactions, the warning is coming from emmeans.

tbl_uvregression() is not covered but it could be complex to update it. I let @ddsjoberg comment on that. I do not know if it should be relevant/possible to pass additional args.

For your last point, you can pass additional args to emmeans:

library(gtsummary)
#> #BlackLivesMatter
library(broom.helpers)
#> 
#> Attachement du package : 'broom.helpers'
#> Les objets suivants sont masqués depuis 'package:gtsummary':
#> 
#>     all_continuous, all_contrasts

mod <- glm(
  response ~ stage + grade + trt,
  gtsummary::trial,
  family = binomial,
)

mod |> 
  tbl_regression(
    exponentiate = TRUE,
    add_pairwise_contrasts = TRUE,
    pairwise_variables = c("stage", "grade"),
    emmeans_args = list(by = "trt")
  ) |> 
  as_kable()
Characteristic OR 95% CI p-value
T Stage
T2 / T1 | Drug A 0.61 0.20, 1.86 0.7
T3 / T1 | Drug A 1.10 0.35, 3.45 >0.9
T3 / T2 | Drug A 1.82 0.56, 5.94 0.6
T4 / T1 | Drug A 0.81 0.27, 2.42 >0.9
T4 / T2 | Drug A 1.33 0.42, 4.18 >0.9
T4 / T3 | Drug A 0.73 0.23, 2.34 0.9
T2 / T1 | Drug B 0.61 0.20, 1.86 0.7
T3 / T1 | Drug B 1.10 0.35, 3.45 >0.9
T3 / T2 | Drug B 1.82 0.56, 5.94 0.6
T4 / T1 | Drug B 0.81 0.27, 2.42 >0.9
T4 / T2 | Drug B 1.33 0.42, 4.18 >0.9
T4 / T3 | Drug B 0.73 0.23, 2.34 0.9
Grade
II / I | Drug A 0.96 0.39, 2.37 >0.9
III / I | Drug A 1.13 0.47, 2.76 >0.9
III / II | Drug A 1.19 0.48, 2.95 0.9
II / I | Drug B 0.96 0.39, 2.37 >0.9
III / I | Drug B 1.13 0.47, 2.76 >0.9
III / II | Drug B 1.19 0.48, 2.95 0.9
Chemotherapy Treatment
Drug A
Drug B 1.24 0.67, 2.30 0.5

Created on 2022-11-25 with reprex v2.0.2

However, I don't think we could cover all possibilities of emmeans.

yuryzablotski commented 1 year ago

Sure, I understand! And I am already very grateful for what you did so far!

The example above though me an error-message:

Error in dplyr::bind_rows(): ! Can't combine ..1$estimate and ..2$estimate <factor>. Run rlang::last_error() to see where the error occurred.

larmarange commented 1 year ago

You should update broom.helpers with the last version of the PR.

Joseph Larmarange

Le ven. 25 nov. 2022 à 16:11, yuryzablotski @.***> a écrit :

Sure, I understand! And I am already very grateful for what you did so far!

The example above though me an error-message:

Error in dplyr::bind_rows(): ! Can't combine ..1$estimate and ..2$estimate . Run rlang::last_error() to see where the error occurred.

— Reply to this email directly, view it on GitHub https://github.com/ddsjoberg/gtsummary/issues/1383#issuecomment-1327614723, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAHL5I6L5LP27RNKSTTHVZDWKDJKXANCNFSM6AAAAAAR3RGMUU . You are receiving this because you commented.Message ID: @.***>